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 --- changed from previous plot
beta1 = seq(0,25,length=B)
beta2 = seq(-.4,.0,length=B)

# contribution of 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])
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]))
}
}

# contribution of likelihood
y=c(1,1,1,1,0,0,0,0,0,0,0,0,1,1,0,0,0,1,0,0,0,0,0)
N=c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)
tau=c(53,57,58,63,66,67,67,67,68,69,70,70,70,70,72,
73,75,75,76,76,78,79,81)
mt=mean(tau)
n=length(y)

etaD=rep(1,B*n)
etaD=matrix(etaD,B,n)
for(j in 1:B)
{
for(k in 1:n)
{
etaD[,k] = beta1 + beta2[j]*(tau[k])
pi[,j] = pi[,j] * ((exp(etaD[,k])/(1+exp(etaD[,k])))**y[k])*
((1-(exp(etaD[,k])/(1+exp(etaD[,k]))))**(N[k]-y[k]))
}
}

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

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