fit negative binomial to Poisson data

Pois_negbin.Rmd

---
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)
```