############################################################################################
#
# Independent bernoulli trials are colleted: y1,...,yn
#
# We model those trials via a Bernoulli model, i.e.
# yi iid Ber(theta) for theta in (0,1), so the likelihood function is
#
#      L(theta;s) = theta^s*(1-theta)^(n-s), 
#
# for s = y1 + y2 + ... + yn, the sum of successes.
#
# The prior for theta is a mixture of betas distributions, i.e.
#
#     theta ~ 0.1*Beta(1,1) + 0.9*Beta(8,2).
#
# The observed data is n=10 and s=2,i.e. 2 sucesses out of 10 trials.
#
############################################################################################
prior = function(theta){
  0.1*dbeta(theta,1,1)+0.9*dbeta(theta,8,2)
}
likelihood = function(theta){
 theta^3*(1-theta)^7
}
posterior = function(theta){
 likelihood(theta)*prior(theta)
}

pdf(file="binomial-mixturebetas-SIR-MH.pdf",width=12,height=9)
# Prior, likelihood and posterior densities
theta = seq(0,1,length=1000)
h = theta[2]-theta[1]
prio = prior(theta)
like = likelihood(theta)
like = like/sum(like)/h
post = posterior(theta)
post = post/sum(post)/h

par(mfrow=c(1,1))
plot(theta,post,xlab="theta",ylab="Density",type="l",xlim=c(0,1),main="",lwd=2,ylim=range(prio,like,post))
lines(theta,prio,col=2,lwd=2)
lines(theta,like,col=4,lwd=2)
legend("topleft",legend=c("Prior: 0.1*Beta(1,1)+0.9*Beta(8,2)","Likelihood: Beta(4,8) (3 successes in 10 trials)","Posterior"),col=c(2,4,1),lty=1,lwd=2,bty="n")

# SIR from U(0,1)
set.seed(2323466)
M = 10000
N = 5000
thetas = runif(M,0,1)
w = posterior(thetas)
w = w/sum(w)
thetas1 = sample(thetas,size=N,replace=TRUE,prob=w)
unique.sir = length(unique(thetas1))/M

par(mfrow=c(2,2))
breaks=seq(0,1,length=20)
plot(thetas,ylab="theta",xlab="Draws",main="Draws from proposal q(theta)=U(0,1)")
hist(thetas,xlab="theta",prob=TRUE,main="Draws from proposal q(theta)=U(0,1)",breaks=breaks);box()
abline(h=1,col=2,lwd=2)
plot(thetas,w,xlab="theta",ylab="weights",main="SIR",type="h",xlim=c(0,1))
hist(thetas1,xlab="theta",prob=TRUE,main="Posterior draws",xlim=c(0,1),breaks=breaks);box()
lines(theta,post,col=2,lwd=3)

sum.sir = c(unique.sir,mean(thetas1),var(thetas1),sqrt(var(thetas1)),quantile(thetas1,c(0.025,0.5,0.975)))
sum.sir

# Metropolis-Hastings algorithm
# The proposal is Beta(b*th.o/(1-th.o),b), so the expectation equals th.o
# It is easy to see that the variance is th.o*(1-th.o)^2/(b+1-th.o)
# Therefore, b works as a tuning parameter of the proposal distribution
mh.proposal = function(theta,thetastar,b){
  dbeta(thetastar,b*theta/(1-theta),b)
}

set.seed(121654)
b      = 20
theta0 = 1/2
theta  = rep(theta0,M)
for (i in 2:M){
  thetastar = rbeta(1,b*theta[i-1]/(1-theta[i-1]),b)
  num = posterior(thetastar)/mh.proposal(theta[i-1],thetastar,b)
  den = posterior(theta[i-1])/mh.proposal(thetastar,theta[i-1],b)
  alpha = min(1,num/den)
  if (runif(1)<alpha){
    theta[i]=thetastar
  }else{
    theta[i]=theta[i-1]
  }
}
unique.mh = length(unique(theta))/M
theta1 = theta[(M-N+1):M]

par(mfrow=c(1,1))
ts.plot(theta,xlab="iteration",main="MH draws\nProposal: Beta(b*theta/(1-theta),b), b=20",ylim=c(0,1))

sum.mh = c(unique.mh,mean(theta1),var(theta1),sqrt(var(theta1)),quantile(theta1,c(0.025,0.5,0.975)))

par(mfrow=c(1,2))
theta = seq(0,1,length=1000)
hist(thetas1,xlab="theta",prob=TRUE,main="SIR",xlim=c(0,1),breaks=breaks,ylim=c(0,2.5));box()
lines(theta,post,col=4,lwd=2)
hist(theta1,xlab="theta",prob=TRUE,main="MH",xlim=c(0,1),breaks=breaks,ylim=c(0,2.5));box()
lines(theta,post,col=4,lwd=2)

par(mfrow=c(1,1))
plot(density(thetas1),xlab="",ylab="Density",main="",lwd=2)
lines(density(theta1),col=2,lwd=2)
lines(theta,post,col=4,lwd=2)
legend("topleft",legend=c("SIR","MH","Exact"),col=c(1,2,4),lty=1,lwd=2)

summary = rbind(sum.sir,sum.mh)
rownames(summary) = c("SIR","MH")
colnames(summary) = c("Unique","Mean","Var","StDev","2.5%","50%","97.5%")
round(summary,4)

#    Unique   Mean    Var  StDev   2.5%    50%  97.5%
#SIR 0.3416 0.4417 0.0265 0.1628 0.1278 0.4538 0.7261
#MH  0.8193 0.4408 0.0266 0.1631 0.1354 0.4481 0.7315 

dev.off()