Stan を利用したポアソン分布のパラメーター推定

ポアソン分布

ポアソン分布のパラメーター平均 λ だけからなる。Stan では poisson 関数を使ってモデルを構築することができる。ここで、λ = 20 であるポアソン分布から乱数を 100 個生成して、次に Stan で乱数が生成されたポアソン分布のパラメーターを推定する例を示す。

x <- rpois(n = 100, lambda = 20)

Stan を利用して x を生成したポアソン分布のパラメーター λ を求める。Stan コードの data ブロックには、これから入力する値である x を定義する。Stan 自身は、x に何個のデータが入っているのかを計算できないので、外部から与える必要がある。このため、x の他に、x の個数も data ブロックで定義する。次に、parameters ブロックには、これから推定したいポアソン分布のパラメーターである λ を定義する。最後に、model ブロックで、入力値変数とパラメーター変数を使ってポアソン分布モデルを構築する。

data {
    int n;
    int x[n];
}
parameters {
    real lambda;
}
model {
    x ~ poisson(lambda);
}

R では、次のようにして Stan コードを呼び出して実行する。推定結果を見ると、λ は 20.85 に推定され、実際に利用した 20 に近いことがわかる。

d <- list(x = x, n = length(x))
fit <- stan(file = 'poissondisp.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
## lambda   20.85    0.01 0.45   19.96   20.55   20.84   21.16   21.75  1624    1
## lp__   4241.31    0.02 0.71 4239.21 4241.14 4241.59 4241.76 4241.81  1686    1
## 
## Samples were drawn using NUTS(diag_e) at Wed Jan 10 15:19:26 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).

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

stan_trace(fit)
stan_hist(fit)
stan のポアソン分布のパラメーター推定結果(サンプリング軌跡)
stan のポアソン分布のパラメーター推定結果(事後確率の分布)