ポアソン分布

5月22日の講義に引き続いて、Excelのソルバーを用いて最尤推定量を求めます。5月22日に欠席してソルバーの準備をしていない人は5月22日の講義のページを見て準備してください。
注:5月22日の講義の時と別のパソコンを使っている人は、改めてソルバーの準備が必要かも知れません。

教科書のデータ

これまでの講義で何度か使った、都市部と農村部の慢性病の数のデータです。 一般化線形モデル入門22ページ表2.1のデータです。

都市群 0 1 1 0 2 3 0 1 1 1 1 2 0 1 3 0 1 2 1 3 3 4 1 3 2 0

農村群 2 0 3 0 0 1 1 1 1 0 0 2 2 0 1 2 0 0 1 1 1 0 2

出典

これらのデータが、どんな母集団分布に従って観測されたのか、特に都市群と農村群で同じ分布と考えるのが妥当なのか、違うと考えるのが妥当なのか、考えます。
母集団分布を寸分の狂いも無く推定するには無限個のデータが必要ですのでそれは諦めて、ポアソン分布を用いて推定してみます。他にももっと適当な分布があるかも知れません。
ポアソン分布の確率関数は f(x|λ)=λ^xe^-λ/x! です。そしてそのグラフをいくつかのλに関して書いてみると

になります。今日はまず都市群のデータ

に一番良く当てはまるλの値を最尤法を用いて求めてみましょう。

パラメータのベクトルをθとすると、対数尤度関数は∑_α=1ⁿlog(f(x_α|λ))です。ポアソン分布の場合はパラメータはλ一つですから、θ=λです。ポアソン分布の場合は対数尤度関数をλで微分して=0とおいた式を解くことでλをx₁x₂,…,x_nを用いて表すことも出来ますが、今回は解けない場合に計算機を使ってどうやって求めるか、という問題を考えていますので、Excelで尤度関数の値を求めて、それを最大にするλを探しましょう。

まず都心部のデータをExcelに貼り付けます。説明を合わせるためにA1セルから貼り付けてください。データはA2からA27に入ります。次にλの値を入力します。後で変化させますがまずはλ=1にしましょう。B1セルにλ、C1セルに1と書いてください。
そして対数尤度関数はlog(f(x_α|λ))を足したものですので、B2セルに、この式のxのところにA2セルの値が、λのところにC1セルの値が来るように式を書きます。但し、この式を後でB3以降にコピーしますので、C1はC$1と書いて、下へコピーしても番号が変わらないようにしてください。
ポアソン分布の確率関数の式はEXCELでは =POISSON.DIST(x, λ, 0) です。EXCEL2007以前では =POISSON(x, λ, 0) と書きます。また、底がeである自然対数はLN( )です。

練習

D1セルに都市群のデータの対数尤度関数を計算し、C1セルの値を変化させて、対数尤度が最大になるλの値をソルバーを用いて探してください。
また、対数尤度関数をλで微分して=0とおいた式を解いて、λの最尤推定量をx₁x₂,…,x_nを用いた式で書き、その式に都市群のデータを代入した値と比較してください。

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

まず農村部のデータを貼り付けます。ひとまとめにして扱うので、都市群のデータの下のA28から貼り付けます。23人分のデータがA29からA51に入ります。
最初は、共通のパラメータを使うモデルを考えましょう。
農村部にもB列にlog(f(x_α|λ))の式を書いてください。共通のパラメータを使うのでλはC1セルの値を使います。
対数尤度はそれらの値を全部足したものですので、B52セルに =SUM(B2:B27)+SUM(B29:B51) と書いてください。

目的セル(E)欄に、最大にしたい対数尤度が書かれているB52、変化させるセル(B)欄にパラメータλが書かれているC1と書きます。欄をクリックしてからそれらのセルをクリックしても構いません。
今回は、パラメータλは正でなければなりませんので 制約のない変数を非負数にする(K) にチェックをつけたままにして解決(S)をクリックします。

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

次に都市群と農村群で別のパラメータを用いるモデルを考えます。
E1セルに都市群のパラメータとしてとりあえず先ほど求めた共通のλの値を、E28セルに農村群のパラメータとしてとりあえず同じ値を入力して、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以下なので断言し辛いですが、今回は練習ですので気にしないことにします。