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

# Grid domain
beta1 = seq(-5,5,length=B)
beta2 = seq(-.5,.3,length=B)

#prior density
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]-69.565)
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]))
}
}

pi=pi*1000#/(sum(pi)/(60*.8))

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