y <- c(74, 99, 58, 70, 122, 77, 104, 129, 328, 119)
n <- length(y)
nsims <- 1000
thetasim <- rgamma(nsims,sum(y),n)
Tobs <- rep(0,nsims)
Trep <- rep(0,nsims)
for (i in (1:nsims)) {
  Tobs[i] <- sum((y-thetasim[i])^2/thetasim[i])
  yrep <- rpois(n,thetasim[i])
  Trep[i] <- sum((yrep - thetasim[i])^2/thetasim[i])
}
plot(Tobs,Trep)
sum(Trep > Tobs)
describe(Trep)

# Stan analysis of naive Poisson model
library(rstan)
traffic.data <- list(M=n, y=y)
stan_fit <- stan(file="H://HAL//Courses//Stat225//poisson.stan",data=traffic.data, iter=1000, chains=4)
print(stan_fit)
d <- as.data.frame(stan_fit)
sum(d$trep > d$tobs)
describe(d$trep)

# Stan analysi of Poisson-Gamma hierarchical model
traffic.data <- list(M=n, y=y)
stan_fit <- stan(file="H://HAL//Courses//Stat225//poissongamma.stan",data=traffic.data, iter=1000, chains=4)
print(stan_fit)
d <- as.data.frame(stan_fit)
sum(d$trep > d$tobs)/length(d$trep)
describe(d$tobs)
describe(d$trep)