#######################################################################
#
# Inspired in http://www.dme.ufrj.br/mcmc/chapter3.html
#
# Example 3.6 (pages 102)
# Figure 3.3 (page 103) 
# Weighted resampling algorithm (SIR) to simulate from
#
# pi(theta) propto (2+theta)^125*(1-theta)^38*theta^34
#
# for theta in (0,1).
#
#  Multinomial data: Genetic linkage 
# References: Rao (1973) Linear statistical inference 
#             and applications. New York: Wiley
#             Tanner (1996, page 66)
#
# Hedibert F. Lopes
# www.hedibert.org
# Date: January 28th 2021. 
#
#######################################################################

# Plotting posterior pi(theta) (up to a normalizing constant)
th = seq(0,1,length=1000)
par(mfrow=c(1,2))
plot(th,(2+th)^(125)*(1-th)^(38)*th^(34),type="l")

# Let's compute the normalizing constant (also known as
# prior predictive) via Monte Carlo integration
M     = 10000
a     = 0
b     = 1
theta = runif(M,a,b)
pred = mean((2+theta)^(125)*(1-theta)^(38)*theta^(34))
pred
post = (2+th)^(125)*(1-th)^(38)*th^(34)/pred

plot(th,post,ylab="Density",xlab=expression(theta),type="l")

# Computing E(theta)
I2 = pred
I1 = mean((2+theta)^(125)*(1-theta)^(38)*theta^(35))
E  = I1/I2

abline(v=E,col=2)

# Computing V(theta) = E(theta^2)-[E(theta)]^2
E2 = mean((2+theta)^(125)*(1-theta)^(38)*theta^(36))/I2
E2


V = E2-E^2
SD = sqrt(V)

abline(v=E+2*SD,col=2,lty=2)
abline(v=E-2*SD,col=2,lty=2)


# Computing Pr(theta>0.5) and Pr(theta>0.7)
ind = rep(0,M)
ind[theta<0.5]=1
lowertail = mean(ind*(2+theta)^(125)*(1-theta)^(38)*theta^(34))/I2
lowertail

ind = rep(0,M)
ind[theta>0.7]=1
uppertail = mean(ind*(2+theta)^(125)*(1-theta)^(38)*theta^(34))/I2
uppertail


# What if I want to draw from this posterior?
# simplest solution: Sampling importance resampling (SIR)
M     = 100000

# SIR step 1: sampling from q(.)=U(0,1)
theta = runif(M)
par(mfrow=c(1,2))
hist(theta,main="Draws from proposal")

# SIR step 2: computing weights (unnormalized)
numerator = (2+theta)^(125)*(1-theta)^(38)*theta^(34)
denominator = 1
w = numerator/denominator
w = w/sum(w)

# SIR step 3: resampling the set {theta[1],...,theta[M]} with 
# with weights {w[1],...,w[M]}
theta1 = sample(theta,replace=TRUE,size=M,prob=w)
hist(theta1,prob=TRUE,main="Draws from target=posterior")
lines(th,post,col=2)


# Comparing MC integration with SIR
MCI = round(c(E,V,SD,lowertail,uppertail),6)  # MC integration approximation
SIR = round(c(mean(theta1),var(theta1),sqrt(var(theta1)),mean(theta1<0.5),
mean(theta1>0.7)),6) # SIR approxiamtion 

tab = rbind(MCI,SIR)
colnames(tab) = c("Mean","Var","StDev","Pr<0.5","Pr>0.7")
tab


# Distribution of gamma=log(theta/(1-theta))
gamma = log(theta1/(1-theta1))
par(mfrow=c(1,1))
hist(gamma)
tails = quantile(gamma,c(0.025,0.5,0.975))
abline(v=tails[1],col=2,lty=2)
abline(v=tails[2],col=2)
abline(v=tails[3],col=2,lty=2)
abline(v=mean(gamma),col=4,lwd=2)