機械学習・自然言語処理の勉強メモ

学んだことのメモやまとめ

Stan:ポアソン混合モデル

はじめに



今回は、アヒル本から離れて「ポアソン混合モデル」をモデリングする。
今回は前回の対数事後確率の計算が必要となる。
kento1109.hatenablog.com

ポアソン混合モデル



例えば、ある非負の値を取る観測データが以下のような分布になったとする。
f:id:kento1109:20180412174708p:plain
ヒストグラムから分布の「多峰性」が見て取れる。
(10あたりにピークを持つ分布と25あたりにピークを持つ分布)
これが「1つのパラメータから成るポアソン分布」から生成されたと考えるのは無理がありそう。
混合モデルでは、「観測データが生成された背景には、複数の分布が存在している」と考え、それをモデリングする。

観測データの生成

まずは、観測データを以下で生成する。

x <- c(rpois(50, 3), rpois(150, 10))
hist(x, breaks = seq(0,35,1) )

f:id:kento1109:20180413180448p:plain
ヒストグラムより多峰性が確認できる。
今回は真のパラメータ\lambda_1=3, \lambda_2=10と所属割合\pi_1=0.25, \pi_2=0.75を推定する。

モデリング



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

\thetaは無情報事前分布としている)

実行コードは以下の通り

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

両方ともそれなりに推定出来ていることが分かる。