rats <- read.table("~/BDA/rats.txt",header=T)
p.hats <- rats$y/rats$N

NNN <- 250

log.a.over.b.grid <- seq(from=-2.3,to=-1.3,length=NNN)
log.a.plus.b.grid <- seq(from=1,to=5,length=NNN)


post.a.b <- function(log.a.over.b,log.a.plus.b,ys,ns)
{
	beta <- exp(log.a.plus.b)/(1+exp(log.a.over.b))
	alpha <- beta*exp(log.a.over.b)
	J <- length(ys)
	log.prior <- -(5/2)*log.a.plus.b
	log.like <- 
		J*(lgamma(alpha+beta)-lgamma(alpha)-lgamma(beta))+
		sum(lgamma(alpha+ys)+lgamma(beta+ns-ys)-
		lgamma(alpha+beta+ns))
	log.posterior <- log.prior+log.like+738
	#  -738 is aproximately the value of the max log like
	post <- exp(log.posterior)
	return(post)
}

post.evals <- matrix(-999,NNN,NNN) # where to keep the posterior

for (i in 1:NNN)
{
	for (j in 1:NNN)
	{
		post.evals[i,j] <- 
			post.a.b(log.a.over.b.grid[i],log.a.plus.b.grid[j],rats$y,rats$N)
	}
}

# draw the contour plot
contour(log.a.over.b.grid,log.a.plus.b.grid,post.evals,xlab="log(alpha/beta)",ylab="log(alpha+beta)",
		levels=round(seq(from=.05,to=.95,length=5)*max(post.evals),1))

# draw the contour plot on alpha and beta scales
beta.grid <- exp(log.a.plus.b.grid)/(1+exp(log.a.over.b.grid))
alpha.grid <- beta.grid*exp(log.a.over.b.grid)

contour(alpha.grid,beta.grid,post.evals,xlab="alpha",ylab="beta",
		levels=round(seq(from=.05,to=.95,length=5)*max(post.evals),1))

post.log.a.over.b <- apply(post.evals,1,sum)
post.log.a.plus.b <- apply(post.evals,2,sum)
keep <- matrix(-999,1000,2)

quartz()
par(mfrow=c(1,2))
plot(log.a.over.b.grid,post.log.a.over.b,type="l")

plot(log.a.plus.b.grid,post.log.a.plus.b,type="l")


for (i in 1:1000)
{
	# Generate random alpha from marginal posterior of alpha.
	# Note: post.alpha doesn't add up to 1. The rmultinom() function
	# automatically normalizes.
	log.a.over.b.star <- 
		log.a.over.b.grid[which.max(as.vector(rmultinom(1,1,post.log.a.over.b)))]

	post.a.plus.b.given.y.l.a.over.b.star <- rep(-999,NNN)
	for (j in 1:NNN)
		post.a.plus.b.given.y.l.a.over.b.star[j] <- 
			post.a.b(log.a.over.b.star,log.a.plus.b.grid[j],rats$y,rats$N)
	
	# generate random beta given alpha.star and y
	log.a.plus.b.star <- 
		log.a.plus.b.grid[which.max(as.vector(rmultinom(1,1,post.a.plus.b.given.y.l.a.over.b.star)))]
	points(log.a.over.b.star,log.a.plus.b.star,pch=16,cex=.5)
	keep[i,] <- c(log.a.over.b.star,log.a.plus.b.star)

}

JJJ <- length(rats$y)
keep.thetas <- matrix(-999,1000,JJJ)

for (j in 1:1000)
{
	beta.star <- exp(keep[j,2])/(1+exp(keep[j,1]))
	alpha.star <- beta.star*exp(keep[j,1])
	keep.thetas[j,] <- rbeta(JJJ,alpha.star+rats$y,beta.star+rats$N-rats$y)
}

post.mean.thetas <- apply(keep.thetas,2,mean)
L.theta <- apply(keep.thetas,2,quantile,p=0.025)
U.theta <- apply(keep.thetas,2,quantile,p=0.975)


plot(p.hats,post.mean.thetas,type="n",xlim=c(0,.5),ylim=c(0,.5))
points(p.hats,post.mean.thetas,pch=16,cex=.5)
abline(c(0,1))
for (j in 1:JJJ)
	lines(c(p.hats[j],p.hats[j]),c(L.theta[j],U.theta[j]))

# compare to independent analysis
quartz()
post.mean.thetas.2 <- (rats$y)/(rats$N)
L.theta.2 <- qbeta(.025,(rats$y),(rats$N-rats$y))
U.theta.2 <- qbeta(.975,(rats$y),(rats$N-rats$y))
plot(p.hats,post.mean.thetas,type="n",xlim=c(0,.5),ylim=c(0,.5))
points(p.hats,post.mean.thetas.2,pch=16,cex=.5)
abline(c(0,1))
# (Note that when p.hat = 0, the posterior isn't defined.)
for (j in 1:JJJ)
	lines(c(p.hats[j],p.hats[j]),c(L.theta.2[j],U.theta.2[j]))