## Scottish lip cancer data (from WinBUGS):
O<-c(9,39,11,9,15,8,26,7,6,20,13,5,3,8,17,9,2,7,9,7,16,31,11,7,19,15,7,10,16,11,5,3,
     7,8,11,9,11,8,6,4,10,8,2,6,19,3,2,3,28,6,1,1,1,1,0,0)

E<-c(1.4,8.7,3.0,2.5,4.3,2.4,8.1,2.3,2.0,6.6,4.4,1.8,1.1,3.3,7.8,4.6,1.1,4.2,5.5,4.4,
     10.5,22.7,8.8,5.6,15.5,12.5,6.0,9.0,14.4,10.2,4.8,2.9,7.0,8.5,12.3,10.1,12.7,9.4,
     7.2,5.3,18.8,15.8,4.3,14.6,50.7,8.2,5.6,9.3,88.7,19.6,3.4,3.6,5.7,7.0,4.2,1.8)

#Summaries of the county-specific MLEs, leaving out county 1
#Model O[i]~Pois(E[i]*lambda[i])

lambda_hat<-O[-1]/E[-1]
mean(lambda_hat)
var(lambda_hat)
hist(lambda_hat,breaks=20)

#A gamma(1.58,1.11) seems to fit these rates fairly well
lambda<-seq(0,5,.01)
prior<-dgamma(lambda,1.58,1.11)
lines(lambda,8*prior/max(prior))


#Now let's analyze the data from County 1 with different priors:

#Prior 1, good guess based on the other counties
a1<-1.58;  b1<-1.11

#Prior 2, Ignores the other counties, huge variance
a2<-0.01;  b2<-0.01

#Prior 3, Strongly centered on the mean of the other counties
a3<-1.43*1.43/0.1;  b3<-1.43/0.1

#Plot these three priors and corresponding posteriors:

lambda<-seq(0,10,.01)

prior1<-dgamma(lambda,a1,b1)
prior2<-dgamma(lambda,a2,b2)
prior3<-dgamma(lambda,a3,b3)

post1<-dgamma(lambda,a1+O[1],b1+E[1])
post2<-dgamma(lambda,a2+O[1],b2+E[1])
post3<-dgamma(lambda,a3+O[1],b3+E[1])

plot(lambda,prior1,type="l",
     col=1,ylim=c(0,1),lty=1,
     xlab="Lambda",ylab="Density",
     main="O[1]/E[1]=6.4")
lines(lambda,prior2,lty=2)
lines(lambda,prior3,lty=3)
lines(lambda,post1,lty=1,lwd=2)
lines(lambda,post2,lty=2,lwd=2)
lines(lambda,post3,lty=3,lwd=2)

legend("topright",c("Prior 1", "Prior 2", "Prior 3", "Post 1", "Post 2", "Post 3"),
   inset=0.05,ncol=2,lwd=c(1,1,1,2,2,2),lty=c(1,2,3,1,2,3))


#Compute posterior summaries:
a<-c(a1,a2,a3)
b<-c(b1,b2,b3)

posterior.mean<-(a+O[1])/(b+E[1])
posterior.sd<-sqrt(a+O[1])/(b+E[1])
posterior.low<-qgamma(0.025,a+O[1],b+E[1])
posterior.high<-qgamma(0.975,a+O[1],b+E[1])