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 X_{i1} \sim Binomial(n_{i1}, \theta_{1}) and X_{i2} \sim Binomial(n_{i2}, \theta_{2}), conditionally independent given \theta_{1}, \theta_{2}, with n_{i1}, n_{i2} specified. Suppose further that the observed random variables are Y_{i} = X_{i1} + X_{i2}, i=1,2,3 so that the likelihood for \theta_{1}, \theta_{2} given Y_{1} = y_{1}, Y_{2} = y_{2}  and Y_{3} = y_{3} takes the form

\prod\limits_{i=1}^3 \left[ \sum\limits_{j_{i}} {n_{i1} \choose j_{k}}{n_{i2} \choose y_{i} - j_{i}}\theta_{1}^{n_{i1}-y_{i}}(1-\theta)^{n_{i1}-y_{i}}\theta_{2}^{y_{i}-j_{i}}(1-\theta_{2})^{n_{2}-y_{i}+j_{i}}\right]

Note that, at least in the version of the paper I have , the likelihood seems to have a typo. In particular, the exponent of (1-\theta_{1}) should be (n_{i1} - j_{i}) as above rather than (n_{i1} - y_{i}) 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 \theta_{1}, \theta_{2} 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
}

 

plot_posterior