#Specify Data (Sufficient statistic for binary data):
x=15;n=35

#Set the support interval for theta:
a=1.0E-04; b=1-1.0E-04

par(mfrow=c(2,3),cex=0.75)
#############################################################################################

#Define likelihood:
likelihood=function(theta){
ifelse(theta<a || theta>b,0,(theta^x)*((1-theta)^(n-x)))}

#Define prior kernel:
#Zellner's prior:
prior.kernel=function(theta){
ifelse(theta<a || theta>b,0,(theta^theta)*((1-theta)^(1-theta)))}
#Jeffreys prior:
#prior.kernel=function(theta){dbeta(theta,0.5,0.5)}

#Compute posterior kernel:
post.kernel=function(theta){
sapply(theta,prior.kernel)*sapply(theta,likelihood)}

#Compute posterior density:
m=integrate(post.kernel,lower=0,upper=1)$value
post.density=function(theta){post.kernel(theta)/m}

#Plot likelihood posterior density
curve(post.density,from=0,to=1,xlab="theta",ylab="posterior density")


#Compute posterior estimates:
theta.mean=function(theta){theta*post.density(theta)}
post.mean=integrate(theta.mean,lower=0,upper=1)$value

theta2.mean=function(theta){(theta^2)*post.density(theta)}
post.var=integrate(theta2.mean,lower=0,upper=1)$value-post.mean^2
post.sd=sqrt(post.var)

#Function of parameters: computing posterior odds
odds=function(theta){(theta/(1-theta))*post.density(theta)}
post.odds=integrate(odds,lower=0,upper=1)$value

#Compute HPD intervals:
require(hdrcde)
theta.val=seq(a,b,l=1000)
hdr.den(den=list(x=theta.val,y=post.density(theta.val)),xlab="theta",ylab="posterior density")
#############################################################################################

#Generate samples from the posterior density:
#Set of number of iid samples to be generated:
N=1000

#Method-I: Use the method of Probability Integral Transform:

#Compute the posterior cdf:
post.cdf=function(theta){integrate(post.density,lower=a, upper=theta)$value}

#Compute the posterior inverse cdf:
post.icdf=function(u){
h=function(theta){post.cdf(theta)-u}
theta=uniroot(h,interval=c(a,b))$root
return(theta)}

#Generate MC samples from the posterior pdf:
set.seed(2357)
u=runif(N)
theta.sample=sapply(u,post.icdf)
hist(theta.sample,prob=T,col="lightblue",breaks=50,xlab="theta",
ylab="posterior density",main="PIF sampling")
theta.ordered=sort(theta.sample)
lines(theta.ordered,post.density(theta.ordered),col="blue")
post.mean1=mean(theta.sample)
post.sd1=sd(theta.sample)
mcse.mean1=post.sd1/sqrt(N)
post.odds1=mean(theta.sample/(1-theta.sample))
mcse.odds1=sd(theta.sample/(1-theta.sample))/sqrt(N)
hdr.den(theta.sample,xlab="theta",ylab="posterior density")
#############################################################################################

#Method-II: Use the method of Rejection sampling:
Max.val=optimize(post.kernel,interval=c(a,b),maximum=T)$objective
rej.sampling=function(N.try=500){
i=1
while(i<=N.try){theta.candidate=runif(1,a,b)
u=runif(1,0,Max.val)
if(u<post.kernel(theta.candidate)){return(theta.candidate)}
i=i+1}
}
theta.sample=sapply(rep(100,N),rej.sampling)
hist(theta.sample,prob=T,col="lightblue",breaks=50,xlab="theta",
ylab="posterior density",main="rejection sampling")
theta.ordered=sort(theta.sample)
lines(theta.ordered,post.density(theta.ordered),col="blue")
post.mean2=mean(theta.sample)
post.sd2=sd(theta.sample)
mcse.mean2=post.sd1/sqrt(N)
post.odds2=mean(theta.sample/(1-theta.sample))
mcse.odds2=sd(theta.sample/(1-theta.sample))/sqrt(N)
hdr.den(theta.sample,xlab="theta",ylab="posterior density")
#############################################################################################

#Summarize posterior estimates:
post.mean
post.mean1
mcse.mean1
post.mean2
mcse.mean2
post.sd
post.sd1
post.sd2
post.odds
post.odds1
mcse.odds1
post.odds2
mcse.odds2