############################################################################################
#
# Example 3.6 (pages 102) – Figure 3.3 (page 103) - Weighted resampling algorithm
#
############################################################################################
#
# Multinomial data: Genetic linkage 
#
# References: Rao (1973) Linear statistical inference and applications. New York: Wiley
#             Tanner (1996, page 66)
#
##############################################################################################################
#
# Author : Hedibert Freitas Lopes
#          Graduate School of Business
#          University of Chicago
#          5807 South Woodlawn Avenue
#          Chicago, Illinois, 60637
#          Email : hlopes@ChicagoGSB.edu
#
###############################################################################################################

# Likelihood function
# -------------------
likelihood  =  function(th){return((2+th)^y[1]*(1-th)^y[2]*th^y[3])}

# Data, prior and likelihood
# --------------------------
y     = c(125,38,34)
th    = seq(0.4,0.8,length=1000)
like  = (2+th)^y[1]*(1-th)^y[2]*th^y[3]
like  = like/sum(like)

# valor da integral através da regra simples
h    = th[2]-th[1]
int  = h*sum(like)
like = like/int

# Weighted resampling
# -------------------
set.seed(92739)
M   = 10000
z   =  rnorm(M,0.63,0.05)
w   =  likelihood(z)/dnorm(z,0.63,0.05)
w   =  w/sum(w)
ind =  sample(1:M,size=M,prob=w,replace=T)
z1  =  z[ind]

# Figure 3.3
# ----------
par(mfrow=c(1,1))
xx = seq(min(z1),max(z1),length=40)
hist(z1,prob=T,xlab="",ylab="",main="",breaks=xx)
lines(th,like,lwd=2,lty=1)
lines(th,dnorm(th,0.63,0.05),lwd=2,lty=2)