Stanによる、ゼロが観察されないポアソン分布のモデリング

ポアソン分布から発生したデータで、ゼロだけ観測できないという状況を想定します。そのようなデータについて元の平均とゼロを含めた観測数を推定します。

準備

グラフ表示用にggplot2と、StanによるMCMCのため今回はcmdstanrパッケージを使用します。

library(ggplot2)
library(cmdstanr)
## This is cmdstanr version 0.8.1
## - CmdStanR documentation and vignettes: mc-stan.org/cmdstanr
## - CmdStan path: /usr/local/cmdstan
## - CmdStan version: 2.36.0
register_knitr_engine(override = TRUE)

データ

データを生成します。観測数は100個で平均2のポアソン分布から生成します。

set.seed(123)
N <- 100
lambda <- 2

x <- rpois(N, lambda)
ggplot(data.frame(x = x), aes(x)) +
  geom_bar()

平均を表示します。

mean(x)
## [1] 2.02

データからゼロを除去します。

x_trunc <- x[x > 0]
ggplot(data.frame(x = x_trunc), aes(x)) +
  geom_bar() +
  xlim(0, NA)

残った観測値の数と平均を表示します。

length(x_trunc)
## [1] 87
mean(x_trunc)
## [1] 2.321839

モデル

Stanのモデルです。コンパイルしてstan_modelというオブジェクトに格納しておきます。

generated quantitiesブロックでは、ゼロになる確率p、ゼロにならなかった数Mから、負の二項分布の乱数によりゼロになった数を生成します。

data {
  int<lower=0> M;          // number of observations
  array[M] int<lower=1> X; // observations
}

parameters {
  real<lower=0> lambda;    // Poisson mean
}

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

generated quantities {
  int N;
  {
    real p = exp(-lambda);  // poisson_lpmf(0 | lambda)
    N = M + neg_binomial_rng(M, (1 - p) / p);
  }
}

あてはめ

CmdStanでモデルをあてはめます。

fit <- stan_model$sample(data = list(M = length(x_trunc),
                                     X = x_trunc),
                         chains = 4,
                         iter_warmup = 1000,
                         iter_sampling = 1000,
                         refresh = 0)
## Running MCMC with 4 sequential chains...
## 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.5 seconds.

結果

結果を表示します。データを生成した値に近い値が得られました。

fit$summary()
## # A tibble: 3 × 10
##   variable   mean median    sd   mad     q5    q95  rhat ess_bulk ess_tail
##   <chr>     <dbl>  <dbl> <dbl> <dbl>  <dbl>  <dbl> <dbl>    <dbl>    <dbl>
## 1 lp__     -21.1  -20.9  0.694 0.310 -22.5  -20.6   1.00    1985.    2555.
## 2 lambda     2.02   2.01 0.168 0.167   1.75   2.31  1.00    1633.    2223.
## 3 N        101.   100    4.78  4.45   94    109     1.00    2762.    3317.