######################################################
#M-H logit Example
######################################################

# Data
D <- c(1,1,1,0,0,0)
Z1<-c(0,3,0,0,0,2)
Z2<-c(0,1,0,0,1,0)
Z3<-c(0,1,0,2,0,0)
Z4<-c(1,-2,1,0,0,0)
X1<-c(1,0,-1,-1,0,1)
Z<-cbind(X1,Z1,Z2,Z3,Z4)
Z
D120<-rep(D,20)
length(D120)
Z120<-rbind(Z,Z,Z,Z,Z,Z,Z,Z,Z,Z)
dim(Z120)
Z120<-rbind(Z120,Z120)
dim(Z120)


beta.logit<-glm(D120~0+Z120,family=quasi(var="mu(1-mu)", link="logit"), start=c(0,0,0,0,0))
summary(beta.logit)
1-beta.logit$deviance/beta.logit$null.deviance	#pseudo-R2


#Functions for samples from the full conditional distributions
#function(y,X,mu,tau.sq,theta,proposal.sd)
sample.theta <- function(y,X,theta,proposal.sd){
                  theta.old <- theta
                  theta.new <- rnorm(ncol(X),theta.old,proposal.sd)
                      	pnew<-plogis(X %*% theta.new)
                      	pold<-plogis(X %*% theta.old)
                      alpha <- min(1,prod((pnew)^y*(1-pnew)^(1-y))/prod((pold)^y*(1-pold)^(1-y)))
                      if(is.nan(alpha)){
                      	print(alpha)
                      	theta<-theta.old
                      	accept<-0
                        return(list(theta=theta,accept=accept))
                      }else{
                      if(runif(1,0,1) < alpha){
                      		theta<-theta.new
                      		accept<-1
                            return(list(theta=theta,accept=accept))
                              }else{
                              	theta<-theta.old
                              	accept<-0
                                return(list(theta=theta,accept=accept))
                                    }
                                    }
}


#Set up MCMC
K <- 11000
accept <- NULL
proposal.sd <- .225

#Initialize parameters
theta.post <- rep(0,ncol(Z))

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

for(k in 1:K){
			 temp <- sample.theta(y=D120,X=Z120,theta.post,proposal.sd)
			 theta.post <- temp$theta
			 accept[k] <- temp$accept
			 theta.samples[k,]<-theta.post
}

samples.keep <- 1001:K

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

#Trace plots
par(mfrow=c(3,2),mar=c(1,2,1,1))
for(j in 1:ncol(Z)){
      plot(theta.samples[samples.keep,j],type='l')
      }

apply(theta.samples[samples.keep,],2,mean)
apply(theta.samples[samples.keep,],2,median)
apply(theta.samples[samples.keep,],2,sd)
apply(theta.samples[samples.keep,],2,quantile,c(.01,.025,.05,.1,.25,.75,.9,.95,.975,.99))


#use MCMCpack
library(MCMCpack)
mcmc.logit<-MCMClogit(D120~0+Z120, burnin = 1000, mcmc = 10000,
   tune=1.1)
summary(mcmc.logit)
#Trace plots
plot(mcmc.logit)