7
Help us understand the problem. What are the problem?

More than 5 years have passed since last update.

posted at

updated at

Organization

Stan で increment_log_prob()

ということで、これを確認するために、以前書いた二番目に簡単な 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 はモデル式を対数尤度関数に変換していると考えると、理解がしやすいように思えます。

以上です。

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Sign upLogin
7
Help us understand the problem. What are the problem?