ソルバーを使って最尤推定量を求める

今日のテーマ

最尤推定量は(対数)尤度を最大にする値としてデータから計算できます。この計算式θ=θ(X₁,X₂,...,X_n)を具体的に書き表すことが出来る問題もあります(以前のレポートで出しました)が、計算式を求めることが出来ない問題もあります。
そこで今日はExcelを使って最尤推定量を求めて、AICを計算し、モデルを選ぶ方法を勉強しましょう。

準備

Excelを立ち上げたら左上の「ファイル」を押して、「オプション」をクリックします。そして左側の「アドイン」をクリックして、下の「管理:Excelアドイン」の横の「設定」をクリックします。

するとアドインの画面に切り替わるので「ソルバーアドイン」にチェックをつけてOKをクリックします。

サンプルデータ

男児 女児 期間(週) 体重(g) 期間(週) 体重(g) 40 2968 40 3317 38 2795 36 2729 40 3163 40 2935 35 2925 38 2754 36 2625 42 3210 37 2847 39 2817 41 3292 40 3126 40 3473 37 2539 37 2628 36 2412 38 3176 38 2991 40 3421 39 2875 38 2975 40 3231

出典:一般化線形モデル入門, Annette J.Dobson 著, 共立出版
これをコピーしてエクセルに貼り付けてください。散布図を書くとこのようになります。

考えるモデル

パラメータが少ない順番に

男女共通モデル

男女ともY=α+βX+ε, ε～N(0,σ²)
つまり妊娠期間Xに対して体重Yの確率分布は
男女ともN(α+βX,σ²)
推定するパラメータはα, β, σの3つ。

傾きだけ共通モデル

男児:Y=α₁+βX+ε, ε～N(0,σ²)
女児:Y=α₂+βX+ε, ε～N(0,σ²)
つまり妊娠期間Xに対して体重Yの確率分布は
男児:N(α₁+βX,σ²)
女児:N(α₂+βX,σ²)
推定するパラメータはα₁, α₂, β, σの4つ。

男女別モデル

男児:Y=α₁+β₁X+ε, ε～N(0,σ²)
女児:Y=α₂+β₂X+ε, ε～N(0,σ²)
つまり妊娠期間Xに対して体重Yの確率分布は
男児:N(α₁+β₁X,σ²)
女児:N(α₂+β₂X,σ²)
推定するパラメータはα₁, α₂, β₁, β₂, σの5つ。

まず最初に、パラメータが一番少ない男女共通モデルについて、α, β, σをExcelに求めてもらいます。
そのためには、まず最初に適当な初期値を与えます。すると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=1ⁿf(X_i))=Σ_i=1ⁿlog(f(X_i))
正規分布N(μ,σ²)の密度関数は =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を変化させる、と指定します。
また、αは正の数とは限らないので「制約のない変数を非負数にする(K)」のチェックを外します。さらにソルバーの設定を変えるのでオプションを押します。

「GRG非線形」のタブを選んで、微分係数を「中央」にします。

OKを押してソルバーの画面に戻って「解決(S)」を押すと、対数尤度が大きくなるようにα, β, σの値が調整されて次のメッセージが表示されます。

OKをクリックすると、対数尤度の値が最大になっています。

AICの計算

このようにして求めた対数尤度の最大値からパラメーターの個数を引いた値が一番大きなモデルが、一番良いモデルなのですが、この方法が考えられた経緯から、-2倍して
-2×(対数尤度の最大値)+2×(パラメータの個数)
をAIC (Akaike information criterion)といい、マイナスを掛けているのでAICが一番小さなモデルが一番良いモデルです。
このモデルのAICは-2×(-159.265)+2×3=324.53です。

レポート課題

傾きだけ共通モデル、及び男女別モデルについて、最尤推定量及びAICの値を求め、どのモデルが一番良いか選びなさい。