Stan:ポアソン混合モデル
ポアソン混合モデル
例えば、ある非負の値を取る観測データが以下のような分布になったとする。
ヒストグラムから分布の「多峰性」が見て取れる。
(10あたりにピークを持つ分布と25あたりにピークを持つ分布)
これが「1つのパラメータから成るポアソン分布」から生成されたと考えるのは無理がありそう。
混合モデルでは、「観測データが生成された背景には、複数の分布が存在している」と考え、それをモデリングする。
観測データの生成
まずは、観測データを以下で生成する。
x <- c(rpois(50, 3), rpois(150, 10)) hist(x, breaks = seq(0,35,1) )
ヒストグラムより多峰性が確認できる。
今回は真のパラメータと所属割合を推定する。
モデリング
Stanでコードを書くと以下の通り。
data { int N; int K; int x[N]; } parameters { real theta[K]; simplex[K] pi; } model { real lp[K]; for(n in 1:N){ for (k in 1:K){ lp[k] = log(pi[k]) + poisson_lpmf(x[n] | theta[k]); } target += log_sum_exp(lp); } }
(は無情報事前分布としている)
実行コードは以下の通り
dat <- list(N=length(x), K=2, x=x) st_model <- stan_model(file="poisson_mixture.stan") fit <- sampling(st_model, data=dat, iter=4000, warmup=2000, thin=1, chain=2)
結果の確認
> fit Inference for Stan model: poisson_mixture. 2 chains, each with iter=4000; warmup=2000; thin=1; post-warmup draws per chain=2000, total post-warmup draws=4000. mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat theta[1] 3.47 0.01 0.48 2.55 3.14 3.46 3.79 4.45 1579 1 theta[2] 9.70 0.01 0.35 9.06 9.47 9.70 9.93 10.40 1708 1 pi[1] 0.26 0.00 0.05 0.17 0.22 0.26 0.30 0.37 1579 1 pi[2] 0.74 0.00 0.05 0.63 0.70 0.74 0.78 0.83 1579 1 lp__ -555.54 0.03 1.23 -558.72 -556.11 -555.22 -554.65 -554.14 1775 1
両方ともそれなりに推定出来ていることが分かる。