cookies=scan("cookies.txt")

#compare data historgram to the Poisson density
hist(cookies)
y=numeric(5)
z=numeric(5)
for(i in 1:5) {
 y[i]=5*i+10
 z[i]=ppois(y[i]+5,mean(cookies))-ppois(y[i],mean(cookies))
}
lines(y+2.5,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,9,1),type="l",ylim=c(0,.5))   
   #type="l" is a line, ylim sets the range on the y-axis
x2=seq(from=20,to=28,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*8e42)
#posterior
lines(x2,dgamma(x2,792,34),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,792,34)
qgamma(.975,792,34)