Stanのlp__とは何なのか

うなどん

2017/12/29

はじめに

RとStanを用いてパラメータ推定を行った結果は、 次のように要約値として出力できます。

library(rstan)
rstan::summary(fit1)$summary
🍣              mean     se_mean        sd       2.5%        25%        50%
🍣 lambda   5.016936 0.005516983 0.2251207   4.599539   4.857984   5.012589
🍣 lp__   307.436659 0.017141353 0.6919190 305.580013 307.268628 307.701961
🍣               75%      97.5%    n_eff     Rhat
🍣 lambda   5.169179   5.468225 1665.051 1.001366
🍣 lp__   307.885998 307.941237 1629.372 1.002923

このとき、毎回必ず出力されてくる“lp__“という値。

今回はその謎について、ポアソン分布とその平均パラメータ\(\lambda\)の観点から迫っていきたいとおもいます。

パラメータ推定の実行

まずは、実際に推定を行います。

平均パラメータ\(\lambda=5\)のポアソン分布から乱数を100個生成してデータYとします。

#ポアソン分布からの乱数
set.seed(123)
Y <- rpois(100, 5)

#可視化
ggplot(data = data.frame(Y), mapping = aes(x = Y)) +
  geom_histogram(colour = "black", alpha = 0.5, bins = 11) +
  theme_bw()

Stanによる推定の実行

続いて、パッケージrstanを用いてポアソン分布の平均\(\lambda\)を推定します。

ファイル名は“lp__test.stan“とします。

//lp__test.stan
data {
  int N;
  int Y[N];
}

parameters {
  real<lower=0> lambda;
}

model {
  for (i in 1:N)
    Y[i] ~ poisson(lambda);
}

Rから推定を実行します。fit1にまずは結果を格納します。

#stanに渡すデータ
datastan = list(Y = Y, N = length(Y))

# ポアソン分布による平均の推定
model1 <- stan_model("lp__test.stan")
fit1 <- sampling(model1, data = datastan, seed = 123)

結果を確認してみましょう。rstanの関数summary()を使って要約した値を出力します。

#結果
round(rstan::summary(fit1)$summary,3)
🍣           mean se_mean    sd   2.5%     25%     50%     75%   97.5%
🍣 lambda   5.017   0.006 0.225   4.60   4.858   5.013   5.169   5.468
🍣 lp__   307.437   0.017 0.692 305.58 307.269 307.702 307.886 307.941
🍣           n_eff  Rhat
🍣 lambda 1665.051 1.001
🍣 lp__   1629.372 1.003

出ました。“lp__“です。平均が307.437。これは何や。

lp__と対数事後確率

“lp__“の正体は、MCMCの各ステップで提案されたパラメータの値(\(\lambda^\ast\))で計算された値です。log posteriorと呼ばれています。対数事後確率に近いものですが、”\(\sim\)確率分布()“を使ったモデル記述方法であるsampling statementを使うと、lp__に「定数項」なるものが加算されていません。

対数事後確率は、事後確率の対数を取ったもの。次の式で得ることができます。

\[ \log p(\lambda^\ast|Y) = \log p(Y|\lambda^\ast)+ \log p(\lambda^\ast)+constant \]

このとき気をつけなければならないのが「定数項(\(constant\))」の問題です。

“lp__“には(\(\sim\)確率分布を使った記法で推定した場合)パラメータに仮定した確率分布の対数尤度関数に含まれる、「定数項」が加味されません(つまり\(\log p(Y|\lambda^\ast)+ \log p(\lambda^\ast)\)がlp__となります)。Stanは推定したいパラメータで偏微分した値をサンプリングに利用しているため、「加算される定数項は、パラメータ推定上は加味しなくても良い」と表現するほうが正確かもしれません。“lp__“に何を足しても、Stanでパラメータ推定値が変化することはありません。たとえ🍣でもです。

ポアソン分布の尤度関数と対数尤度関数

「定数項」とは何か。ポアソン分布を例に見ていきます。ポアソン分布を例とする理由は、パラメータが一つでシンプルだからです。

ポアソン分布の確率質量関数は次のように表されます。確率質量関数を「尤度関数」と呼び変えても、ここでは差し支えありません。

\[ p(Y|\lambda) = \frac{\lambda^Ye^{-\lambda}}{Y!} \]

この尤度関数を対数尤度関数として表現すると、次になります。

\[ \log p(Y|\lambda) = Y \log \lambda-\sum_{k}^{Y}\log k-\lambda \]

このとき、推定したいパラメータ\(\lambda\)とは直接関係のない「定数項」が\(-\sum_{k}^{Y}\log k\)となります。

繰り返しとなりますが、StanではY[i] ~ poisson(lambda);のsampling statementで記述されたモデルを推定する際には、この定数項が省かれています。

これを、Stanで実際に確かめましょう。sanpling statementのかわりに、target記法を利用します。

定数項のないポアソン分布の対数尤度関数で推定

target記法とは、target+=を利用して、lp__を足し上げていく方法です。target記法を使えば、lp__に加算するものを自前で設定できます。Y[i] ~ poisson(lambda);をtarget記法で書き換えましょう。ポアソン分布の対数尤度関数を直接書きます。定数項\(-\sum_{k}^{Y}\log k\)は除きます。

ファイル名は“lp__test2.stan“とします。

//lp__test2.stan
data {
  int N;
  int Y[N];
}

parameters {
  real<lower=0> lambda;
}

model {
  for (i in 1:N)
    target += Y[i] * log(lambda) - lambda;
}

Rで推定を実行します。

#定数項の無い対数尤度関数を使ったモデル
model2 <- stan_model("lp__test2.stan")
fit2 <- sampling(model2, data = datastan, seed = 123)

結果が、こちらです。

#結果
round(rstan::summary(fit2)$summary,3)
🍣           mean se_mean    sd    2.5%     25%     50%     75%   97.5%
🍣 lambda   5.009   0.006 0.225   4.584   4.854   5.008   5.163   5.453
🍣 lp__   307.434   0.017 0.714 305.505 307.264 307.705 307.892 307.941
🍣           n_eff  Rhat
🍣 lambda 1536.609 1.002
🍣 lp__   1758.205 1.002

“lp__“の結果が一致していますね。

適当にlp__に値を足した場合

最後に、lp__に適当な値を加算して推定を実行してみます。

試しに、10000のsushiを足してみます。わざわざsushiに10000を代入しているのは、ただの遊び心です。特に意味はありません。

// lp__test3.stan
data {
  int N;
  int Y[N];
}

parameters {
  real<lower=0> lambda;
}

model {
  int sushi;
  sushi = 10000;
  for (i in 1:N)
    target += Y[i] * log(lambda) + sushi - lambda;
}

推定実行です。

#定数項としてsushiを加算した場合
model3 <- stan_model("lp__test3.stan")
fit3 <- sampling(model3, data = datastan, seed = 123)

結果がこちら。

#結果
round(rstan::summary(fit3)$summary,3)
🍣               mean se_mean    sd        2.5%         25%         50%
🍣 lambda       5.012   0.006 0.226       4.585       4.855       5.011
🍣 lp__   1000307.432   0.019 0.734 1000305.411 1000307.277 1000307.715
🍣                75%       97.5%    n_eff  Rhat
🍣 lambda       5.161       5.459 1495.023 1.002
🍣 lp__   1000307.892 1000307.941 1539.236 1.003

10000がデータ100個分だけlp__に加算されていますが、 \(\lambda\)パラメータの推定に影響は出ていません。

target記法でpoisson_lpmfを使った場合

まだ行きます。Stanでは、確率分布の対数尤度関数が定義された関数を使うことも出来ます。ポアソン分布の場合、poisson_lpmf()がそれです。ちゃんと定数項が入った場合のlp__、すなわち対数事後確率も計算しておきましょう。

// lp__test4.stan
data {
  int N;
  int Y[N];
}

parameters {
  real<lower=0> lambda;
}

model {
  for (i in 1:N)
    target += poisson_lpmf(Y[i] | lambda);
}

Rで実行です。

#定数項としてsushiを加算した場合
model4 <- stan_model("lp__test4.stan")
fit4 <- sampling(model4, data = datastan, seed = 123)

結果はこうなります。

#結果
round(rstan::summary(fit4)$summary,3)
🍣            mean se_mean    sd     2.5%      25%      50%      75%    97.5%
🍣 lambda    5.011   0.006 0.229    4.558    4.853    5.015    5.168    5.447
🍣 lp__   -215.275   0.020 0.732 -217.359 -215.462 -214.992 -214.800 -214.748
🍣           n_eff  Rhat
🍣 lambda 1389.894 1.002
🍣 lp__   1275.066 1.004

lp__の平均は-215あたりですね。

target+=で定数項のある対数尤度関数を書いて推定

最後です。poisson_lpmfの結果とlp__が一致するか、定数項のあるポアソン分布の対数尤度関数を書いて確認します。\(\sum_{k=1}^{Y} \log k\)\(\log(1) + \log(2)+ ... \log(Y)\)となりますので、for文を使って足し上げ、constantに格納してから定数項として組み入れます。

// lp__test5.stan
data {
  int N;
  int Y[N];
}

parameters {
  real<lower=0> lambda;
}

model {
  real constant;
  for (i in 1:N){
    constant = 0;
    for(k in 1 : Y[i]){
      constant = constant + log(k);
    }    
    target += Y[i] * log(lambda) - constant - lambda;
  }
}

Rでいざ実行です。

model5 <- stan_model("lp__test5.stan")
fit5 <- sampling(model5, data = datastan, seed = 123)

結果はこちら。

#結果
round(rstan::summary(fit5)$summary,3)
🍣            mean se_mean    sd     2.5%      25%      50%      75%    97.5%
🍣 lambda    5.017   0.006 0.224    4.583    4.865    5.021    5.166    5.460
🍣 lp__   -215.248   0.018 0.711 -217.242 -215.403 -214.974 -214.804 -214.748
🍣           n_eff  Rhat
🍣 lambda 1545.076 1.003
🍣 lp__   1596.317 1.000

-215で一致しました。

さいごに

stanでsampling statementを使ったモデル記述を行った場合、lp__には定数項の値が加算されていません。仮定する確率分布次第で定数項は変わりますので、lp__をモデル比較に用いる際には、とことん注意しましょう。

Enjoy!!