cookies=scan("cookies.txt")

#compare data histogram to the Poisson density
hist(cookies)
y=numeric(10)
z=numeric(10)
for(i in 1:10) {
 y[i]=2*i+12
 z[i]=ppois(y[i]+2,mean(cookies))-ppois(y[i],mean(cookies))
}
lines(y+1,z*length(cookies))

#plot prior, likelihood, and posterior on same graph
summary(cookies)
sum(cookies)
x=seq(from=.5,to=40,by=.5)
#prior
plot(x,dgamma(x,144/25 ,12/25 ),type="l",ylim=c(0,.7))   
   #type="l" is a line, ylim sets the range on the y-axis
x2=seq(from=17,to=25,by=.01)
#likelihood is a product of iid Poissons
likelihood=numeric(length(x2))
for(i in 1:length(x2)) likelihood[i] = prod(dpois(cookies,x2[i]))
#rescale so that it shows up on the plot
lines(x2,likelihood*6e70)
#posterior
lines(x2,dgamma(x2,1191+144/25,56+12/25),lty=2)  #lty=2 makes the dashed line

#compare 95% confidence intervals and credible intervals
mean(cookies)-1.96*sqrt(mean(cookies)/length(cookies))
mean(cookies)+1.96*sqrt(mean(cookies)/length(cookies))
qgamma(.025,1191+144/25,56+12/25)
qgamma(.975,1191+144/25,56+12/25)