rm(list = ls())
# plot the density function pi()
B=1000  #size of square grid
# initialize pi as a B by B matrix of 1s
pi=rep(1,B*B)
pi=matrix(pi,B,B)

# Grid domain
beta1 = seq(-15,25,length=B)
beta2 = seq(-.4,.2,length=B)

# evaluate prior density on grid
rk=2  #rank of X
eta=rep(1,B*rk)
eta=matrix(eta,B,rk)
tautilde=c(55,75)
ytilde=c(1,.577)
Ntilde=c(1.577,1.577)

for(j in 1:B)
{
for(k in 1:rk)
{
eta[,k] = beta1 + beta2[j]*(tautilde[k])
pi[,j] = pi[,j] * ((exp(eta[,k])/(1+exp(eta[,k])))**ytilde[k])*
((1-(exp(eta[,k])/(1+exp(eta[,k]))))**(Ntilde[k]-ytilde[k]))
}
}

# rescale because we only care about
# shape of density
pi=pi*1000 #/(sum(pi)/(60*.8))

contour(beta1,beta2,pi,nlevels=8,
ylab=expression(beta["1"]),xlab=expression(beta[0]))