ゼロがなくなったポアソン分布のパラメーターを切断ポアソン分布で推定する

先日のモデルですが、Stanでは切断ポアソン分布をつかえるので、それをつかえばよいのでした。

マニュアル13.7節を参考にしました。

data {
  int<lower = 0> M;
  int<lower = 1> X[M];
}

parameters {
  real<lower = 0> lambda;
}

model {
  for (m in 1:M) {
    X[m] ~ poisson(lambda) T[1, ];
  }
}

generated quantities {
  int N;
  {
    real p = poisson_cdf(0, lambda);  // 失敗確率
    N = M + neg_binomial_rng(M, (1 - p) / p);
  }
}

Rコードです。データは前のエントリと同じです。

fit2 <- stan("truncPois.stan",
             data = list(M = length(x2), X = x2))

結果です。

Inference for Stan model: truncPois.
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   2.03    0.00 0.17   1.70   1.91   2.02   2.14   2.39  1366    1
N      100.56    0.09 4.69  93.00  97.00 100.00 103.00 111.00  2833    1
lp__   -21.14    0.02 0.72 -23.26 -21.30 -20.86 -20.68 -20.63  1850    1

タグ：STAn R