#Define the prior kernel:
#Zellner's prior:
cat("Zellner's prior:",fill=T)
k=function(theta){prod(theta^theta*(1-theta)^(1-theta))}
#Jeffrey's prior:
#cat("Jeffrey's prior:",fill=T)
#k=function(theta){prod(1/sqrt(theta*(1-theta)))}
#Laplace prior:
#cat("Laplace's prior:",fill=T)
#k=function(theta){prod(theta^0)}

#Enter sufficient statistics:
s=c(15,8);n=c(20,25)
#Define the likelihood function:
L=function(theta){prod(theta^s*(1-theta)^(n-s))}

#Define the posterior kernel:
K=function(theta){L(theta)*k(theta)}

#Compute the normalizing constant (denominator):
library(adapt)
denominator=adapt(ndim=2,functn=K,lower=rep(0,2),upper=rep(1,2))$value

#Define the parameter of interest:
eta=function(theta){theta[1]/theta[2]}
#eta=function(theta){theta[1]*(1-theta[2])/(theta[2]*(1-theta[1]))}

#Define the integrand for numerator and compute it:
num.fn=function(theta){eta(theta)*K(theta)}
numerator=adapt(ndim=2,functn=num.fn,lower=rep(0,2),upper=rep(1,2))$value

#Compute the posterior mean of eta:
post.mean.eta=numerator/denominator
cat(c("Posterior mean of eta:",round(post.mean.eta,4)),fill=T)

#Compute the maximum likelihood estimate (mle):
mle.theta=optim(par=c(0.5,0.5),fn=L,control=list(fnscale=-1))$par
mle.eta=eta(mle.theta)
cat(c("MLE of eta:",round(mle.eta,4)),fill=T)