数理統計学(第六回） 最尤推定とは？
Download
Report
Transcript 数理統計学(第六回） 最尤推定とは？
数理統計学(第六回）
最尤推定とは？
浜田知久馬
数理統計学第６回
1
数理統計学第６回
2
数理統計学第６回
3
数理統計学第６回
4
ダーウィンの例 母数推定の前提
自家受精群と他家受精群 に別々の正規分布を
あてはめ
ｎ個（ｎ＝１５）の確率変数Ｙｉが互いに独
立に同一の正規分布にしたがう
Ｙ１ ,Ｙ２ ,Ｙ３ ,・・・，Ｙｎ ～Ｎ(μ,σ２)
i.i.d.（independent identically distributed）
数理統計学第６回
5
正規分布の確率密度関数
2
1
y
f ( y)
exp
2
2
は既知
2
2
σ２
ｎ個Ｙ１ ，・・・，Ｙｎ のｎ個のデータの得られる確率ｆ
f＝ｆ(y１) ・ｆ(y2) ・・・ｆ(yｎ) ＝Πｆ(yi)
n
f
i 1
1
2 2
yi 2
exp
2
2
2
n
y
1
i
exp
2
2
2
2 数理統計学第６回
i 1
n
6
対数尤度(log likelihood)
2
y
1
i
l log f log
exp
2
2
2
i 1
2
2
n
yi
n
2
log(2 )
2
2
2
i 1
n
n
y
yi
2
2
i 1
y
n
i 1
2
y
2
2
i
n
2
i 1
y
2
2
i
y
2
2
n y
2
2数理統計学第６回
2
n
i 1
2
2
7
対数尤度(log likelihood)
n y
l
2
2
2
Lはμについての2次関数
尤度ｆの最大化⇒ 対数尤度Lの最大化
⇒ｄL/ｄμ＝0となるμを探す.
数理統計学第６回
8
スコア統計量
Ｒａｏの有効スコア統計量
（efficiency score function）
対数尤度をパラメータで微分した統計量
d log L
n y
U ( )
d
2
U ( ) y
自家受精群： 17.708 ，他家受精群： 20.192
数理統計学第６回
9
尤度の計算プログラム
data mle;set mle;
do m=15 to 25 by 0.1;s=3;
f=1/(2*3.141728*s**2)**.5*exp(-(y1-m)**2/s**2/2);
l=log(f);
output;end;
proc sort;by m;
proc summary data=mle ;var l;output out=out
sum=;by m;
data out;set out;f=exp(l);
proc gplot;
plot (f l)*m/href=17.708;symbol1 i=spline;run;
数理統計学第６回
10
尤度関数(自家受精群）
数理統計学第６回
11
対数尤度関数(自家受精群）
n y
l C
2
2
数理統計学第６回
2
12
スコア統計量の性質を導くために
利用すること
＾
E[U]=0，V
E[-U']=I =1/V[θ]
1) 確率密度関数の和は１ ∫ｆ(ｙ,θ) ｄｙ=1
2) d log f ( y, ) 1 df ( y, ) df ( y, ) f ( y, ) d log f ( y, )
[U]=E[U2]=
d
３）
４）
５）
６）
f ( y, )
d
d
d
微分と積分の交換可能性
（ｆ・ｇ）'＝ ｆ'・ｇ+ｆ・ｇ'
（ｆ/ｇ）'＝ (ｆ'・ｇ-ｆ・ｇ')/ｇ2
V[X]＝E[X2]－E[X]2
数理統計学第６回
13
スコア統計量の性質
ｆ（ｙ；θ）：確率変数Yの確率密度関数
L（θ；ｙ）＝logｆ（ｙ；θ）：対数尤度関数
dL( ; y)
d log f ( y; )
U
d
d
1
df ( y; )
: スコア関数
f ( y; )
d
dL( ; y)
E[U ]
f ( y; )dy
d
df ( y; )
dy
d
14
数理統計学第６回
スコア統計量の性質
2
E[U]=0，V [U]=E[U ]= E[-U’]=I
f ( y; )dy 1 : 確率の和は１
両辺をで微分
d f ( y; )dy
df ( y; )
dy 0
d
d
E[U ] 0 傾きの期待値は０
数理統計学第６回
15
Uボートで尤度の
山を一周したら
数理統計学第６回
16
スコア統計量の性質
E[U ] 0なので
dL( ; y)
d
f ( y; )dy
dE[U ]
d
0
d
d
2
d L( ; y)
dL( ; y) df ( y; )
f ( y; )dy
dy
2
d
d
d
数理統計学第６回
17
スコア統計量の性質
d L( ; y)
d 2 f ( y; )dy U ' f ( y; )dy E[U ' ]
2
dL( ; y) df ( y; )
dL( ; y) dL( ; y)
d d dy d d f ( y; )dy
2
2
U f ( y; )dy E[U ]
Ｅ[U’] +Ｅ[U 2 ]＝0
V [U]＝E[U2]‐ E[U] 2数理統計学第６回
＝ E[U2]＝ E[-U’]：情報量18
どちらの対数尤度の山の方が登りやすい
だろうか？
情報量が小さい場合
情報量が
大きい場合
数理統計学第６回
19
2項分布の場合
y
n y
L n C y (1 )
log L log n C y y log (n y) log(1 )
d log L d ( y log (n y) log(1 ))
U
d
d
y n y y n
1 (1 )
y n E[ y n ]
E[U ] E
0
(1 ) (1 )
数理統計学第６回
20
2項分布の場合
V [U ] I E[U ]
2
2
y n
E[( y n ) ]
2
2
E[U ] E
2
(1 ) (1 )
V [ y]
n (1 )
2
2
2
2
(1 ) (1 )
n
1
(1 ) V [ ]
2
数理統計学第６回
21
2項分布の場合
y n
U
(1 )
( y n )' (1 ) ( y n ) (1 )'
U'
2 (1 ) 2
n (1 ) ( y n )(1 2 )
2 (1 ) 2
n n 2 y n 2y 2n 2
2 (1 ) 2
n 2 y 2y
2 (1 ) 2 数理統計学第６回
22
2項分布の場合
n 2 y 2y
E[U ' ] E
2
2
(1 )
2
2
n n 2n
2
2
(1 )
n n
n (1 )
2
2
2
2
(1 )
(1 )
2
V [U ] E[U ] I
2
数理統計学第６回
23
＾
＾
V[θ]≒1／I[θ]
f (a)(x a) 2
f ( x) ≒ f (a) f (a)(x a)
2
U ( )を の周りでテーラー展開 して一次近似
U ( ) ≒U ( ) U ' ( )( ) I ( )
U ( )
E[U ( )]
( ) ≒
, E[ ] ≒
0
I
I
V [U ( )] I 1
V [ ] ≒
2
2
I
I
I
数理統計学第６回
中心極限定理と
大数の法則が
適用できる．
分散は
1/情報量24
平均値の法則
d log f i
d log f i
U[ ]
d
d
大数の法則
平均値はｎを大きくすると，真の値に収束する.
平均値→E(X)＝μ （ｎ→∞）
limP（|平均値－μ|≧ε）＝0
ｎ→∞
中心極限定理
ｎを大きくすると，平均値の分布は正規分布になる.
数理統計学第６回
25
まとめ スコア統計量とMLEの性質
d log f
, E[U ] 0
U [ ]
d
2
V [U ] E[U ] E[U ' ] I
U ( ) ≒ I ( )
E[ ] ,V [ ] 1 / I
Uも も漸近的に正規分布
数理統計学第６回
26
MLEの性質
１）ｎが大きくなれば，MLEは真値に近づ
く(一致性） ← 大数の法則
２）最尤推定量の分布は，漸近的に正規分
布にしたがう(漸近正規性）
← 中心極限定理
３）最尤推定量の分散は，漸近的にFisher
の情報量の逆数となり,Cramer-Raoの下
限を達成する（有効性）
数理統計学第６回
27
ダイオキシン
ベトナム戦争時に使用された枯葉剤（２，４，
５－Ｔ）に、不純物として含まれていた毒性物質
として有名。ベトナムにおける流産や先天異常児
の多発、ベトナム帰還米兵のガンの多発などは、
この影響であると指摘されている。
ダイオキシンは廃棄物焼却や金属精錬の際に発
生したり、農薬等各種の化学物質を製造する際に
副生するなどの例が知られているが、中でも廃棄
物焼却が最大の原因であると言われている。とこ
ろが、最近になって、これまで知られていた発ガ
ン性を示すレベルの量よりも更に微量のダイオキ
シンが、子宮内膜症とそれに伴う不妊症の原因で
はないかと疑われるようになり、社会的な関心が
28
一層高まっている。数理統計学第６回
ベトナム戦争で散布された枯葉剤に
より枯れたマングローブの林
数理統計学第６回
29
ダイオキシン類の特徴
４塩化物以上に毒性が強い
2,3,7,8－ＴＣＤＤの毒性が最も強い
水に難溶性。脂肪に溶けやすい
極めて安定 ７００℃でも分解しない
副産物として生成される
除草剤、殺菌剤の製造過程、焼却、
漂白、自動車の排ガス等
数理統計学第６回
30
2,3,7,8-Tetrachlorodibenzo-p-dioxin
(2,3,7,8-TCDD) の構造
塩素のつく数と位置による同族体と異性体が
多数あり、このうち毒性発現から７種、またポリ
塩化ジベンゾフラン10種、PCB12種に対して毒
性等価係数が付与されている）
数理統計学第６回
31
測定単位
• mg （ﾐﾘｸﾞﾗﾑ ＝千分の１ｸﾞﾗﾑ）
μg （ﾏｲｸﾛﾘｸﾞﾗﾑ ＝百万分の１ｸﾞﾗﾑ）
ng （ﾅﾉｸﾞﾗﾑ＝十億分の１ｸﾞﾗﾑ)
pg （ﾋﾟｺｸﾞﾗﾑ＝一兆分の１ｸﾞﾗﾑ）
• 1g中に1 μg → ｐｐｍ
1g中に1 ng → ｐｐｂ
1g中に1 pg → ｐｐｔ
(100m×100m×100mの水槽に１ｇ）
数理統計学第６回
32
ダイオキシンの毒性
毒性
症
状
実験動
物
ＮＯＡＥＬ＊
発ガン性
肝細胞への影響
肝臓ガン発生
ラット
ラット
1 ng/1kg 体重/1日
10ng/1kg 体重/1日
免疫毒
性
感染防御機能低
下
マウス
５ng /1kg 体重/1日
催奇形
性
口蓋裂、腎盂拡
張
マウス
100ng /1kg 体重/1
日
その他の
慢性毒
性
体重増加抑制等
ラット
1 ng/1kg 体重/1日
数理統計学第６回
33
性比の経年変化 男/女
Fig1
Yearly Trend of R atios in Japan
108.0
107.5
Sex Ratio
107.0
106.5
106.0
105.5
105.0
Sex Ratio
104.5
Moving Average
104.0
1945
1950
1955
1960
1965
1970
1975
Year
数理統計学第６回
1980
1985
1990
1995
34
問題の背景
1) ダイオキシンの人間に対する毒性用量は
判明している.
2）ダイオキシン濃度がある濃度を越えた個
体の割合を推定したい.
3)1998年については個別データがあるが，
後の年については，平均値しかわからない
4）ダイオキシン濃度の分布は対数正規分布
にしたがうことが経験的に判明している．
数理統計学第６回
35
セベソ事件
セベソ事件とは，北イタリア，ミラノの北約 20 Km
のところにある農薬工場が，2,4,5-トリクロロフェ
ノール製造中に爆発事故を起こし，セベソを中心
とする周辺地域をダイオキシンを含む大量の化
学物質が汚染した．
・1976 年 7 月 10 日爆発事故発生
・最終的に，２０万人以上の住民が被害を受ける.
･性比を1%減らす濃度(129 pg/g）
5%(146)
10%(160）
数理統計学第６回
36
モーメント法による
対数正規分布の母数推定
N=467 （単位：pg/g fat)
平均値：25.5
SD：14.5
変動係数：14.5/25.5＝56.9％
対数正規分布のパラメータ
μ,σ
数理統計学第６回
37
対数正規分布
対数変換すると正規分布にしたがう.
y
f ( y)
exp
2
2
2
2
2
1
ｘ＝log(ｙ) ⇒ ｄｘ＝1/ｙｄｙ
ｆ(ｙ) ＝ｄＦ(ｙ) /ｄｙ＝ｄＦ(ｙ) /ｄｘ
×ｄｘ/ｄｙ
2
1
log
y
f ( y)
exp
2
2
2
y 2
μ,σを推定したい
数理統計学第６回
38
対数正規分布
E［ｙ］＝exp(σ２/2＋μ）
V［ｙ］＝E［ｙ］２ ( exp(σ２)－1)
SD ＝ E［ｙ］( exp(σ２)－1)0.5
CV＝SD/E［ｙ］＝ ( exp(σ２)－1)0.5
変動係数はμに関係なく一定
数理統計学第６回
39
モーメント法による推定
E［ｙ］＝exp(σ２/2＋μ)＝標本平均
V［ｙ］＝E［ｙ］２ ( exp(σ２)－1)＝標本分散
になるように， σとμを推定
標本平均：25.5 標本SD：14.5
exp(σ２)－1＝14.5２/25.5２⇒σ=0.5293
exp(σ２/2＋μ)＝25.5 ⇒ μ＝3.0986
数理統計学第６回
40
推定された対数正規分布の
確率密度関数
数理統計学第６回
41
対数正規乱数の発生
data data;
cv=0.5686;mean=25.5;
s=(log(1+cv**2))**.5;
m=log(mean)-s**2/2;
do i=1 to 10000;
x=exp(m+s*rannor(4989));y=log(x);output;
end;
proc means;var x y;
proc gchart;vbar x/space=0;
数理統計学第６回
42
正規乱数のヒストグラム
数理統計学第６回
43
対数正規分布の％点
log(ｙ)が正規分布するのを利用して
α％点：ｅxp(μ+Ｚασ)
Ｚα：正規分布のα％点
Ｚ50：0
Ｚ25：-0.6745 Ｚ75：0.6745
数理統計学第６回
44
％点の比較
実測値
対数正規分布
25％点 15.2
15.5
50％点 21.7
22.2
75%点
31.3
31.7
%点からは対数正規分布がよくあてはまってい
ると考えられる.
変動係数が一定になるように各年度のＳＤを推
定
数理統計学第６回
45
カットオフ値を越える確率
（単位：pg/g）
year
1973
1974
1975
1976
1977
1978
1979
1980
1981
1982
1983
平均
57.1
66.2
58.2
48.9
48.7
51.9
46.6
40.0
38.1
38.6
40.6
σ
0.5293
0.5293
0.5293
0.5293
0.5293
0.5293
0.5293
0.5293
0.5293
0.5293
0.5293
μ
3.90472
4.05260
3.92381
3.74970
3.74560
3.80924
3.70152
3.54880
3.50014
3.51317
3.56369
推定SD BMD1%
（129）
32.4687 0.035582
37.6433 0.063623
33.0942 0.038498
27.8060 0.017982
27.6923 0.017643
29.5119 0.023581
26.4981 0.014322
22.7452 0.006627
21.6648 0.005102
21.9491 0.005477
23.0864 0.007168
数理統計学第６回
BMD5%
（146)
0.020759
0.039294
0.022628
0.009871
0.009669
0.013253
0.007713
0.003356
0.002534
0.002734
0.003652
BMD10%
（160）
0.013507
0.026684
0.014805
0.006136
0.006003
0.008385
0.004726
0.001965
0.001462
0.001583
0.002148
46
カットオフ値を越える確率
（単位：pg/g）
year
1984
1985
1986
1988
1989
1990
1991
1992
1993
1994
1995
1996
平均
37.3
32.4
30.5
34.4
33.0
31.9
29.2
26.2
26.5
27.1
23.5
24.1
σ
0.5293
0.5293
0.5293
0.5293
0.5293
0.5293
0.5293
0.5293
0.5293
0.5293
0.5293
0.5293
μ
推定SD
3.47891
3.33808
3.27765
3.39798
3.35643
3.32253
3.23409
3.12568
3.13707
3.15945
3.01692
3.04213
21.2099
18.4236
17.3432
19.5609
18.7648
18.1393
16.6040
14.8981
15.0687
15.4099
13.3628
13.7040
BMD1%
（129）
0.004541
0.002020
0.001399
0.002874
0.002253
0.001840
0.001065
0.000526
0.000567
0.000658
0.000249
0.000297
数理統計学第６回
BMD5%
（146)
0.002236
0.000939
0.000634
0.001369
0.001055
0.000850
0.000474
0.000224
0.000243
0.000284
0.000101
0.000122
BMD10%
（160）
0.001282
0.000516
0.000342
0.000766
0.000583
0.000464
0.000252
0.000115
0.000125
0.000148
0.000050
0.000061
47
問題
ある大学で，100人の学生の血液型について，調
べた結果，A型41人，B型22人，O型27人，AB
型10人となった．血液型の人数の分布は次の
ように確率Lが示される多項分布に
したがうと考えられる．
( XA XB XO XAB)! XA XB XO XAB
L
A B O AB
XA! XB! XO! XAB!
XA,･･･,XAB：各血液型の人数である．
XA+XB+XO+XAB＝N
数理統計学第６回
48
問題
Lを母数, πA , πB , πO , πABの関数を考えてL
が最大になるようにパラメータを推定したい．
ただしπA+πB +πO + πAB=1
という制約が存在する．
1)対数尤度を示すこと．
2)対数尤度をπAで微分すること．
3) πA , πB , πO , πABの最尤推定値を計算する
こと
数理統計学第６回
49
Lagrange の未定乗数法
• p 次元ベクトル x について，等式制約：g(x) = 0
の下で，ある目的関数：f(x) の最大（小）値を求
める問題の一つの解法に，ラグランジュの未定
乗数法がある．
Q f ( x) g ( x) とおき，連立方程式
dQ
dQ
0, i 1,2,...p;
0
dxi
d
を解くと，この連立方程式の解の中に求める解が
50
数理統計学第６回
ある．