Åý·×²Ê³ØƱ±é½¬ Âè6²ó ¡Ê½¤ÀµÈÇ¡Ë

統計科学同演習 第 6 回 (修正版)
清 智也 + TA’s
2014 年 5 月 16 日(金)
線形モデル
最小二乗法
正規線形モデル
独立性
1 / 13
線形モデル
I
変量 Y を変量 X で「説明」したい.
Example (cars データ)
60
40
20
0
cars$dist
80
100
120
> plot(cars$speed, cars$dist)
> lm.cars = lm(dist ~ speed, data=cars)
> abline(lm.cars)
5
10
15
20
25
cars$speed
2 / 13
lm クラス
lm = linear model = 線形モデル(後述)
> lm.cars = lm(dist ~ speed, data=cars)
> attributes(lm.cars)
$names
[1] "coefficients" "residuals"
"effects"
[4] "rank"
"fitted.values" "assign"
[7] "qr"
"df.residual"
"xlevels"
[10] "call"
"terms"
"model"
$class
[1] "lm"
> lm.cars$coef # 切片と傾き
(Intercept)
speed
-17.579095
3.932409
切片は -17.579095,傾きは 3.932409
→ 求め方の理論は後で.
3 / 13
線形モデル
I
モデル = 模型.
I
一般に,2 組の確率変数 (X , Y ) に対して,X = x が与えられ
た下での Y の条件付き分布のモデルを回帰モデルという.
I
X を説明変量という.説明変量は複数個(多次元)であって
もよいが,今日は簡単のため 1 個の場合を考える.
I
加法的な誤差のあるモデル:
Y = f (x) + ε
I
線形モデル (線形回帰モデル)
Y = a + bx + ε
誤差 ε の平均は E [ε] = 0,分散は Var [ε] = σ 2 .
I
直線の式 y = a + bx を回帰直線という.
4 / 13
最小二乗法
I
a, b, σ 2 は未知なので,データ {(xi , yi )}ni=1 から推定する必要
がある.
I
基本は最小二乗法
Minimize
a,b
n
∑
{yi − (a + bxi )}2
i=1
演習
データが中心化されているとき,すなわち
n
∑
i=1
xi = 0,
n
∑
yi = 0
i=1
のとき,上の最小化問題の解を求めよ.
5 / 13
最小二乗法
中心化されているとは限らない場合の最小二乗法の解
∑n
(x − x¯)(yi − y¯ )
∑n i
a = y¯ − b¯
x , b = i=1
¯)2
i=1 (xi − x
ただし
1∑
xi ,
x¯ =
n
n
i=1
I
1∑
y¯ =
yi .
n
n
i=1
性質:y¯ = a + b¯
x .標本平均のペア (¯
x , y¯ ) は回帰直線上に
ある.
6 / 13
最小二乗法
より抽象的には,
I Rn において,ベクトル y = (y1 , . . . , yn )T を 2 つのベクトル
1 = (1, . . . , 1)T , x = (x1 , . . . , xn )T の張る線形空間に直交射
影している.
I 直交射影
Minimize ky − Xβk2
β
ただし
X = (1, x) ∈ R
I
I
n×2
,
( )
a
∈ R2 .
β=
b
行列の計算から,解は
β = (XT X)−1 XT y.
(1)
Xβ = X(XT X)−1 XT y.
(2)
説明変量が複数個の場合(後日説明)も同じ式になる.
演習
式 (1), (2) を確かめよう!
7 / 13
最小二乗法の背景:正規線形モデルのあてはめ
I
(正規)線形モデル
Yi = a + bxi + εi ,
i = 1, . . . , n.
ε1 , . . . , εn ∼ N(0, σ 2 ) 独立
I
Yi の確率密度関数
(
)
(yi − (a + bxi ))2
f (yi |xi ) = √
exp −
2σ 2
2πσ 2
1
I
同時確率密度
n
∏
i=1
f (yi |xi ) =
( ∑n
2)
1
i=1 (yi − (a + bxi ))
.
exp
−
2σ 2
(2πσ 2 )n/2
これを最大化する (a, b, σ 2 ) を最尤推定量といい,a, b は最小
二乗法の解となる.
8 / 13
なぜ正規分布を使うのか?
(
)
(x − µ)2
f (x) = √
exp −
,
2σ 2
2πσ 2
1
I
正規分布:誤差の分布
I
根拠:中心極限定理
x ∈ R.
Theorem (中心極限定理)
Z1 , . . . , Zn が独立同一分布に従い,E [Zi ] = 0, σ 2 = Var [Zi ] < ∞
ならば,
Z1 + · · · + Zn
√
n
は正規分布 N(0, σ 2 ) に収束する.
収束(位相)の定義や証明についてはここでは触れない.
9 / 13
あてはめ結果
I
ˆ と書く.
以下,最小二乗法で求めた (a, b) を (ˆ
a, b)
I
残差 (residuals)
ˆ i ).
εˆi = yi − (ˆa + bx
I
Rn で考えると,
εˆ ⊥ 1, x
I
正規線形モデルの仮定の下では,εˆi も正規分布に従う.ただ
し独立にはならない.
I
残差の標準化については来週補足説明する予定.
I
残差の正規性をチェックすることにより,モデルの妥当性を
検討することができる.
10 / 13
独立性 (追記)
確率変数 X , Y が独立
⇐⇒ 任意の集合 A, B ⊂ R に対して
P(X ∈ A, Y ∈ B) = P(X ∈ A)P(Y ∈ B).
⇐⇒ 任意の集合 A, B ⊂ R に対して
P(Y ∈ B | X ∈ A) = P(Y ∈ B).
視覚的に独立性を検証するための一つの方法:
I
集合 A として,X の四分位点で区切った 4 つの区間を用いる.
I
X ∈ A のもとでの条件付き分布を boxplot で描く.
I
boxplot がほぼ同じだったら独立,違ったら独立でない,と
判断(次のスライド)
11 / 13
独立性 (追記)
dist
0
20
40
60
80
100
120
> attach(cars)
> plot(dist ~ cut(speed, quantile(speed)))
(4,12]
(12,15]
(15,19]
(19,25]
cut(speed, quantile(speed))
位置(や形状)が異なるので独立とは言えない.
12 / 13
大レポート
1. 東日本大震災以外の地震データを取得し,演習で行ったのと
同様な解析を行い結果を比較しなさい.データの出典は明記
すること.
2. トヨタ自動車以外の株価データを取得し,演習で行ったのと
同様な解析を行い結果を比較しなさい.データの出典は明記
すること.
1, 2 のいずれかを行えばよい.
締切:6 月 6 日(金)の演習開始時
注:地震データについては,東大地震研のページ
http://eoc.eri.u-tokyo.ac.jp/harvest/
などから入手してください.
13 / 13