今日のテーマ
最尤推定量は(対数)尤度を最大にする値としてデータから計算できます。この計算式θ=θ(X1,X2,...,Xn)を具体的に書き表すことが出来る問題もあります(以前の小テストで出しました)が、計算を間違えることもありますし、どんなに頑張っても計算式を求めることが出来ない問題もあります。そこで今日はExcelを使って最尤推定量を求めて、AICを計算し、モデルを選ぶ方法を勉強しましょう。
準備
Excelを立ち上げたら左上の丸いボタンを押して、「Excelのオプション」をクリックします。そして左側の「アドイン」をクリックして、下の「管理:Excelアドイン」の横の「設定」をクリックします。
するとアドインの画面に切り替わるので「ソルバーアドイン」にチェックをつけてOKをクリックします。 「インストールしますか」と聞かれたら「はい」と答えます。
サンプルデータ
| |||||||||||||||||||||||||||||||||||||||||||||||||||||||||
これをコピーしてエクセルに貼り付けてください。散布図を書くとこのようになります。
考えるモデル
パラメータが少ない順番に- 男女共通モデル
- 男女ともY=α+βX+ε, ε～N(0,σ2)
つまり妊娠期間Xに対して体重Yの確率分布は
男女ともN(α+βX,σ2)
推定するパラメータはα, β, σの3つ。 - 傾きだけ共通モデル
- 男児:Y=α1+βX+ε, ε～N(0,σ2)
女児:Y=α2+βX+ε, ε～N(0,σ2)
つまり妊娠期間Xに対して体重Yの確率分布は
男児:N(α1+βX,σ2)
女児:N(α2+βX,σ2)
推定するパラメータはα1, α2, β, σの4つ。 - 男女別モデル
- 男児:Y=α1+β1X+ε, ε～N(0,σ2)
女児:Y=α2+β2X+ε, ε～N(0,σ2)
つまり妊娠期間Xに対して体重Yの確率分布は
男児:N(α1+β1X,σ2)
女児:N(α2+β2X,σ2)
推定するパラメータはα1, α2, β1, β2, σの5つ。
そのためには、まず最初に適当な初期値を与えます。するとExcelのソルバーはそこから出発して真の値に近づいていきます。初期値をどのように与えればよいでしょうか。最初から真の値が分かるならExcelのソルバーに頼る必要はないのですから、真の値に近くてしかも簡単に求められる値を初期値として使います。
直線y=α+βxのαとβの両方を簡単に知ることは出来ませんが、水平な線、つまりβ=0ならば、その高さαは平均値にしておけば大体データに近そうです。この時σは平均値からのずれ、つまり標準偏差にしておきましょう。
そこで早速体重の平均値と標準偏差を求めましょう。Excelに次のように貼り付けます。
A15に平均、A16に標準偏差と書いてから、B15に平均を求める式、B16に標準偏差を求める式を書きます。
B15にはB3からB14までとD3からD14までの平均ということで =AVERAGE(B3:B14,D3:D14) と書きます。先頭の等号を忘れないでください。
B16にはB3からB14までとD3からD14までの標準偏差ということで =STDEV(B3:B14,D3:D14) と書きます。
するとこのようになります。
セルの幅に合わせて四捨五入されますので、人によって表示される値が少し違うこともあります。
この値を参考にしてα, β, σの初期値を決めます。
まずA17, A18, A19セルにα, β, σと書いて
B17セルにはαの初期値として平均値を書きます。適当なところで四捨五入してかまいません。
B18セルにはβの初期値を書きます。傾きが良く分からなくて水平線を引くので0にします。
B19セルにはσの初期値として標準偏差を書きます。これも適当なところで四捨五入してかまいません。
対数尤度の計算
確率関数や密度関数にデータを代入した値が尤度、それにlogをつけたのが対数尤度です。独立な複数の観測値に関して尤度は掛け算になっているので、logをつけて掛け算を足し算にします。log(Πi=1nf(Xi))=Σi=1nlog(f(Xi))
正規分布N(μ,σ2)の密度関数は =NORMDIST(x, μ, σ, 0) です。最後の0を1に変えると分布関数になります。
今回の問題では、xのところが体重yで、平均値μがα+βxです。それにlogを付けますが、Excelでは自然対数は =LN(x) ですのでE3セルに =LN(NORMDIST(B3,B$17+B$18*A3,B$19,0)) と書きます。対応するセルに色がつくので確認してください。
α, β, σに相当するB17, B18, B19の数字の前には$を付けています。これはこの式をコピーして下に貼り付けても番号が変わらないようにするためです。
E4セルからE14セルにも同様に数式を書きます。一つ一つ書いていると手間がかかりますのでE3セルをコピーして貼り付けます。数式の中のセルの番号がどのように変わるか確認してください。
次にF3セルに女児の一人目の対数尤度の式 =LN(NORMDIST(D3,B$17+B$18*C3,B$19,0)) を書きます。
これもF4からF14へコピーして、E3からF14まで全部足したものが対数尤度です。
E15セルに対数尤度と書いて、F15セルに合計を求める式 =SUM(E3:F14) を書きます。
ソルバー
ここが今回のポイントです。これからα, β, σに色々な値を代入して対数尤度が最大になればよいのですが、試行錯誤するのは大変ですのでExcelのソルバーを使います。エクセルの上の列の「データ」をクリックして、右側の「ソルバー」をクリックします。「ソルバー」が見当たらなかったら、このページの上の準備のところを見てください。
するとこのような画面が表示されるので、F15を最大にするためにB17からB19を変化させる、と指定します。
「実行」をクリックすると対数尤度が大きくなるようにα, β, σの値が調整されて次のメッセージが表示されます。
OKをクリックすると、対数尤度の値が最大になっています。
AICの計算
このようにして求めた対数尤度の最大値からパラメーターの個数を引いた値が一番大きなモデルが、一番良いモデルなのですが、この方法が考えられた経緯から、-2×(対数尤度の最大値)+2×(パラメータの個数)
をAIC (Akaike information criterion)といい、マイナスを掛けているのでAICが一番小さなモデルが一番良いモデルです。
このモデルのAICは-2×(-159.265)+2×3=324.53です。
