Stan：ベルヌーイ分布（対数事後確率計算） - 機械学習・自然言語処理の勉強メモ

はじめに

前回は、ベルヌーイ分布を用いて事前確率のことを書いた。
kento1109.hatenablog.com

今回は、高度なモデルを作る際に避けては通れない「対数事後確率計算」についてまとめる。

対数事後確率計算

対数事後確率とは何なのか。
これはサンプル結果に出現する「lp__」のことだ。

> model
Inference for Stan model: bern.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

        mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
theta   0.71    0.00 0.11   0.48   0.64   0.71   0.78   0.89  1414    1
lp__  -10.80    0.02 0.78 -12.95 -10.92 -10.51 -10.35 -10.30  1629    1

以下、「アヒル本」の内容を引用しながら説明する。

まず、パラメータ $\theta$ の事後分布であるが、ベイズの定理を用いて以下のように表記できる。
$p(\theta|X)\propto p(X|\theta)p(\theta)$
$p(X|\theta)$ は尤度関数、 $p(\theta)$ は事前分布。

これを対数を取った対数事後確率は以下のように計算できる。
$\log p(\theta|X)\propto \log[p(X|\theta)p(\theta)]$
$\log p(\theta|X)\propto \log p(X|\theta)+\log p(\theta)$

次にある観測データ $X$ の確率分布を下記でモデル化する。
$Y[n]\sim{\rm Normal}(\mu,1)$
$\mu\sim{\rm Normal}(0,100)$

上記モデルの事後確率は以下の式で計算できる。
$\begin{align} p(\theta|X)\propto [\prod_{n=1}^N{\rm Normal}(\mu,1)] \times{\rm Normal}(\mu|0,1) \end{align}$
$\begin{align} \log p(\theta|X)\propto [\sum_{n=1}^N \log {\rm Normal}(\mu,1)]+\log {\rm Normal}(\mu|0,1) \end{align}$

標準的なStanコードは、

model{
  for (n in 1:N){
    Y[n] ~ normal(mu, 1);
  }
  mu ~ normal(0, 1); 
}

と書けたが、対数尤度を明示的に書くと、

model{
  for (n in 1:N){
    target += normal_lpdf(Y[n] | mu, 1);
  }
  target += normal_lpdf(mu | 0, 1); 
}

のような書き方となる。
normal_lpdf(Y[n] | mu, 1)は、 $\log {\rm Normal}(\mu,1)$ を意味する関数である。

簡単なモデルは、「～」の形で書けるが、高度なモデルではこの書き方が必要となる場合がある。

ベルヌーイ分布での計算

前回、尤度関数がベルヌーイ分布で表現されるモデルを以下のように書いた。

data { 
  int<lower=0> N;
  int<lower=0,upper=1> x[N];
} 
parameters {
  real<lower=0,upper=1> theta;
} 
model {
  for (n in 1:N)
    x[n] ~ bernoulli(theta);
}

明示的に尤度関数を書くと以下の記述となる。

data { 
  int<lower=0> N;
  int<lower=0,upper=1> x[N];
} 
parameters {
  real<lower=0,upper=1> theta;
} 
model {
  for (n in 1:N)
    target += bernoulli_lpmf(x[n] | theta);
}

結果を標準的な「～」と比べる。

１）x[n] ~ bernoulli(theta)

> model
Inference for Stan model: bern.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

       mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
theta   0.7    0.00 0.11   0.48   0.63   0.71   0.78   0.89  1411    1
lp__  -10.8    0.02 0.73 -12.90 -10.96 -10.52 -10.35 -10.30  1536    1

２）target += bernoulli_lpmf(x[n] | theta);

> model
Inference for Stan model: bern.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

        mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
theta   0.70    0.00 0.11   0.47   0.63   0.71   0.78   0.88  1107    1
lp__  -10.81    0.02 0.72 -12.80 -10.99 -10.54 -10.36 -10.30  1983    1

２つが同じ意味であることが分かる。

対数尤度関数を以下のように自分で記述しても問題ない。
$\log{\rm Bern}(x[n]|\theta)=x[n]\times\log(\theta)+(1-x[n])\times\log(1-\theta)$

data { 
  int<lower=0> N;
  int<lower=0,upper=1> x[N];
} 
parameters {
  real<lower=0,upper=1> theta;
} 
model {
  for (n in 1:N)
    target += x[n] * log(theta) + (1 - x[n]) * log(1-theta);
}

> model
Inference for Stan model: bern.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

        mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
theta   0.71    0.00 0.10   0.48   0.64   0.72   0.78   0.89  1533 1.01
lp__  -10.79    0.02 0.73 -12.88 -10.94 -10.50 -10.34 -10.30  1572 1.00

この場合も結果は同じとなる。