Stan：ベルヌーイ分布（事前分布の導入）

はじめに

---

前回は単回帰分析を見ながら基本的なことを確認した。
kento1109.hatenablog.com

今回は、確率分布の問題としてなじみのある「ベルヌーイ分布」を見ていく。
また、事前分布の違いについても確認する。

無情報事前分布

---

まずは、以下のようにデータを作成する。

x <- rbinom(n = 15, size=1, prob = 0.7)
head(x)
[1] 1 0 0 1 1 1
mean(x)
[1] 0.7333333

目的は、データからベルヌーイ分布のパラメータを推定する。
Stanコードは以下の通り

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);
}

の事前分布は何も設定していない。
（設定しない場合、一様分布となる。）

以下でMCMC実行

model <- stan(file="bern.stan", data=list(N=length(x), x=x))

結果は以下の通り

> 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

平均はほぼ真の分布と一致している。

事後分布はこんな感じ。

stan_hist(model)

サンプルの分散は以下の通り

> ms <- rstan::extract(model)
> var(ms$theta)
[1] 0.01134403

共役事前分布（１）

---

ベルヌーイ分布の共役事前分布として知られる「ガンマ分布」を事前分布とする。
ただし、ハイパーパラメータa,bは固定する。
Stanコードは以下の通り

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

今回はずるいが、a=7,b=3を設定する。

Stanを回した結果は以下の通り

> model
> 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.72    0.00 0.09   0.53   0.66   0.72   0.78   0.87  1607    1
lp__  -15.32    0.02 0.69 -17.34 -15.50 -15.05 -14.88 -14.82  2009    1

stan_hist(model)

> ms <- rstan::extract(model)
> var(ms$theta)
[1] 0.007695371

共役事前分布（２）

---

次にガンマ分布のハイパーパラメータも最適化していく。
コードは以下の通り

data { 
  int<lower=0> N;
  int<lower=0> a;  // hyper parameter
  int<lower=0> b;  // hyper parameter
  int<lower=0,upper=1> x[N];
} 
parameters {
  real<lower=0,upper=1> theta;
} 
model {
  theta ~ beta(a + sum(x),b + N - sum(x));
  for (n in 1:N)
    x[n] ~ bernoulli(theta);
}

の更新式の導出に関しては下記参照。
qiita.com

MCMC実行

model <- stan(file="bern.stan", data=list(N=length(x), a=1, b=1, x=x))

実行結果

> 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.72    0.00 0.08   0.55   0.67   0.72   0.77   0.85  1486    1
lp__  -19.52    0.02 0.70 -21.38 -19.70 -19.25 -19.07 -19.01  1697    1

stan_hist(model)

> ms <- rstan::extract(model)
> var(ms$theta)
[1] 0.00612496

観測データに合った情報を事前分布に設定することで、事後分布の分散が小さくなっている（精度が高くなっている）ことが分かる。