Ê7/25 ½¤ÀµÈÇ¡Ë

統計科学同演習第 14 回（7/25 修正版）
清智也 + TA’s
2014 年 7 月 18 日（金）
線形モデルのつづき
モデルの検証
シミュレーション
授業のまとめ
1 / 24
線形モデル（復習）
I
線形モデル
yi = a +
p
∑
bj xij + i ,
i = 1, . . . , n
j=1
I
ベクトルと行列による書き換え
y = Xβ + I
最小二乗法
Minimize ky − Xβk2
β
I
残差：最小二乗法で求めた解 βˆ に対して，ˆ = y − Xβˆ のこと．
2 / 24
補足：残差を検討する意味
（先週の演習のとき，アナウンスのページにも書きました．
）
I
変量 yi , xi1 , xi2 があるとする．
I
例：yi が賃料，xi1 が徒歩時間，xi2 が建築年．
I
最初に徒歩時間 xi1 についてあてはめたとする．
結果を yi = ˆ
a + bˆ1 xi1 + ˆi とする．
I
I
ˆ i1 ) と変量 xi2 の散布図を描いてみ
もし残差 ˆi = yi − (ˆ
a + bx
て，これらの間に線形関係（相関関係）が見られれば，
yi = a + b1 xi1 + b2 xi2 + i
というモデルを検討する余地がある．
I
このとき，元のモデルの回帰係数 a, b と新しいモデルの a,
b1 の推定値は一般に異なる．
3 / 24
補足：交互作用
I
交互作用項を含む線形モデル
yi = a + b1 xi1 + b2 xi2 + b12 xi1 xi2 + i
I
傾きが他の変数に依存するという解釈ができる．
yi = a + b1 xi1 + (b2 + b12 xi1 )xi2 + i
I
例：y は筋肉，x1 は運動量，x2 はプロテイン摂取量．
4 / 24
補足：対比
I
I
I
因子変量 x が K 個の水準 l1 , . . . , lK を持つとする．
例：l1 =”日吉”, l2 =”武蔵小杉”, l3 =”田園調布”．
式 y ~ -1+x は次のモデルに対応する
yi = b1 1{xi =l1 } + · · · + bK 1{xi =lK } + i
I
回帰係数は b1 , . . . , bK の K 個．
式 y ~ x は次のモデルに対応する
yi = a + b2 1{xi =l2 } + · · · + bK 1{xi =lK } + i
I
I
回帰係数は a, b2 , . . . , bK の K 個
以上の 2 つはモデルとしては全く同等である．しかし，
「回帰
係数が 0 か否か」などを考えるとき，その意味は異なる．
y ~ -1+x+is.na(bus)*walk は次のモデルに対応．
yi =
K
X
bk 1{xi =lk } +c1 (is.na(bus))i +c2 (walk)i +c3 (is.na(bus))i ×(walk)i +i
k=1
5 / 24
多項式や基底関数を用いた回帰
I
線形回帰モデルは，パラメータに関して線形であればよい．
I
例えば次の多項式モデルは線形回帰モデルの一例である：
yi = a + b1 xi + b2 xi2 + b3 xi3 + i
I
多項式以外の基底関数：局所的な変化を表現
I
I
I
I
スプライン（滑らかな区分多項式系）
ウェーブレット（多重解像度を持つ直交関数系）
2
2
動径基底関数（ガウス型の基底 {e −kx−µj k /(2σ ) }j など）
あまり柔軟に表現しすぎると過適合（over fitting）を起こす．
そのため，ペナルティ項を付けて対処する．
I
例えばこんな感じ → 基底関数を v1 (x), . . . , vk (x) として
{
}
∫ u
Minimize
kyi − φ(xi )k2 + λ
φ00 (x)2 dx
φ∈span(v1 ,...,vk )
l
λ > 0 は別途決定する（後述の CV など）．
6 / 24
正規線形モデル
ここまでは特に確率変数を明示しなかった．
I
正規線形モデル（単に線形モデルと言うことも多い）
yi ∼ N(µi , σ ),
2
µi = a +
p
∑
bj xij
j=1
I
復習：平均 µ，分散 σ 2 の正規分布 N(µ, σ 2 ) の確率密度関数は
I
1
2
2
e −(y −µ) /(2σ )
f (y ; µ, σ 2 ) = √
2
2πσ
∏n
y1 , . . . , yn の同時密度関数は i=1 f (yi ; µi , σ 2 ) で与えられる．
7 / 24
最尤法
I
一般に，同時密度関数が最大となるようなパラメータを求め
る方法を最尤法 (maximum likelihood method) という．
I
演習：正規線形モデルに対する最尤法から自然に最小二乗法
が導かれることを示せ．
（ただし σ 2 は定数としてよい）
8 / 24
モデルの検証
目的
I
説明変数を選択する．
I
過適合を防ぐ．
方法
I
検定
I
情報量規準
I
ブートストラップ
I
クロスバリデーション (CV)
I
正則化法（ペナルティ項を付ける）
以下，正規線形モデルに限って，これらの方法の概観を述べる．
9 / 24
検定（統計的仮説検定）
I
I
検定は，背理法の考え方に近い．
検定の手順：モデル yi = a + b1 xi1 + b2 xi2 + i の場合．
1. 母集団において b1 = 0 であると仮定する（背理法の仮定）．
2. この仮定のもとで，推定された回帰係数 bˆ1 の起こり具合
（p 値）を計算する．p 値の正確な定義は省略する．
3. p 値が小さければ（たとえば 0.05 を下回れば）矛盾と考え，
変数 x1 は説明力を持つと結論する．このとき「有意」という．
4. x2 や，x1 , x2 のペアについても同様のことを行う．
I
尤度比検定：尤度関数（同時密度関数）の最大値を比較する．
I
分散分析 (analysis of variance; ANOVA) ：正規線形モデルに
対する尤度比検定．標本分散に関連する平方和 ky − y¯ 1k2 の
分解によって解釈できるため，この名が付いている．
I
数理統計学第２やデータサンプリングの講義でも扱われる予定．
10 / 24
例：サッカーのシュート数
注意：このデータにはポアソン回帰モデル（後述）などを用いた方が適切．
> lm1 = lm(shoot ~ playing * position, data=soccer)
> lm1$coef
(Intercept)
playing
positionFW
-1.365374152
0.007257349
2.303964300
positionMF playing:positionFW playing:positionMF
6.464985065
0.015970817
0.004702370
> anova(lm1)
Analysis of Variance Table
Response: shoot
Df Sum Sq Mean Sq F value
Pr(>F)
playing
1 15441.8 15441.8 117.838 < 2.2e-16 ***
position
2 16739.7 8369.9 63.871 < 2.2e-16 ***
playing:position
2 4564.8 2282.4 17.417 9.814e-08 ***
Residuals
214 28043.0
131.0
--Signif. codes: 0 ‘ *** ’ 0.001 ‘ ** ’ 0.01 ‘ * ’ 0.05 ‘ . ’ 0.1 ‘
I
Pr(>F) と書かれている部分が p 値．
I
モデル：出場時間 x ，ポジション z ，シュート数 y として
’1
yi = a + b1 xi + b2 1{zi =FW } + b3 1{zi =MF } + (b4 xi 1{zi =FW } + b5 xi 1{zi =MF } ) + i
11 / 24
情報量規準
I
赤池情報量規準 (AIC)
AIC = −2
n
∑
log f (yi |ˆ
µi , σ
ˆ 2 ) + 2(p + 1)
i=1
ここで p は説明変数の個数．
I
AIC が小さくなるように説明変数を選ぶ．
> lm1 = lm(shoot ~ playing * position, data=soccer)
> AIC(lm1)
[1] 1704.864
> lm0 = lm(shoot ~ playing + position, data=soccer)
> AIC(lm0)
[1] 1734.042
12 / 24
ブートストラップ
I
第 12 回の講義で補足した手法．
I
「データをリサンプリングして推定する」という操作を繰り
返すことにより，推定値の標準誤差が求まる．
I
正規線形モデルに対しては，解析的に標準誤差が求まるので
あまり使わない（と思う）．
I
クラスタリングや主成分分析のような手法にも適用できる．
13 / 24
クロスバリデーション
I
データを「訓練データ」と「テストデータ」に分ける．
I
訓練データだけで得られた回帰係数を使って，テストデータ
のあてはまり具合（予測性能）を評価する．
I
この操作を交差的に行う．例えば，各 i ∈ {1, . . . , n} に対し
て，i 番目の観測値を除いたデータを訓練データとして回帰
係数を求め，それを使って (xi , yi ) の予測性能をはかる．
1∑
(yi − yˆi,−i )2 ,
n
n
CV =
yˆi,−i = ˆa−i + bˆ−i xi ,
i=1
(ˆa−i , bˆ−i ) = argmin
a,b
∑
(yj − (a + bxj ))2 .
j6=i
※簡単のため σ 2 の推定は考慮に入れていない．
14 / 24
例
> demo.cv
function(){
soccer = read.csv("soccer.csv")
n = nrow(soccer)
er = matrix(NA, n, 2)
for(i in 1:n){
lm0 = lm(shoot ~ playing + position, data=soccer[-i,])
er[i,1] = soccer$shoot[i] - predict(lm0, soccer[i,])
lm1 = lm(shoot ~ playing * position, data=soccer[-i,])
er[i,2] = soccer$shoot[i] - predict(lm1, soccer[i,])
}
print(mean(er[,1]^2))
print(mean(er[,2]^2))
}
> demo.cv()
[1] 153.4587
[1] 133.4336
15 / 24
正則化法
I
最小二乗法にペナルティ項を付けることで過適合を防ぐ．
{
}
Minimize ky − Xβk2 + λ pen(β)
β
I
ペナルティの例
I
I
∑p
pen(β) = j=1 βj2 (ridge)
∑p
pen(β) = j=1 |βj | (lasso)
16 / 24
シミュレーション
少し話題を変えて
Example (待ち行列, queue (6= matrix))
I
駅に 4 つの券売機がある．一気にまとまった人が切符を買い
にきて，自分は 9 番目に到着した．
I
1 人が切符を買うのにかかる時間 [秒] は次の確率密度関数に
従うと仮定する：
 1
 30 if 20 < x ≤ 40,
1
f (x) =
if 40 < x ≤ 60,
 60
0 otherwise.
I
次の各状況について，待ち時間の平均と標準偏差を求めてみ
よう．
1. 客が券売機ごとに分かれて並ぶ（自分の前には 2 人）．
2. 顧客が 1 列に並び，空いた券売機に順番に移動する（自分の
前には 8 人）．
17 / 24
シミュレーション
105 回の乱数実験．関数の中身は各自確認のこと．
> source("e14.R")
> set.seed(20140718) # 乱数の種 (seed)
> queue.par() # 「並列」の場合
$mean
[1] 73.3804
$sd
[1] 15.60276
> queue.ser()
$mean
[1] 58.80944
# 「直列」の場合
$sd
[1] 7.989037
18 / 24
解析的に計算する場合
I
「並列」の場合の期待値，分散，標準偏差
E [X + Y ] = E [X ] + E [Y ] (期待値の線形性)
(∫ 40
)
∫ 60
x
x
= 2E [X ] = 2
dx +
dx
20 30
40 60
220
=
= 73.33 · · ·
3
V [X + Y ] = V [X ] + V [Y ] (独立性から)
)
(
= 2V [X ] = 2 E [X 2 ] − E [X ]2
(∫ 40 2
)
∫ 60 2
x
x
2
=2
dx +
dx − (110/3)
20 30
40 60
2200
=
= 244.4 · · ·
√
√9
V [X + Y ] = 2200/9 = 15.63 · · ·
I
注意：公式 V [X ] = E [X 2 ] − E [X ]2 は，数値計算上は使わない方が
よい（大きな正の数どうしの引き算は「桁落ち」が生ずる）．
19 / 24
解析的に計算する場合
I
「直列」の場合…
I
8 人が切符を買う時間をそれぞれ X1 , . . . , X8 とおく．
待ち時間を表す確率変数 W は次のように定まる：
I
1. P1 = X1 , P2 = X2 , P3 = X3 , P4 = X4 とおく．
2. for i ∈ {5, 6, 7, 8}
I
Pj = min(P1 , P2 , P3 , P4 ) なる j に対して Pj ← Pj + Xi とおく．
3. W = min(P1 , P2 , P3 , P4 ).
I
原理的には，場合分けによって W の期待値，分散，標準偏
差を解析的に計算できる．（でも計算してません）
I
客の到着時刻や処理時間に確率分布を仮定して，待ち時間の
分布・特性を調べる理論は一般に待ち行列理論 (queueing
theory) と呼ばれる．
20 / 24
この授業で扱ったこと
I
扱ったデータ：地震，株価，遺伝子，気象，指の長さ，家賃．
I
前処理：経度・緯度の変換，欠損値の補間，文字列の置換
など．
I
可視化：散布図，Q-Q プロット，経験分布，3D プロット，箱
形図，matplot，dotplot，biplot など．
I
分析手法：クラスター分析，主成分分析，偏相関，線形モデ
ル（回帰分析）など．
I
R の使い方：クラス，オブジェクト，データフレーム，関数
など．
データ分析は一辺倒ではないので，例えば家賃データに主成分分
析を用いることもあり得るし，気象データにクラスタリングを使
うこともあり得る．いろいろ試してみるとよいと思う．
21 / 24
この授業で触れなかったこと
分析手法のみ列挙してみる
I 一般化線形モデル (generalized linear model; GLM)
I
I
I
ロジスティック回帰
ポアソン回帰 yi ∼ Poisson(λi ), log λi = a + bxi
判別分析 (discrimination analysis)
I
I
線形判別分析 (LDA)
サポートベクターマシン (SVM)
I
多次元尺度法 (multi-dimensional scaling; MDS)
I
数量化 III 類 (corresponding analysis)
平滑化 (smoothing)
I
I
I
I
スプライン平滑化
カーネル平滑化
他にもいろいろ
22 / 24
期末試験について
I
試験日時は 7 月 25 日（金）3 限です．
I
詳しくは学事のページを参照してください．
I
初回に説明した通り，レポートの合計点が合格点に達してい
れば期末試験は受けなくても合格です（A, B, C の基準につ
いても初回の資料を参照してください）．
I
今日返却するレポートまでの合計点（第 2 回大レポは除く）
はお伝えします．
I
第 2 回大レポと，今日の演習の点数は，採点次第，各自の
keio.jp のアカウントにお知らせする予定です（試験前日の夜
くらいかも）．
I
FD アンケートにもご協力ください．
https://fd-enquete.st.keio.ac.jp/sfc-sfs-st/index.cgi
23 / 24
復習事項
I
地震のマグニチュードの従う分布としてどんな分布が考えられる
か答えなさい．
I
対数収益率の定義を述べなさい．
I
RNA (DNA) データの演習ではどんな距離行列を用いたか答えな
さい．
I
主成分分析の概要を述べなさい．
I
（標本）相関係数の定義を述べなさい．
I
累積分布関数とは何か説明しなさい．
I
R におけるデータフレームと行列の違いを一つ述べなさい．
I
R 関数における引数の省略時既定値の役割を説明しなさい．
I
授業で説明した，個体空間と変量空間とは何か答えなさい．
I
固定欄形式と自由欄形式の違いを述べなさい．
24 / 24

Download Report