Help us understand the problem. What is going on with this article?

Stan で increment_log_prob()

More than 5 years have passed since last update.

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

以上です。

hoxo_m
ホクソエム (hoxo_m) は架空のデータ分析者であり、日本の若手のデータ分析者集団のペンネームである。当初このデータ分析者集団は秘密結社として活動し、ホクソエムを一個人として活動させ続けた。
https://blog.hoxo-m.com/
hoxom
Machine Learning and Data Analysis Company for Your Smiles :)
http://hoxo-m.com/
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
Comments
No comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  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
ユーザーは見つかりませんでした