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
共役事前分布(1)
ベルヌーイ分布の共役事前分布として知られる「ガンマ分布」を事前分布とする。
ただし、ハイパーパラメータ
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
共役事前分布(2)
次にガンマ分布のハイパーパラメータも最適化していく。
コードは以下の通り
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
観測データに合った情報を事前分布に設定することで、事後分布の分散が小さくなっている(精度が高くなっている)ことが分かる。