grads<-read.csv("http://madison.byu.edu/bayes/graduates.csv",header=T)
grads<-grads[!is.na(grads$Number.Of.Semesters),]
grads<-grads[grads$Act.Comp!=0,]
attach(grads)
xx<-Act.Comp
yy<-Number.Of.Semesters

# MCMC settings

length <- 2000
burn <- 500
candsig.b0 <- 0.01
candsig.b1 <- 0.0007
beta0 <- numeric(length+burn)
beta1 <- numeric(length+burn)

# starting values

beta0[1] <- 2.6
beta1[1] <- -.01

# priors

m.b0 <- 2.7
s.b0 <- 1
m.b1 <- 0
s.b1 <- .1


for(i in 2:(length+burn)){
# update for beta0
	beta0[i] <- beta0[i-1]
	old <- beta0[i-1]
	new <- rnorm(1,old,candsig.b0)
	llo <- -sum(exp(old+beta1[i-1]*xx))+
		sum(yy*(old+beta1[i-1]*xx))+
	       -((old-m.b0)^2)/(2*s.b0^2)
	lln <- -sum(exp(new+beta1[i-1]*xx))+
		sum(yy*(new+beta1[i-1]*xx))+
	       -((new-m.b0)^2)/(2*s.b0^2)
	uu<-runif(1,0,1)
	if(log(uu)<(lln-llo)){beta0[i]<-new}
# update for beta1
	beta1[i] <- beta1[i-1]
	old <- beta1[i-1]
	new <- rnorm(1,old,candsig.b1)
	llo <- -sum(exp(beta0[i]+old*xx))+
		sum(yy*(beta0[i]+old*xx))+
	       -((old-m.b1)^2)/(2*s.b1^2)
	lln <- -sum(exp(beta0[i]+new*xx))+
		sum(yy*(beta0[i]+new*xx))+
	       -((new-m.b1)^2)/(2*s.b1^2)
	uu<-runif(1,0,1)
	if(log(uu)<(lln-llo)){beta1[i]<-new}
}

par(mfrow=c(2,1))
plot(beta0)
plot(beta1)
good<-(burn+1):(length+burn)
library("MASS")
x11()
par(mfrow=c(1,1))
contour(kde2d(beta0[good],beta1[good]))
x11()
par(mfrow=c(2,1))
plot(density(exp(beta0[good])))
plot(density(beta1[good]))
x11()
par(mfrow=c(1,1))
plot(density(exp(beta0[good]+beta1[good]*20)),col="red",lwd=4,xlim=c(9,11))
lines(density(exp(beta0[good]+beta1[good]*25)),col="green",lwd=4)	
lines(density(exp(beta0[good]+beta1[good]*30)),col="blue",lwd=4)	
lines(density(exp(beta0[good]+beta1[good]*35)),col="purple",lwd=4)