install.packages("MCMCpack", repos="http://lib.stat.cmu.edu/R/CRAN")
install.packages("coda", repos="http://lib.stat.cmu.edu/R/CRAN")
install.packages("VR", repos="http://lib.stat.cmu.edu/R/CRAN")
#Documentation at  http://mcmcpack.wustl.edu/files/MCMCpack.pdf

library(MCMCpack)

post.fun = function(theta) {
  return( theta[2]^0.5 * exp(-0.5*40/(.106^2)*(theta[1]-5.494)^2
            - 0.5*theta[2]*(theta[1]-5.5)^2 - theta[2]) )
}

post.samp = MCMCmetrop1R(post.fun,theta.init=c(5.5,1),mcmc=5000,logfun=F)
plot(post.samp)

#what we're used to seeing
par(mfrow=c(2,2))
plot(as.numeric(post.samp[,1]),pch=".",main="Mu")
hist(post.samp[,1])
plot(as.numeric(post.samp[,2]),pch=".",main="Tau")
hist(post.samp[,2])
par(mfrow=c(1,1))

summary(post.samp)
acf(post.samp)

#numerically usually better to use log-posteriors
logpost.fun = function(theta) {
  return( 0.5*log(theta[2]) - 0.5*40/(.106^2)*(theta[1]-5.494)^2
            - 0.5*theta[2]*(theta[1]-5.5)^2 - theta[2] )
}
post.samp = MCMCmetrop1R(logpost.fun,theta.init=c(5.5,1),mcmc=5000,logfun=T)
plot(post.samp)
acf(post.samp)

#tune proposal step size
post.samp = MCMCmetrop1R(logpost.fun,theta.init=c(5.5,1),mcmc=5000,tune=2)
plot(post.samp)
acf(post.samp)


#Linear regression
beer = read.table("beer.txt")
names(beer) = c("brand","price","calories","alcohol","type","domestic")
post.samp = MCMCregress(calories ~ alcohol, data=beer)
summary(post.samp)
plot(post.samp)

summary(lm(calories ~ alcohol, data=beer))
13.97^2 #residual variance
acf(post.samp)

attach(beer)
plot(alcohol,calories)

detach()
beer=beer[1:78,]
post.samp = MCMCregress(calories ~ alcohol, data=beer)
summary(post.samp)
plot(post.samp)


#second regression example
cereal = read.table("cereal.txt",header=T)
post.samp = MCMCregress(calories ~ protein + fat + carbo, data=cereal)
summary(post.samp)
plot(post.samp,ask=T)

summary(lm(calories ~ protein + fat + carbo, data=cereal))

#specify priors
post.samp = MCMCregress(calories ~ protein + fat + carbo, data=cereal, b0=c(0,4,9,4), B0=1)
summary(post.samp)
plot(post.samp,ask=T)


#logistic regression
post.samp = MCMClogit(domestic ~ price + calories + alcohol, data=beer)
summary(post.samp)
plot(post.samp)
acf(post.samp)

#do some ad-hoc variable selection
post.samp = MCMClogit(domestic ~ price + calories, data=beer)
summary(post.samp)
plot(post.samp)
acf(post.samp)