a.al<-1
b.al<-1
a.be<-1
b.be<-1


Al<-rgamma(10000,a.al,b.al)
Be<-rgamma(10000,a.be,b.be)
hist(Al/(Al+Be))

mvec<-c(18,9,10,6,8,11,13,11)
yvec<-c(13,2,2,0,3,3,6,1)

n<-length(yvec)

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

jt.like<-function(yv,mv,thv,alv,bev){

	return(dbinom(yv,mv,thv,log=T)+dbeta(thv,alv,bev,log=T))

}

prior.func<-function(alv,bev,aa,ab,ba,bb){

	return(dgamma(alv,aa,ab,log=T)+dgamma(bev,ba,bb,log=T))
}

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

nburn<-20
nthin<-1
nsamp<-2000
nits<-nburn+nthin*nsamp

old.al<-old.be<-1
old.th<-(yvec+old.al)/(mvec+old.al+old.be)

old.likei<-jt.like(yvec,mvec,old.th,old.al,old.be)
old.like<-sum(old.likei)
old.prior<-prior.func(old.al,old.be,a.al,b.al,a.be,b.be)

ico<-0
theta.samp.1<-matrix(0,nrow=nsamp,ncol=n)
albe.samp.1<-matrix(0,nrow=nsamp,ncol=2)

for(iter in 1:nits){

	#Update the thetas in parallel
	new.th<-old.th+rnorm(n)*0.1
	new.th[new.th <= 0 | new.th >= 1]<-old.th[new.th <= 0 | new.th >= 1]
	new.likei<-jt.like(yvec,mvec,new.th,old.al,old.be)
	acc<-log(runif(n)) < new.likei-old.likei
	old.th[acc]<-new.th[acc]
	old.likei[acc]<-new.likei[acc]

	#Update alpha and beta using a Gibbs sampler
	new.al<-abs(rnorm(1,old.al,0.2))
	new.likei<-jt.like(yvec,mvec,old.th,new.al,old.be)
	if(log(runif(1)) < (sum(new.likei)+dgamma(new.al,a.al,b.al,log=T))-(sum(old.likei)+dgamma(old.al,a.al,b.al,log=T))){
		old.al<-new.al
		old.likei<-new.likei
	}

	new.be<-abs(rnorm(1,old.be,0.2))
	new.likei<-jt.like(yvec,mvec,old.th,old.al,new.be)
	if(log(runif(1)) < (sum(new.likei)+dgamma(new.be,a.be,b.be,log=T))-(sum(old.likei)+dgamma(old.be,a.be,b.be,log=T))){
		old.be<-new.be
		old.likei<-new.likei
	}

	

	if(iter > nburn & iter %% nthin == 0){
		ico<-ico+1
		theta.samp.1[ico,]<-old.th
		albe.samp.1[ico,]<-c(old.al,old.be)
	}

}

pdf("BB-Trace-1.pdf",paper="USr",height=9,width=11)
par(mfrow=c(1,2))
for(j in 1:n){
	plot(theta.samp.1[,j],ylab=substitute(theta[j],list(j=j)))
	acf(theta.samp.1[,j],main=substitute(theta[j],list(j=j)))
}
dev.off()

pdf("BB-Mu-1.pdf",paper="USr",height=9,width=11)
hist(albe.samp.1[,1]/apply(albe.samp.1,1,sum),breaks=seq(0,1,by=0.05),
	main=expression(paste("Posterior for ",mu)),xlab=expression(mu))
dev.off()

cvec<-c("black","red","blue","green","cyan","orange","purple","yellow")
plot(theta.samp.1[,1],pch=19,cex=.5,ylim=range(0,1))
for(j in 1:n){
	points(theta.samp.1[,j],pch=19,cex=0.5,col=cvec[j])
}
abline(h=mean(albe.samp.1[,1]/apply(albe.samp.1,1,sum)),lwd=2)

#######################################################################
#Compute, up to proportionality, the marginal posterior for (alpha,beta)

betabinom.like<-function(av,bv,yv,mv,aa,ab,ba,bb){

	alv<-rep(av,length(yv))
	bev<-rep(bv,length(yv))
	t1<-lgamma(alv+bev)-lgamma(alv)-lgamma(bev)	
	t2<-lgamma(mv+1)-lgamma(yv+1)-lgamma(mv-yv+1)	
	t3<-lgamma(mv+alv+bev)-lgamma(yv+alv)-lgamma(mv-yv+bev)
	t4<-dgamma(av,aa,ab,log=T)+dgamma(bv,ba,bb,log=T)	
	return(sum(t1+t2-t3)+t4)
}

a.al<-1
b.al<-1
a.be<-1
b.be<-1
betabinom.like(1,1,yvec,mvec,a.al,b.al,a.be,b.be)
al.vec<-be.vec<-seq(0.1,5.0,by=0.1)


bb.mat<-outer(al.vec,be.vec,"+")
for(i in 1:length(al.vec)){
	for(j in 1:length(be.vec)){
		bb.mat[i,j]<-betabinom.like(al.vec[i],be.vec[j],yvec,mvec,a.al,b.al,a.be,b.be)
	}
}

image(al.vec,be.vec,bb.mat,xlab=expression(alpha),ylab=expression(beta))
contour(al.vec,be.vec,bb.mat,add=T,levels=seq(-60,-20,by=1))
points(albe.samp.1,cex=0.5,pch=19)
title(main=expression(paste("Marginal posterior for ", "(",alpha, ",", phantom(), beta,")")),
sub = "Points drawn from full joint posterior")