Working out an example of a Bayesian weighted bootstrap problem
From time to time I like to go back and work a problem or example on a topic to refresh my memory, clarify the concept for myself and hopefully learn some new implementation techniques.
The following example comes from the well known paper “Bayesian Statistics Without Tears: A Sampling-Resampling Perspective” by Smith and Gelfand (1992). The problem appears earlier in the classic text Generalized Linear Models by McCullagh and Nelder, and Smith and Gelfand rework it in a Bayesian context using a weighted bootstrap to obtain an estimate of the parameter values.
The Problem
As directly described in Smith and Gelfand:
For i=1,2,3, suppose that and
, conditionally independent given
, with
specified. Suppose further that the observed random variables are
so that the likelihood for
given
and
takes the form
Note that, at least in the version of the paper I have , the likelihood seems to have a typo. In particular, the exponent of should be
as above rather than
as specified in the paper. This corrected version agrees with the original from the McCullagh and Nelder book.
The solution and R code for the implementation of the weighted bootstrap
The solution in the paper takes the joint prior for to be uniform over the unit square. The weighted bootstrap approach will sample from the prior distribution using a weighting proportional to the likelihood for each draw from the prior. So essentially, a point from the prior sample which has a likelihood of twice that of another will have twice the chance of being included in the posterior sample.
This code could be vectorized to be more efficient, but for our purposes works OK.
What follows is the R code for the calculation of the likelihoods and a plot of the posterior generated from sampling 2,000 posterior sample from an initial prior sample of 20,000.
The main code code will eventually call two functions which follow
#Populate the data from the problem
y<-c(7,5,6)
ni1<-c(5,6,4)
ni2<-c(5,4,6)
dta<-cbind(y,ni1,ni2)
# Populate the priors (these are the thetas)
U1<-runif(20000)
U2<-runif(20000)
for (i in 1:20000){
lkl[i] = calc_lkh(dta,U1[i],U2[i])
}
lkl_nrm=lkl/sum(lkl) # normalize likelihood
#
# Take a weighted bootstrap sample of indices using lkl as weight
bsamp=sample(c(1:20000),size=2000,replace=TRUE,prob=lkl_nrm)
plot(U1[bsamp],U2[bsamp],xlim=(0:1),ylim=(0:1),xlab=expression(paste(theta[1])),ylab=expression(paste(theta[2])),main="Sampled Posterior Distribution")
Function calc_lkh loops through calculating all components of the likelihood and the likelihood itself
function(p_dta,th1,th2){
lksum=c(0,0,0)
for (i in 1:3){
jmin=max(0,p_dta[i,1]-p_dta[i,3])
jmax=min(p_dta[i,2],p_dta[i,1])
for (j in jmin:jmax){
lk=lklhood(th1,th2,p_dta[i,1],p_dta[i,2],p_dta[i,3],j)
lksum[i]=lksum[i]+lk
}
}
lksum[1]*lksum[2]*lksum[3]
}
Function lklhood does the evaluation for each component in the sum of the stated likelihood
function(th1,th2,y,ni1,ni2,j){
chs1<-choose(ni1,j)
chs2<-choose(ni2,(y-j))
lpt1<-th1^j
lpt2<-(1-th1)^(ni1-j)
lpt3<-th2^(y-j)
lpt4<-(1-th2)^(ni2-y+j)
lk<-chs1*chs2*lpt1*lpt2*lpt3*lpt4
lk
}