最尤推定の手計算:転載
◆最尤推定はイメージがしづらい.手計算で数字を見ながら様子が実感できるシートを作成した.glm 重回帰では、純粋に手計算で実感するのは無理、表計算なら幾分の実感を持てる.
◆対数尤度は、
lnL = Σln conbine(N,yi) + Σ{ yi ln(q) + (N-yi) ln(1-q) } [1].
さらに、qについて係数ごとに微分し、各回(個体)について計算値を総和し、連立方程式を解くことになるが、非常に大変[2].(目的は、最尤推定を実感することにある)
上式の右辺をq微分して、式を簡潔にする.これが、計算式の組み立てを楽にすることに気付く.
線形予測子とその指数関数からなる式が導出できるから、これを計算シートで実現する.
・生起y、暴露xを数値表とし、線形予測子(切片α、xiの係数としてβi)の係数をそれぞれ参照させ、個体ごとの指数関数計算式を設定する.これを全個体について総和する.
・係数を参照する1組のセル群を設定する.いったんは、予めglmにより決めた係数を入れる.これを加減して入力したものの総和を観察すれば、与える係数によって微分対数尤度が非常に敏感に変化する様子を実感できる.
・例によって、重回帰ロジスティックとし、暴露は01の2値による場合(2×2表)とした.
データのソースは、とりあえず観光事例から5因子をもとにした.
glmから得てある(真)値に近づくと、この値がゼロに近づくが、小数点以下での与える係数変化により総和が微妙に変化する.ニュートンラフソン法を手探りする感覚だ.
・MHによりいったんORを概略計算し、係数を得ておき、このシートで細密に推定値を手探りできるという使い道もある(としてもglmが楽!).
---------------------------------------------------
memo
尤度;L
L=Πf()=(1-pi)pi・・=(1-1/(1+exp(-z))・1/(1+exp(-z’))・・・
ここでz:線形予測子 z(')=β。+β1x1+β2x2・・
対数尤度
Lを対数とすると、線形予測子の係数のうち、生起のものの数(a1,a2・・)だけ表にでてくる
logL = - log(1+exp(-z)) -log( )- ・・・+β。+β1a1+β2a2・・
微分すると、log(1+exp(-z))は、係数を背負うが、2値の係数なので暴露ありのpi和になる
δlogL/δβx = -pi - pi -pi ・・ + δβ。/δβx+δβ1a1/δβx+δβ2a2/δβx・・ = 0
= -exposureΣpi + (微分する項のy=1となるxの個数)=0
---------------------------------------------------
参考文献
[1] 久保氏 「統計数理研究所公開講座」資料2011
[2] 土居氏 一般化線形モデル入門の入門