Excelのソルバー

今日は二通りの統計モデルを考えて、最尤推定量を求めます。
都市群と農村群で共通のパラメータを用いるモデルと、別のパラメータを用いるモデルです。

まず先週のデータを貼り付けます。ひとまとめにして扱うので、都市群のデータの下に農村群のデータを貼り付けます。

最初は、共通のパラメータを使うモデルを考えましょう。
C1セルにλの値としてとりあえず1を代入し、先週の最後のようにB2セルにlog(f(x_α|λ))の式を書いてください。ここでポアソン分布の確率関数はf(x|λ)=λ^xe^-λ/x!で、xの値はA2に、λの値はC1にあります。この式を下へコピーしますので、C1の方は番号がずれないようにC$1と書いてください。

書いたら、下へコピーします。農村群の部分にもコピーしてください。対数尤度はそれらの値を全部足したものですので、B52セルに =SUM(B2:B27)+SUM(B29:B51) と書いてください。

B52セルの値が一番大きくなるようにC1セルの値を調整します。色々入れて大きくなるのを探すのは面倒ですので、Excelのソルバーを使います。
画面の上の「データ」をクリックして、右端の「ソルバー」をクリックします。もしソルバーがなかったら、先週のページを見てインストールします。

目的セル(E)欄に、最大にしたい対数尤度が書かれているB52、変化させるセル(B)欄にパラメータλが書かれているC1と書きます。欄をクリックしてからそれらのセルをクリックしても構いません。
制約条件(U)欄には、本来ならパラメータλは正でなければなりませんのでその条件を書くべきですが、とりあえず実行してみてうまくいかなかったらその時制約条件を入力することにして、実行(S)をクリックします。

うまく最適値が見つかりました。解を記入する(K)を選んでOKをクリックすると、対数尤度の最大値がB52に、そして対数尤度が最大になるλの値、つまり最尤推定値がC1に入力されます。

次に都市群と農村群で別のパラメータを用いるモデルを考えます。
E1セルに都市群のパラメータとしてとりあえず1を、E28セルに農村群のパラメータとしてとりあえず1を入力して、D2セルにデータxとしては先ほど同様A2を、パラメータλとしてE1セルを使う式を書いてD27までコピーします。

今度はD29セルに、データxとしてはA29を、パラメータλとしてE28セルを使う式を書いてD51までコピーします。

D52セルにそれらの合計、つまり =SUM(D2:D27)+SUM(D29:D51) を書いて、この値を最大にするようにE1セルとE28セルの数字を変化させます。
SUM(D2:D27)を最大にするE1と、SUM(D29:D51)を最大にするE28を探せばよいのですが、これもソルバーを使いましょう。
最大にする目的セルはD52, 変化させるセルはE1,E28です。

対数尤度の最大値が求まりました。この大きさだけをみて、大きいほうが良い、と考えてはいけません。

AICを求めます。AIC=-2×対数尤度の最大値+2×パラメータ数です。これが小さいほうが良いモデルです。

別のパラメータを使うモデルのほうがAICが小さいので良いモデルといえます。実はこの場合は値の差が1以下なので断言し辛いですが、今回は練習ですので気にしないことにします。

練習問題

男の子	女の子
期間(週) 体重(g) 40 2968 38 2795 40 3163 35 2925 36 2625 37 2847 41 3292 40 3473 37 2628 38 3176 40 3421 38 2975 平均 3024.0 標準偏差 272.1	期間(週) 体重(g) 40 3317 36 2729 40 2935 38 2754 42 3210 39 2817 40 3126 37 2539 36 2412 38 2991 39 2875 40 3231 平均 2911.3 標準偏差 268.5	出典

このデータに対して次の三つのモデルを考えます。AICを用いてモデルを選びなさい。但しxは妊娠期間です。
モデル1:体重の分布は男女とも正規分布N(α+βx,σ²)に従う。
モデル2:男の子の体重の分布は正規分布N(α₁+βx,σ²), 女の子の体重は正規分布N(α₂+βx,σ²)に従う。βとσは男女共通で、切片αだけ異なる。つまり平行線を当てはめる。
モデル3:男の子の体重の分布は正規分布N(α₁+β₁x,σ²), 女の子の体重は正規分布N(α₂+β₂x,σ²)に従う。σは男女共通で、切片αと傾きβが異なる。つまり全く別の線を当てはめる。

正規分布N(μ,σ²)の密度関数はf(x)=1/√(2πσ²) exp(-(x-μ)²/2σ²)です。
πはExcelでは=PI()
√xはExcelでは=SQRT(x)です。

注意点

σ^2の値をパラメータにすると、ソルバーが負の数を代入するとエラーになります。範囲に制限をつけても良いのですが、簡単にするために、σの値をパラメータにします。するとソルバーが負の値を代入しても、2乗して正の数になるので大丈夫です。
Excelで -A1^2 と書くと(-A1)*(-A1)の意味になります。-A1*A1の意味にするには -(A1^2) と書いてください。
A1/B1*C1 と書くと、A1をB1で割ってC1を掛けることになります。B1*C1でA1を割るには A1/(B1*C1) と書いてください。
パラメータの値として最初に代入する値があまりにもデータと離れすぎていると、尤度が殆ど0になって対数尤度がマイナス無限大になってしまいます。大雑把に水平な線を当てはめることを考えて、α=3000, β=0, σ=100を最初に代入してください。