rm(list=ls()); set.seed(20210323);
# Compare two Poisson rates
# Mr. October data
n0 <- 2820; y0 <- 563;
n1 <- 27; y1 <- 10;
library(rstan)
stan_model <- "
data{
int
int
int
int
}
parameters{
real
real
}
transformed parameters{
real
ratio <- theta1 / theta0;
}
model{
y0 ~ poisson(n0*theta0);
y1 ~ poisson(n1*theta1);
}
"
fit <- stan(model_code=stan_model,
data=c("n0", "y0", "n1", "y1"), iter=2000, chains=2)
round(summary(fit)$summary[,4:8], 2)
sims <- extract(fit)
names(sims)