正規分布は、平均 μ と標準偏差 σ によって表される分布である。多くのデータが正規分布に従うために、よく利用される。ここで、平均 12.2、標準偏差 1.5 の正規分布から乱数を 100 個生成して、Stan でこれらの乱数が生成された正規分布のパラメーターを推定する例を示す。

x <- rnorm(n = 100, mean = 12.2, sd = 1.5)

Stan を利用して x を生成した正規分布の平均と標準偏差を推定する。Stan コードの data ブロックには、これから入力する値である x を定義する。Stan 自身は、x に何個のデータが入っているのかを計算できないので、外部からその個数を与える必要がある。このため、data ブロックには x の他に、x の要素数も合わせて定義する。次に、parameters ブロックには、これから推定したい正規分布のパラメーターである平均 μ と標準偏差 σ を定義する。最後に、model ブロックで、入力値変数とパラメーター変数を使って正規分布モデルを構築する。

data {
    int n;
    real x[n];
}
parameters {
    real mu;
    real sigma;
}
model {
    x ~ normal(mu, sigma);
}

R で Stan コードを呼び出して実行する。推定された平均と標準偏差はそれぞれ 12.17 および 1.59 であり、乱数を生成するパラメーターである 12.2 および 1.5 に非常に近いことがわかる。

d <- list(x = x, n = length(x))
fit <- stan(file = 'normdisp.stan', data = d)
fit
## Inference for Stan model: tmp.
## 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
## mu     12.17    0.00 0.16  11.85  12.06  12.17  12.27  12.53  3210    1
## sigma   1.59    0.00 0.11   1.47   1.51   1.59   1.67   1.96  2966    1
## lp__  -96.00    0.03 1.00 -98.64 -96.41 -95.70 -95.28 -95.01  1396    1

各回 MCMC サンプリングの結果をプロットする場合は、stan_trace 関数を利用し、パラメーターの事後確率のヒストグラムをプロットする場合は stan_hist 関数を利用する。

stan_trace(fit)
stan_hist(fit)

MCMC サンプリング結果は、そのままパラメーターが取りうる値とみなすことができる。そのため、α/2 分位点から 1 - α/2 分位点までの区間は、そのパラメーターにおける (1 - α)% のベイズ信用区間とみなすことができる。このことを利用して、平均と標準偏差の 95% 信頼区間を求めるには、次のように計算すればいい。

ms <- rstan::extract(fit)
quantile(ms$mu, probs = c(0.025, 0.975))
##     2.5%    97.5%
## 11.85349 12.53448

quantile(ms$sigma, probs = c(0.025, 0.975))
##     2.5%    97.5%
## 1.470064 1.959621