## Week 5 Problems Solution

## Logistic Regression with Metropolis Hastings

library(Zelig)
data(turnout)


logit.mh <- function(y,X,beta.start,prior.mean=0,prior.var=10000,jump.var,n.sims=10000,burnin=0){
	
	library(mvtnorm)
	library(coda)

	k <- ncol(X)
	beta.cur <- t(beta.start)
	
	log.post.func <- function(beta,...){
		pi.i <- 1/(1+exp(-X %*% t(beta)))
		log.like <- sum(dbinom(y, size=1, prob=pi.i, log=T))
		log.prior <- dmvnorm(beta, mean=rep(prior.mean,k), sigma=diag(prior.var, k), log=T)
		log.post <- log.like + log.prior
		return(log.post)		
		}
	
	beta.update <- function(beta,...){
 		beta.cand <- rmvnorm(1, mean=beta, sigma=jump.var)		
 		r <- exp(log.post.func(beta.cand) - log.post.func(beta))
 		if(runif(1) <= r)
 			beta.cand
 		else beta
 		}
 		
 		
	beta.update2 <- function(beta,...){
		
		for(i in 1:k){
			beta.cand <- rnorm(1, mean=beta[i], sd=sqrt(jump.var[i]))		
			beta.old <- beta
			beta[i] <- beta.cand
			r <- exp(log.post.func(beta) - log.post.func(beta.old))
			if(runif(1) <= r)
				beta <- beta
			else beta <- beta.old
		}
		return(beta)
		}
		
	
	draws <- matrix(NA, nrow=n.sims+burnin, ncol=k)
	for(i in 1:(n.sims+burnin)){
		draws[i,] <- beta.cur <- beta.update(beta.cur)
		#draws[i,] <- beta.cur <- beta.update2(beta.cur)
		print(i)
		}
	res <- mcmc(draws[(burnin+1):(n.sims+burnin),])
	cat("Acceptance Rate:", 1-rejectionRate(res), "\n")
	return(res)
	}
	
	
y <- turnout$vote
X <- cbind(1,turnout$age, turnout$income)
mle <- glm(vote~age+income, data=turnout, family=binomial)
start.val <- c(0,0,0)

## part 1

jump.var <- vcov(mle)
system.time(posterior <- logit.mh(y=y,X=X,beta.start=start.val, jump.var=jump.var, n.sims=5000, burnin=500))
plot(posterior)

## part 2
## change beta.update to beta.update2

posterior2 <- logit.mh(y=y,X=X,beta.start=start.val, jump.var=c(.01,.0001,.001), n.sims=5000, burnin=500)
plot(posterior2)



check <- MCMClogit(vote~age+income, data=turnout)