set.seed(3122)
library(MASS)
library(coda)
n<-200
p<-3
Sig<-rWishart(1,p,diag(1,p))[,,1]
X<-cbind(rep(1,n),mvrnorm(n,mu=rep(0,p),Sigma=Sig))
p<-ncol(X)

be<-c(-2,1,0.5,-2.5)

muvec<-1/(1+exp(-X  %*%be))

Y<- rbinom(n,1,muvec)

lr.loglike<-function(Yv,Xv,bv){
	muvec<-1/(1+exp(-Xv  %*% bv))
	return(sum(dbinom(Yv,1,muvec,log=T)))
}

	

#####################################################
nburn<-20
nsamp<-2000
nthin<-1
nits<-nburn+nsamp*nthin
fg<-glm(Y~X-1,family=binomial)
beta.prec<-solve(summary(fg)$cov.unscaled)
beta.var<-solve(beta.prec)
beta.hat<-coef(fg)

#Prior variance
prior.sd<-10;prior.var<-prior.sd^2

##############################################################

#Algorithm 1

sigma.q<-sqrt(mean(diag(beta.var)))

old.beta<-beta.hat
old.like<-lr.loglike(Y,X,beta.hat)
old.prior<--0.5*sum(old.beta^2/prior.var)
old.term<-0

ico<-0
beta.samp1<-matrix(0,ncol=p,nrow=nsamp)
t0<-proc.time()[3]
for(iter in 1:nits){

	new.beta<-old.beta+rnorm(p)*sigma.q
	new.like<-lr.loglike(Y,X,new.beta)
	new.prior<--0.5*sum(new.beta^2/prior.var)
	new.term<-0

	if(log(runif(1)) < (new.like+new.prior-new.term)-(old.like+old.prior-old.term)){
		old.beta<-new.beta
		old.prior<-new.prior
		old.term<-new.term
	}
	if(iter > nburn & iter %% nthin == 0){
		ico<-ico+1
		beta.samp1[ico,]<-old.beta
	}

}
(t1<-proc.time()[3]-t0)

effectiveSize(beta.samp1)

plot(beta.samp1[,1])
pairs(beta.samp1)


pdf("Alg1-Trace.pdf",paper="USr",height=9,width=11)
par(mfrow=c(2,2))
for(j in 1:4){plot(beta.samp1[,j],xlab="Iteration",ylab=substitute(beta[jj], list(jj=j)))}
mtext("Algorithm 1: Joint Metropolis", side = 3, line = -2, outer = TRUE)
dev.off()
pdf("Alg1-ACF.pdf",paper="USr",height=9,width=11)
par(mfrow=c(2,2))
for(j in 1:4){acf(beta.samp1[,j],xlab="Iteration",main=substitute(beta[jj], list(jj=j)))}
mtext("Algorithm 1: Joint Metropolis", side = 3, line = -2, outer = TRUE)
dev.off()



##############################################################

#Algorithm 2


old.beta<-beta.hat
old.like<-lr.loglike(Y,X,beta.hat)
old.prior<--0.5*sum(old.beta^2/prior.var)
old.term<-(-0.5*t(old.beta-beta.hat) %*% beta.prec %*% (old.beta-beta.hat))[1,1]

ico<-0
beta.samp2<-matrix(0,ncol=p,nrow=nsamp)
t0<-proc.time()[3]
for(iter in 1:nits){

	new.beta<-mvrnorm(1,mu=beta.hat,Sigma=beta.var)
	new.like<-lr.loglike(Y,X,new.beta)
	new.prior<--0.5*sum(new.beta^2/prior.var)
	new.term<-(-0.5*t(new.beta-beta.hat) %*% beta.prec %*% (new.beta-beta.hat))[1,1]

	if(log(runif(1)) < (new.like+new.prior-new.term)-(old.like+old.prior-old.term)){
		old.beta<-new.beta
		old.prior<-new.prior
		old.term<-new.term
	}
	if(iter > nburn & iter %% nthin == 0){
		ico<-ico+1
		beta.samp2[ico,]<-old.beta
	}

}
(t1<-proc.time()[3]-t0)

effectiveSize(beta.samp2)

plot(beta.samp2[,1])
pairs(beta.samp2)


pdf("Alg2-Trace.pdf",paper="USr",height=9,width=11)
par(mfrow=c(2,2))
for(j in 1:4){plot(beta.samp2[,j],xlab="Iteration",ylab=substitute(beta[jj], list(jj=j)))}
mtext("Algorithm 2: Joint Metropolis-Hastings", side = 3, line = -2, outer = TRUE)
dev.off()
pdf("Alg2-ACF.pdf",paper="USr",height=9,width=11)
par(mfrow=c(2,2))
for(j in 1:4){acf(beta.samp2[,j],xlab="Iteration",main=substitute(beta[jj], list(jj=j)))}
mtext("Algorithm 2: Joint Metropolis-Hastings", side = 3, line = -2, outer = TRUE)
dev.off()

##############################################################

#Algorithm 3

sigma.qv<-sqrt(diag(beta.var))

old.beta<-beta.hat
old.like<-lr.loglike(Y,X,beta.hat)
old.prior<--0.5*sum(old.beta^2/prior.var)
old.term<-0

ico<-0
beta.samp3<-matrix(0,ncol=p,nrow=nsamp)
t0<-proc.time()[3]
for(iter in 1:nits){

	for(j in 1:p){
		new.beta<-old.beta
		new.beta[j]<-old.beta[j]+rnorm(1)*sigma.qv[j]
		new.like<-lr.loglike(Y,X,new.beta)
		new.prior<--0.5*sum(new.beta^2/prior.var)
		new.term<-0

		if(log(runif(1)) < (new.like+new.prior-new.term)-(old.like+old.prior-old.term)){
			old.beta<-new.beta
			old.prior<-new.prior
			old.term<-new.term
		}
	}
	if(iter > nburn & iter %% nthin == 0){
		ico<-ico+1
		beta.samp3[ico,]<-old.beta
	}

}
(t1<-proc.time()[3]-t0)

effectiveSize(beta.samp3)

plot(beta.samp3[,1])
pairs(beta.samp3)

pdf("Alg3-Trace.pdf",paper="USr",height=9,width=11)
par(mfrow=c(2,2))
for(j in 1:4){plot(beta.samp3[,j],xlab="Iteration",ylab=substitute(beta[jj], list(jj=j)))}
mtext("Algorithm 3: Metropolis-within-Gibbs", side = 3, line = -2, outer = TRUE)
dev.off()
pdf("Alg3-ACF.pdf",paper="USr",height=9,width=11)
par(mfrow=c(2,2))
for(j in 1:4){acf(beta.samp3[,j],xlab="Iteration",main=substitute(beta[jj], list(jj=j)))}
mtext("Algorithm 3: Metropolis-within-Gibbs", side = 3, line = -2, outer = TRUE)
dev.off()


##############################################################

#Algorithm 4
beta.prec.post<-solve(summary(fg)$cov.unscaled)+diag(rep(1/prior.var,p))
beta.var.post<-solve(beta.prec.post)
beta.hat.post<-beta.var.post %*% (solve(summary(fg)$cov.unscaled) %*% coef(fg))

old.beta<-beta.hat.post
old.like<-lr.loglike(Y,X,old.beta)
old.prior<--0.5*sum(old.beta^2/prior.var)
old.term<-(-0.5*t(old.beta-beta.hat.post) %*% beta.prec.post %*% (old.beta-beta.hat.post))[1,1]

ico<-0
beta.samp4<-matrix(0,ncol=p,nrow=nsamp)
t0<-proc.time()[3]
for(iter in 1:nits){

	new.beta<-mvrnorm(1,mu=beta.hat.post,Sigma=beta.var.post)
	new.like<-lr.loglike(Y,X,new.beta)
	new.prior<--0.5*sum(new.beta^2/prior.var)
	new.term<-(-0.5*t(new.beta-beta.hat.post) %*% beta.prec.post %*% (new.beta-beta.hat.post))[1,1]

	if(log(runif(1)) < (new.like+new.prior-new.term)-(old.like+old.prior-old.term)){
		old.beta<-new.beta
		old.prior<-new.prior
		old.term<-new.term
	}
	if(iter > nburn & iter %% nthin == 0){
		ico<-ico+1
		beta.samp4[ico,]<-old.beta
	}

}
(t1<-proc.time()[3]-t0)

effectiveSize(beta.samp4)

plot(beta.samp4[,1])
pairs(beta.samp4)

pdf("Alg4-Trace.pdf",paper="USr",height=9,width=11)
par(mfrow=c(2,2))
for(j in 1:4){plot(beta.samp4[,j],xlab="Iteration",ylab=substitute(beta[jj], list(jj=j)))}
mtext("Algorithm 4: Improved Joint Metropolis-Hastings", side = 3, line = -2, outer = TRUE)
dev.off()
pdf("Alg4-ACF.pdf",paper="USr",height=9,width=11)
par(mfrow=c(2,2))
for(j in 1:4){acf(beta.samp4[,j],xlab="Iteration",main=substitute(beta[jj], list(jj=j)))}
mtext("Algorithm 4: Improved Joint Metropolis-Hastings", side = 3, line = -2, outer = TRUE)
dev.off()


##############################################################

library(coda)
effectiveSize(beta.samp1)
effectiveSize(beta.samp2)
effectiveSize(beta.samp3)