ポアソン分布のデータに負の二項分布をあてはめてみるテスト

ポアソン分布のデータに負の二項分布をむりやりあてはめてみたテストです。

R Notebookです。

---
title: "fit negative binomial to Poisson data"
output: html_notebook
---

## Data

```{r}
library(ggplot2)
library(magrittr)
library(rstan)
set.seed(20180902)

N <- 100
lambda <- 3
Y <- rpois(N, 3)

mean(Y)
var(Y)
```

```{r}
data.frame(Y) %>%
  ggplot() +
  geom_bar(aes(Y))
```

## Poisson

```{stan output.var=model1}
data {
  int<lower = 0> N;
  int<lower = 0> Y[N];
}
parameters {
  real<lower = 0> lambda;
}
model {
  Y ~ poisson(lambda);
}
```

```{r}
stan_data <- list(N = N, Y = Y)
fit1 <- sampling(model1, data = stan_data)
print(fit1)
```

## Negative binomial

```{stan output.var=model2}
data {
  int<lower = 0> N;
  int<lower = 0> Y[N];
}
parameters {
  real<lower = 0> alpha;
  real<lower = 0> beta;
}
model {
  Y ~ neg_binomial(alpha, beta);
}
```

```{r}
fit2 <- sampling(model2, data = stan_data)
print(fit2)
```

## Negative binomial (alternative parameterization)

```{stan output.var=model3}
data {
  int<lower = 0> N;
  int<lower = 0> Y[N];
}
parameters {
  real<lower = 0> mu;
  real<lower = 0> phi;
}
model {
  Y ~ neg_binomial_2(mu, phi);
}
```

```{r}
fit3 <- sampling(model3, data = stan_data)
print(fit3)
```

このようなサンプルとなりました。

サンプルの平均と分散です。このサンプルではやや分散がおおきくなっています。

[1] 3.07
[1] 3.439495

ポアソン分布にあてはめた結果です。当然、きっちりあてはめができました。

Inference for Stan model: stan-271669d441a5.
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  3.08    0.00 0.17  2.76  2.97  3.08  3.20  3.42  1446    1
lp__   37.99    0.02 0.70 36.11 37.82 38.26 38.43 38.48  1200    1

負の二項分布にあてはめた結果です。この程度の分散のおおきさではやはり識別不可能になるようです。

Inference for Stan model: stan-2716b014910.
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%
alpha 1.227166e+308     Inf  Inf 3.253547e+307 9.399872e+307 1.304465e+308 1.573187e+308
beta  4.020641e+307     Inf  Inf 1.010099e+307 3.063435e+307 4.248914e+307 5.126989e+307
lp__   1.454390e+03    0.07 1.16  1.451350e+03  1.453880e+03  1.454690e+03  1.455260e+03
              97.5% n_eff Rhat
alpha 1.776193e+308  4000  NaN
beta  6.023931e+307  4000  NaN
lp__   1.455710e+03   284    1

別のパラメタライゼーション(neg_binomial_2)の結果です。平均は推定できましたが、phiはやはり識別不可能でした。

Inference for Stan model: stan-27166ea4e7a8.
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
mu       3.09    0.03     0.18    2.76     2.96     3.08     3.21     3.43    35 1.13
phi  48705.30 2850.91 27822.69 4560.51 24091.45 46934.99 72696.75 97038.48    95 1.04
lp__    48.51    0.08     1.02   46.03    47.95    48.66    49.31    49.86   162 1.03

タグ：STAn