最尤推定法

先週の練習問題の注意点

• σ^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, σは標準偏差に近い値として270を最初に代入してください。

モデル1の最尤推定量を求めるとこのようになります。

レポート課題1

次の3つのモデルの中から、AICを用いて最も適切なモデルを選びなさい。

モデル1

体重の分布は男女とも正規分布N(α+βx,σ²)である。

モデル2

男の子の体重の分布は正規分布N(α₁+βx,σ²), 女の子の体重は正規分布N(α₂+βx,σ²)である。βとσは男女共通で、切片αだけ異なる。

モデル3

男の子の体重の分布は正規分布N(α₁+β₁x,σ²), 女の子の体重は正規分布N(α₂+β₂x,σ²)である。σは男女共通で、切片αと傾きβが異なる。

解答例(わざと見辛くしています。自分で計算してください。)

教科書のデータ

表6.3 20名の男性インスリン依存性糖尿病患者の 炭水化物摂取量, 年齢, 身長から計算された理想体重に対する実際の体重の比率%, 蛋白摂取量	表7.2 かぶと虫の死亡データ ガス濃度はlog10CS2mgl-1
炭水化物 年齢 体重 蛋白 y x1 x2 x3 33 33 100 14 40 47 92 15 37 49 135 18 27 35 144 12 30 46 140 15 43 52 101 15 34 62 95 14 48 23 101 17 30 32 98 15 38 42 105 14 50 31 108 17 51 61 85 19 30 63 130 19 36 40 127 20 41 50 109 15 42 64 107 16 46 56 117 18 24 61 100 13 35 48 118 18 37 28 102 14	ガス濃度 かぶと虫の数 死亡数 x n y 1.6907 59 6 1.7242 60 13 1.7552 62 18 1.7842 56 28 1.8113 63 52 1.8369 59 53 1.861 62 61 1.8839 60 60 出典

レポート課題2

表6.3の炭水化物Yの分布として次の8つのモデルを考えて、AICを用いて選択しなさい。

1. 正規分布N(β₀+β₁x₁+β₂x₂+β₃x₃,σ²)
2. 正規分布N(β₀+β₁x₁+β₂x₂,σ²)
3. 正規分布N(β₀+β₁x₁+β₃x₃,σ²)
4. 正規分布N(β₀+β₂x₂+β₃x₃,σ²)
5. 正規分布N(β₀+β₁x₁,σ²)
6. 正規分布N(β₀+β₂x₂,σ²)
7. 正規分布N(β₀+β₃x₃,σ²)
8. 正規分布N(β₀,σ²)

ロジスティック回帰

表7.2から、ガス濃度毎の死亡割合を計算して図にしました。

このデータに基づいて、死亡確率をガス濃度の関数として表したいと思います。但し確率は0から1の間ですので、直線や多項式を用いると、ガス濃度によっては0から1の範囲をはみ出てしまうので不適です。そこで0から1の範囲をはみ出さないモデルとして次の3つがしばしば使われます。
表の上からi番目の行は、ガス濃度がx_iであるときに、かぶと虫n_i匹の中のx_i匹が死んだということを表しますので、各々のかぶと虫に関して死んだら1、生き残ったら0となるデータをn_i個観測したら、x_i個が1だった、ということになります。従って二項分布B(1,p_i)に従う確率変数をn_i個観測したと考えて、このp_iをx_iのどのような関数で表すか、という問題になります。二項分布の確率関数p(x)=_nC_xp^x(1-p)^n-xで、Excelでは_nC_x=COMBIN(n,x)ですが、₁C₀も₁C₁も1なので省いて構いません。

モデル1:ロジスティックモデル

死亡数の分布は二項分布B(n_i, exp(α+βx_i)/(1+exp(α+βx_i)))である。

モデル2:プロビットモデル

死亡数の分布は二項分布B(n_i, Φ((x_i-μ)/σ))である。但しΦは標準正規分布の分布関数で、Excelでは Φ((x-μ)/σ))=NORMDIST(x,μ,σ,1)で書けます。

モデル3:極値モデル

死亡数の分布は二項分布B(n_i, 1-exp(-exp(α+βx_i))である。

まずロジスティックモデルで確率を計算しましょう。

次に対数尤度の計算です。確率関数はp^x(1-p)^n-xでしたから、x=1ならば尤度はp、x=0ならば尤度は1-pです。x=1がx_i個、x=0がn_i-x_i個ですので、対数尤度はx_ilog p_i+(n_i-x_i)log (1-p_i)をiを1から8まで全部足したものになります。

F12セルの対数尤度を最大にするようなF1、F2セルの値を、ソルバーを用いて求めてください。
モデルによって推定した死亡確率と、実際の死亡割合を比較するとこのようになります。

点は実際の死亡割合、曲線がモデルです。

プロビットモデルに関しても同様にまず確率を求めて、それから対数尤度を計算してみましょう。μは最初はデータの中ごろということで1.8、σはとりあえず1にしてソルバーを使ってください。

すると、分母のσが0になってエラーになりました。
0にならないように変数の範囲に制約をつけましょう。

制約条件(U)の追加(A)をクリックします。

制約をつけるセルはH2、制約式に ＞0 と書きたいのですが、≧しか使えないのでとりあえず ≧0.1 と書きます。
OKをクリックすると思いがちですが、何故か追加(A)、キャンセルの順番でクリックしないと怒られます。

制約条件が追加されているのを確認したら実行します。

すると最適化は実行されましたが、σの値が先ほど書いた制約式の端っこである0.1になっています。
制約を変えると最適値も変わるかも知れません。

制約式を変えるために、先ほどの制約式をクリックして変更(C)をクリックします。
制約条件の0.1を0.01に書き換えて、先ほど同様追加(A)とキャンセルをクリックして、制約式が変わっているのを確かめてから実行します。

すると0.01では小さすぎてエラーになりました。

仕方ないので、今度は0.02にして実行してみましょう。

すると制約式の端っこではない最適値が求まりました。

最後に極値モデルに関しても最尤推定量を求めてください。

どのモデルが一番良いでしょうか？

レポート課題3

それぞれのモデルの最尤推定量を求めて、どのモデルが一番良いか選択しなさい。