数理統計学(第六回）最尤推定とは？

数理統計学(第六回）
最尤推定とは？
浜田知久馬
数理統計学第６回
1
数理統計学第６回
2
数理統計学第６回
3
数理統計学第６回
4
ダーウィンの例母数推定の前提
自家受精群と他家受精群に別々の正規分布を
あてはめ
ｎ個（ｎ＝１５）の確率変数Ｙｉが互いに独
立に同一の正規分布にしたがう
Ｙ１ ,Ｙ２ ,Ｙ３ ,・・・，Ｙｎ～Ｎ(μ,σ２)
i.i.d.（independent identically distributed）
数理統計学第６回
5
正規分布の確率密度関数
f ( y) 
2


1
y  
exp 
2
2
2
2





σ２は既知
ｎ個Ｙ１，・・・，Ｙｎのｎ個のデータの得られる確率ｆ
f＝ｆ(y１) ・ｆ(y2) ・・・ｆ(yｎ) ＝Πｆ(yi)
n
f 
i 1
1
2 2
  y i   2
exp 
2
2

2
n




yi   
1

exp  

2
2
2
 2  数理統計学第６回
 i 1
n








6
対数尤度(log likelihood)
2



yi   
1

l  log f   log
exp 
2
2
2


i 1

 2
n
 yi   
n
2
 log(2 )  
2
2
2
i 1
n
n

 yi   
2
i 1
n

i 1
y
2

2
y
2
2
i
n

2
i 1
y

y
2
2
i


2





2
n

i 1
y   
2
2
2
2
n y

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
  U f ( y; )dy  E[U ]
2
2
Ｅ[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  logn 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
2
 (1   )


 n  n
 n
2
 y  n  2y  2n
 2 (1   ) 2
 y  2y
2
2
 (1   ) 数理統計学第６回
2
2
22
2項分布の場合
  n 2  y  2y 
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
平均値の法則
U [ ] 
d  log f i
d
d log f i

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
Y early T rend 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
Y ear
数理統計学第６回
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
対数正規分布
対数変換すると正規分布にしたがう.
f ( y) 
2


1
y  
exp 
2
2
2
2

ｘ＝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
数理統計学第６回
ある．

Download Report