Stan で model にモデル式追加するのと increment_log_prob() で対数尤度関数追加するのは同じことって認識でいいのかな？
— hoxo_m (@hoxo_m) 2014, 7月 4
ということで、これを確認するために、以前書いた二番目に簡単な rstan コードを increment_log_prob() で書き直してみます。
書き直すのは model のこの部分。
model{
for(i in 1:N)
y[i] ~ normal(alpha + beta * x[i], s);
s ~ inv_gamma(0.001, 0.001);
}
これを次のように increment_log_prob() で対数尤度関数に書き換えてみます。
model{
for(i in 1:N)
increment_log_prob(-0.5 * log(2 * pi() * s^2) -0.5 * ((y[i] - (alpha + beta * x[i])))^2/s^2);
s ~ inv_gamma(0.001, 0.001);
}
set.seed(314) を最初に実行して、データをそろえた場合のモデル式バージョンと対数尤度関数バージョンの結果を比較してみます。
モデル式バージョンの結果
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
alpha 8.69 0.06 1.18 6.27 7.90 8.72 9.43 10.99 371 1.00
beta 0.83 0.00 0.02 0.78 0.81 0.83 0.84 0.88 373 1.00
s 7.01 0.01 0.17 6.66 6.91 7.01 7.12 7.33 478 1.01
lp__ -2449.55 0.07 1.34 -2453.09 -2450.15 -2449.17 -2448.58 -2448.10 362 1.03
対数尤度関数バージョンの結果
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
alpha 8.77 0.05 1.18 6.40 7.97 8.76 9.61 11.08 480 1.00
beta 0.83 0.00 0.02 0.78 0.81 0.83 0.84 0.87 477 1.01
s 7.02 0.01 0.15 6.71 6.92 7.01 7.12 7.31 593 1.01
lp__ -3368.42 0.05 1.17 -3371.37 -3369.02 -3368.13 -3367.55 -3367.04 527 1.00
ほぼ一致しました。
つまり Stan は尤度関数の形さえ分かってしまえば、どんなモデルでも自由に MCMC できるということですね。
この事実を知ったことで、自分の中でかなり Stan の理解が進んだように思えます。
追記
上記2つの結果で lp__ だけ大きく違いますが、これは Stan のモデル式が対数尤度の定数を計算しないためだと考えられます。
そこで、increment_log_prob() から定数を取り除いてみましょう。
model{
for(i in 1:N)
increment_log_prob(-0.5 * log(s^2) -0.5 * ((y[i] - (alpha + beta * x[i])))^2/s^2);
s ~ inv_gamma(0.001, 0.001);
}
定数項を除いた対数尤度関数バージョンの結果
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
alpha 8.83 0.05 1.11 6.66 8.09 8.86 9.56 10.93 413 1.01
beta 0.83 0.00 0.02 0.78 0.81 0.82 0.84 0.87 409 1.01
s 7.02 0.01 0.16 6.73 6.91 7.02 7.13 7.35 597 1.00
lp__ -2449.44 0.05 1.13 -2452.24 -2450.02 -2449.18 -2448.59 -2448.09 545 1.01
おー、lp__ もほぼ一致しましたね。
Stan はモデル式を対数尤度関数に変換していると考えると、理解がしやすいように思えます。
以上です。