「データ解析のための統計モデリング入門」8章 ふらふら試行錯誤

「データ解析のための統計モデリング入門」8章 MCMCとメトロポリス法について。

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

• 作者: 久保 拓弥
• 出版社/メーカー: 岩波書店
• 発売日: 2012/05/19
• メディア: 単行本

まずは、対数尤度を計算するストリームを作る必要があります。

対数尤度の計算式は、

data * log(q) + (8 - data) * log(1 - q) + log(二項係数)

となるのですが、あらかじめ、二項係数を作っておく必要があります。

【二項係数】

# Rの場合は、choose(8, x)

if data = 0 then 1
elseif data = 1 then 8
elseif data = 2 then 28
elseif data = 3 then 56
elseif data = 4 then 70
elseif data = 5 then 56
elseif data = 6 then 28
elseif data = 7 then 8
elseif data = 8 then 1
else -1 endif

尤度を計算するノードとして、
function
data * log(q) + (8 - data) * log(1 - q) + log(二項係数)

function_new
data * log(q_new) + (8 - data) * log(1 - q_new) + log(二項係数)

各レコードの値を計算し、最後に集計をすれば、完成です。

完成したストリームはこちら。


スクリプトは、

var x
set x = 0.0
　
for I from 1 to 100
execute 'outputfile'
execute 'tablenode'
　
set x = value:tablenode.output at 1 1
　
set q.formula_expr = x
　
endfor

q = 0.30 とした場合、
-46.378
となります。

ふらふら試行錯誤のロジックとして、
1. ある初期値からスタートする。（例：q = 0.30）
2. 0.01、あるいは、-0.01だけ変化させる。
3. 対数尤度が改善されれば、その値を採用し、そうでない場合は、元の値を採用する。

これをストリームで書いて、後は、100回ほどループをさせれば、OK！

35回くらいで、最大値 0.46 に張り付いているのがわかります。
初期値を 0.30 にしても、0.60 にしても同じ結果になります。