ベイズ推計を手に入れた

ベイズ推定に使うメトロポリス法のプログラム、若干の修正が必要でしたがちゃんと動くようになりました。主な修正点は移動先の候補を

　X + (0.5 - Rnd) * 0.1

で求めるようにしたことと、移動先の確率密度が負になるときには0にするIF文を追加したことです。プログラムは次の通りで、これで事後分布の平均や95%信頼区間を求めることができます。

--------------------------------------------

Sub メトロポリス1218()

' メトロポリス法 Macro
' マクロ記録日 : 2010/12/18
'

'初期値設定と乱数初期化
  X = 0.5
  Randomize

'メインルーチン開始
  For T = 1 To 10000

'移動先候補
  Y = X + (0.5 - Rnd) * 0.1

'確率密度計算
  P = X ^ 250 * (1 - X) ^ 750
  Q = Y ^ 250 * (1 - Y) ^ 750
 
  If Q < 0 Then
   Q = 0
  End If

'移動判定
  If Q > P Then
   X = Y
  Else
   If Rnd < Q / P Then
     X = Y
   End If
End If

'出力ルーチン
    
   Sheets("sheet1").Cells(T + 1, 1).Value = X

  Next T
   
End Sub

--------------------------------------------

ためしに1000人中250人が「はい」と答えるという架空の調査結果に対して、「はい」と答える確率の事後分布を求めてみました。

個々の調査対象者が「はい」と答える確率がｘのとき、1000人中250人が「はい」と答える事象が発生する確率は

　　C(1000,250)*x^250*(1-x)^750

です。ここでC(1000,250)は1000人から250人を選ぶときの場合の数を表します。これが「はい」と答える確率がｘのときの尤度（ゆうど）となります。

事後分布は尤度と事前分布の積に比例するので、事前分布を一様分布とすると

　事後分布　∝　x^250*(1-x)^750

となります。上のメトロポリス法のプログラムによって事後分布に比例した確率でランダムウォークさせると次の図のようになりました。

これからヒストグラムをつくると次のようになります。

この結果から事後分布の平均を求めると25%で、調査対象者が「はい」と答える確率の95%信頼区間は、22%～28%であることがわかりました。それぞれ妥当な結果ですので、プログラムが期待通りの動きをしていると考えることができます。

これは簡単な例ですが、もっと複雑な式になる場合でも同様にパラメーターの平均や信頼区間を求めることができますので、一応ベイズ推定ができるようになったということができるでしょう。