va = read.table("va.txt",skip=9)
N.iter = 5000
N.burnin = 500
alpha = numeric(N.iter+N.burnin)
beta = numeric(N.iter+N.burnin)
n = length(va[,1])
theta = matrix(0,ncol=n,nrow=N.iter+N.burnin)
alpha[1] = 10                           # start at prior means
beta[1] = 10
theta[1,] = c(rep(.5,n))
accept.alpha = 0
accept.beta = 0
for(i in 1:(N.iter+N.burnin-1)) {
  log.theta.sum = sum(log(theta[i,]))
  log.theta.minus.sum = sum(log(1-theta[i,]))
  alpha.star = rnorm(1,alpha[i],sd=1)   
  if(alpha.star > 0) {
    prob = exp( n*lgamma(alpha.star+beta[i])-n*lgamma(alpha.star)
       +alpha.star*log.theta.sum-alpha.star/10 - (n*lgamma(alpha[i]+beta[i])
       -n*lgamma(alpha[i])+alpha[i]*log.theta.sum-alpha[i]/10) )
    if(runif(1) < prob) {
      alpha[i+1] = alpha.star
      accept.alpha = accept.alpha + 1
    }
    else alpha[i+1] = alpha[i]
  }
  else alpha[i+1] = alpha[i]

  beta.star = rnorm(1,beta[i],sd=1)   
  if(beta.star > 0) {
    prob = exp( n*lgamma(alpha[i+1]+beta.star)-n*lgamma(beta.star)
        +beta.star*log.theta.minus.sum-beta.star/10
        - (n*lgamma(alpha[i+1]+beta[i])-n*lgamma(beta[i])
           +beta[i]*log.theta.minus.sum-beta[i]/10) )
    if(runif(1) < prob) {
      beta[i+1] = beta.star
      accept.beta = accept.beta + 1
    }
    else beta[i+1] = beta[i]
  }
  else beta[i+1] = beta[i]

  #sample theta.j ~ Beta(y.j + alpha[i+1], n.j - y.j + beta[i+1])
  theta[i+1,] = rbeta(n, va[,1]+alpha[i+1], va[,2]-va[,1]+beta[i+1])
#  temp = rbeta(n, va[,1]+alpha[i+1], va[,2]-va[,1]+beta[i+1])
#  theta[i+1,] = temp
}
accept.alpha/(N.iter+N.burnin)
accept.beta/(N.iter+N.burnin)
par(mfrow=c(3,2))
plot(alpha,pch=".")
plot(beta,pch=".")
for(j in 1:4) plot(theta[,j],pch=".",ylab=paste("theta[",j,"]"))

alpha = alpha[(N.burnin+1):(N.iter+N.burnin)]
beta = beta[(N.burnin+1):(N.iter+N.burnin)]
theta = theta[(N.burnin+1):(N.iter+N.burnin),]
hist(alpha)
hist(beta)
for(j in 1:4) hist(theta[,j],main=paste("Histogram of theta[",j,"]"))

acf(alpha)
acf(beta)
for(j in 1:4) acf(theta[,j],main=paste("Series theta[",j,"]"))

mean(alpha)
var(alpha)
mean(beta)
var(beta)
mean(theta)

round(apply(theta,2,mean),3)       #separate means for each hospital (column)
round(va[,1]/va[,2],3)             #compare to MLE's

par(mfrow=c(1,1))
plot(va[,1]/va[,2],apply(theta,2,mean),xlab="MLE",ylab="posterior mean")
abline(0,1)

plot(va[,1]/va[,2],apply(theta,2,mean),col=1+(va[,2]<100))  
abline(0,1)               #color 1 is black, color 2 is red
abline(h=mean(theta),col="green")

#probability hospital 1 had a higher dissatisfaction rate than hospital 2
mean(theta[,1]>theta[,2])

#probability hospital 1 rates in the top half of dissatisfaction rates
theta.rank = apply(theta,1,rank)   #compute rank within rows (iterations)
dim(theta.rank)
mean(theta.rank[1,])
sd(theta.rank[1,])
mean(theta.rank[1,] < n/2)