######################################################
#M-H Example
#Allergy Data
#2/2/06
######################################################

#Allergy Data
y <- c(17,12,22,35,14)
N <- c(29,18,28,45,30)
J <- 5

logit <- function(p) return(log(p/(1-p)))

#Functions for samples from the full conditional distributions
sample.theta.j <- function(y.j,N.j,mu,tau.sq,theta.j,proposal.sd){
                  theta.old <- theta.j
                  theta.new <- rnorm(1,theta.old,proposal.sd)
                  if(theta.new <= 0 | theta.new >= 1){
                      return(c(theta.old,0))
                      }else{
                      alpha <- (y.j*log(theta.new)+(N.j-y.j)*log(1-theta.new)-.5*(logit(theta.new)-mu)^2/tau.sq)-(y.j*log(theta.old)+(N.j-y.j)*log(1-theta.old)-.5*(logit(theta.old)-mu)^2/tau.sq)
                      if(log(runif(1,0,1)) < alpha){
                              return(c(theta.new,1))
                              }else{
                                    return(c(theta.old,0))
                                    }
                                    }
}

sample.mu <- function(theta,tau.sq,J)  return(rnorm(1,mean(logit(theta)),sqrt(tau.sq/J)))

sample.tau.sq <- function(theta,mu,J,alpha,beta){
                 alpha.post <- (J/2)+alpha
                 beta.post <- sum(.5*sum((logit(theta)-mu)^2)+beta)
                    return(1/rgamma(1,shape=alpha.post,scale=1/beta.post))
                    }

#Set up MCMC
K <- 10000
accept <- matrix(0,K,J)
proposal.sd <- .1

#Set hyperprior parameters
alpha <- .01
beta <- .01
                                     
#Initialize parameters
theta.post <- rep(.5,J)
mu.post <- 0
tau.sq.post <- 1

#Storage of samples
theta.samples <- matrix(0,K,J)
mu.samples <- rep(0,K)
tau.sq.samples <- rep(0,K)

for(k in 1:K){
      #print(k)
      for(j in 1:J){
            temp <- sample.theta.j(y[j],N[j],mu.post,tau.sq.post,theta.post[j],proposal.sd)
            theta.post[j] <- temp[1]
            accept[k,j] <- temp[2]
            }
      theta.samples[k,] <- theta.post      
      mu.post <- sample.mu(theta.post,tau.sq.post,J)
      mu.samples[k] <- mu.post
      tau.sq.post <- sample.tau.sq(theta.post,mu.post,J,alpha,beta)
      tau.sq.samples[k] <- tau.sq.post
}

samples.keep <- 1:K

#Compute the acceptance rate
accept.rate <- apply(accept[samples.keep,],2,mean)
accept.rate

#Trace plots
par(mfrow=c(J+2,1),mar=c(1,2,1,1))
for(j in 1:J){
      plot(theta.samples[samples.keep,j],type='l')
      }
plot(mu.samples[samples.keep],type='l')
plot(tau.sq.samples[samples.keep],type='l')      

#Boxplots of the thetas and inv.logit.mu
par(mfrow=c(1,1),mar=c(5,4,4,2))
boxplot(as.data.frame(cbind(theta.samples[samples.keep,],exp(mu.samples[samples.keep])/(1+exp(mu.samples[samples.keep])))),names=c(as.character(1:J),"inv.logit.mu"),main="Posterior of the thetas and the inv.logit.mu",ylim=c(0,1))
points(1:J,y/N,pch=19,col='green')