Post 5: Refactoring Part I: a generic Metropolis-Hastings sampler

[This article was first published on Markov Chain Monte Carlo in Item Response Models, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

A best practice in software engineering is to re-use code instead of duplicating it. One reason for this is that it makes finding and fixing bugs easier. You only have to find the bug once instead of finding it every place you duplicated the code.

Since the person ability, item discrimination, and item difficulty samplers will all use the same Metropolis-Hastings (MH) algorithm, it would be useful to define each of those samplers in terms of a generic MH sampler.

In this post, we refactor the person ability sampler from Post 4 into: (1) a generic MH sampler, (2) a function to draw proposals for the person ability parameters, and (3) a function to calculate the complete conditional density qfor the person ability parameters. We then implement a new version of the person ability sampler and compare its output to the sampler from Post 4.

General refactoring strategy

The Metropolis-Hastings algorithm is roughly as follows:

  1. Propose a new value.
  2. Calculate the probability of accepting the new value.
  3. Flip a coin to determine if we accept the proposal or not.
  4. Update the state of the chain.

The sample.th code from the previous post implements these steps for the specific context of the person ability parameter. Below, we examine the code for each step of the algorithm, identify pieces that are specific to the person ability context, and then refactor the code into a general part and a specific part.

Step 1: Propose a new value

There are two places in the sample.th code that are specific to proposing new values in the person ability context. The first is in drawing the proposal values:

th.old    

The second is in calculating the last two terms of the acceptance probability:

 log.prop.star 

One way to combine these two tasks is to write a specific proposal function for the person ability parameters which will do two things (1) make the proposal draws and (2) return the log-density of the proposal distribution. For example:

## Proposal function for the person ability parameters
prop.th.abl 

Then the generic sampler would look something like this:

## A snippet of the generic sampler 
mh.sample 

where ... snip ... stands in for parts of the code which are not yet fleshed out.

Step 2: Calculating the acceptance probability

Now having drawn the proposals we can calculate the probability of accepting them. In the theta.abl sampler, this occurs in one line:

log.alpha.star 

where log.cc.star,log.cc.old, log.prop.star, and log.prop.old represent the value of complete conditional density at the proposal, the value of the complete conditional density at the old state of the chain, the value of the proposal density at the proposal, and the value of the “reversed” proposal density respectively.

Since the proposal function we defined above returns the proposal density, the only other information needed is the complete conditional density. We can implement this for the person ability parameter as:

## Complete conditional for the person ability parameters 
th.abl.cc 

Then we can use it as follows in the generic MH sampling function:

## A snippet of the generic sampler 
mh.sample 

Note that I have removed the call to pmin that was in the sample.th code. Details are in the next section.

Step 3: Flipping the coin

Conceptually, choosing which proposals are accepted is independent of the specific context. In the sample.th sampler, we implemented it so that it depended on P.persons:

 ## Accept or reject the proposals by flipping coins with the
 ## calculated acceptance probabilities 
 acc.new <- ifelse(log(runif(P.persons))<=log.alpha.star, 1, 0)

To make it general, we replace P.persons with the length of the log.alpha.star value:

    ## Step 3: Accept proposal draws
    ## ... Flip coins with the calculated probabilities  
    acc.new <- ifelse(
        log(runif( length(log.alpha.star) )) <= log.alpha.star,
        1, 0 )

which can be included in the generic mh.sample function immediately before the final ## ... snip ... line.

An aside about pmin

Recall that I removed the call to pmin that was in the sample.th code. In the chapter, the call to pmin paralleled the explanation of the MH algorithm that was given in the text. Specifically, the acceptance probability must be less than or equal to one. Since the fraction alpha.star can be greater than 1, we use the pmin function to set it equal to 1 in those cases. Or, on the log scale, the fraction log.alpha.star may be greater than zero, so we use pmin to set it equal to 0 in those cases.

Because we implement the coin flipping with the inverse CDF trick, it does not matter if log.alpha.star has been forced to be zero for all instances where its value is greater than zero. The inequality holds either way. Therefore we remove the call to pmin because it is unnecessary.

Step 4: Updating the chain

In the sample.th sampler, the updates were handled by referring to the th object within the cur object via the cur$th syntax:

 ## Update the chain 
 cur          

In R, an additional way to access the element of list is with the [[...]] syntax, where the ... is a variable that stores the name of the list element. For example, if param.name = 'th', then cur$th and cur[[param.name]] both reference the same object.

The [[...]] syntax can be useful when element names within a list are related to each other. For example, in the cur$ACC object there are two elements th and th.n. If param.name = 'th', then we can construct the th.n name with the paste function:

    name.n 

where the value of name.n will be th.n. Therefore, cur$ACC[[name.n]] will refer to cur$ACC$th.n and cur$ACC[[param.name]] will refer to cur$ACC$th.

See the next section for how these ideas are implemented in the generic MH sampler.

The final refactored code

The final generic MH sampler is then:

## A generic MH sampler
mh.sample 

The person ability proposal function is:

## Proposal function for the person ability parameters
prop.th.abl 

The person ability complete conditional density is:

## Complete conditional for the person ability parameters 
th.abl.cc 

The refactored person ability sampler is:

## Define the person ability sampler 
sample.th.refactor 

Testing the refactored code

In the previous post, we demonstrated that sample.th was implemented correctly. Therefore, to check that the refactored sampler is working correctly, it is sufficient to show that it returns identical output. To do so, we run a short chain with sample.th, run an otherwise identical chain with sample.th.refactor, and then compare their output.

## Load the necessary code from this and previous posts
source('http://mcmcinirt.stat.cmu.edu/setup/post-5.R')

## Check that we are using the Post 4
## version of sample.th
head(sample.th)
## 
## 1 function (U.data, old)                        
## 2 {                                             
## 3     th.old 

That the output is identical suggests that the code was refactored correctly.

To leave a comment for the author, please follow the link and comment on their blog: Markov Chain Monte Carlo in Item Response Models.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)