like = function(gamma,theta){
  sum(dbinom(0:3,3,gamma)*dbinom(8:5,9,theta))
}
pri = function(gamma,theta)
  dbeta(gamma,3,3)*dbeta(theta,2,2)
}
post = function(gamma,theta)
  pri(gamma,theta)*like(gamma,theta)
}

N          = 100
theta      = seq(0,1,length=N)
gamma      = seq(0,1,length=N)
prior      = matrix(0,N,N)
likelihood = matrix(0,N,N)
posterior  = matrix(0,N,N)
for (i in 1:N)
  for (j in 1:N){
    prior[i,j]= pri(gamma[i],theta[j])
    likelihood[i,j]= like(gamma[i],theta[j])
    posterior[i,j]= prior[i,j]*likelihood[i,j]
}
  
    
par(mfrow=c(1,3))
contour(gamma,theta,prior,xlab=expression(gamma),ylab=expression(theta))
title("Prior: gamma~Beta(3,3)- theta~Beta(2,2)")
contour(gamma,theta,likelihood,xlab=expression(gamma),ylab=expression(theta))
title("Likelihood: Pr(x1+x2=8|gamma,theta)\nx1~Bin(3,gamma) - x2~Bin(9,theta)")
contour(gamma,theta,posterior,xlab=expression(gamma),ylab=expression(theta))
title("Posterior: p(gamma,theta|x1+x2)=8")




# SIR (prior as proposal)
M = 10000
N = 10000
gammas.draw = rbeta(M,3,3)
thetas.draw = rbeta(M,2,2)
w = rep(0,M)
for (i in 1:M)
  w[i] = like(gammas.draw[i],thetas.draw[i])
k = sample(1:M,size=N,replace=TRUE,prob=w)
gammas.post = gammas.draw[k]
thetas.post = thetas.draw[k]

par(mfrow=c(1,2))
plot(gammas.draw,thetas.draw,xlab=expression(gamma),ylab=expression(theta),main="Prior draws",
xlim=c(0,1),ylim=c(0,1))
contour(gamma,theta,prior,add=TRUE,col=2,lwd=2)

plot(gammas.post,thetas.post,xlab=expression(gamma),ylab=expression(theta),main="Posterior draws",
xlim=c(0,1),ylim=c(0,1))
contour(gamma,theta,posterior,xlab=expression(gamma),ylab=expression(theta),add=TRUE,col=2,lwd=2)
abline(v=0.5,col=3,lwd=4)
abline(h=0.5,col=3,lwd=4)


# A few posterior quantities

mean((gammas.post<0.5)&(thetas.post<0.5))
mean((gammas.post<0.5)&(thetas.post>0.5))
mean((gammas.post>0.5)&(thetas.post<0.5))
mean((gammas.post>0.5)&(thetas.post>0.5))
c(mean(gammas.post),sqrt(var(gammas.post)),quantile(gammas.post,c(0.025,0.975)))
c(mean(thetas.post),sqrt(var(thetas.post)),quantile(thetas.post,c(0.025,0.975)))