[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.

In this post we refactor the proposal function from the previous post into a generic normal proposal function. This allows us to implement normal proposals for the (yet to be developed) samplers for the item parameters without duplicating code.

Our approach is to write a function which returns a function. We begin with a toy example to explain why that would work, implement a function that returns a normal proposal function, and then check that the refactored function works.

# The toy example

Writing functions that return functions sounds strange. Here we give a toy example to explain how it works.

Pretend that we want to write two similar functions for comparing numbers:

## A function to check if the parameter is less than three
is.three.less.than <- function( b.value ) {
return( 3 < b.value )
}

## A very similar function
is.two.less.than <- function( b.value ) {
return( 2 < b.value )
}

## Testing them out
c( is.three.less.than(2), is.three.less.than(4),
is.two.less.than(1), is.two.less.than(3) )
##  FALSE  TRUE FALSE  TRUE
##
## As expected.

R is designed to allow us to make a function which can build other functions. For example, we can make an "is X less than" function that takes X as a variable and returns a function like the is.three.less.than function above:

## We can write a function in R to build those other functions
is.X.less.than <- function( X ) {
## Return
return(
## a function
function( b.value ) {
## that depends on the value of X
return( X < b.value )
}
)
}

Now we can use is.X.less.than to define our functions from before, which still work the same:

## Define them
is.three.less.than <- is.X.less.than( 3 )
is.two.less.than <- is.X.less.than( 2 )

## Use them
c( is.three.less.than(2), is.three.less.than(4),
is.two.less.than(1), is.two.less.than(3) )
##  FALSE  TRUE FALSE  TRUE
##
## As expected.

### For the curious: Why does that even work?

This works because functions in R are implemented as a kind of paired object called a closure. A closure contains both the "text" of the function to be run and an "environment" which can contain additional variables.

In the example above, the is.X.less.than function returned a closure which contained both the "text" of function( b.value ) { return( X < b.value )} and the value of X. To see this for our example, we simply call the generated is.three.less.than function without parenthesis:

is.three.less.than
## function( b.value ) {
##              return( X < b.value )
##           }
## < environment: 0x3489ee8 >

And peek inside the environment with the ls(...) function, to see that it does contain X:

ls( envir=environment(is.three.less.than) )
##  "X"

And even see what X is set to in the environment of is.three.less.than:

environment(is.three.less.than)$X ##  3 Similarly, for the other function, we get environment(is.two.less.than)$X
##  2

as expected.

# A function to return a proposal function

Recall that in post 5 we defined the person ability proposal function to be:

## Proposal function for the person ability parameters
prop.th.abl <- function( state ) {

## Extract parameters from the state
## ... Previously:
##         th.old    <- old$th ## MH.th <- old$MH$th ## ... Now: th.old <- state$th
MH.th  <- state$MH$th

## Draw the proposal to update the state
## ... Previously:
##         P.persons <- length(th.old)
##         th.star   <- rnorm(P.persons,th.old,MH.th)
## ... Now:
state$th <- rnorm( length( th.old ), mean=th.old, sd=MH.th ) ## Build a function to return the density of the proposal ## ... Previously: ## log.prop.star <- log(dnorm(th.star,th.old,MH.th)) ## log.prop.old <- log(dnorm(th.old,th.star,MH.th)) ## ... Now, we define the function log.density <- function( state.to, state.from ) { dnorm( x = state.to$th,
mean = state.from$th, sd = MH.th, log = TRUE ) } ## Return the proposal object ## ... Note that we'll need the name of the parameter later, ## so return it too return( list( param.name='th', state=state, log.density=log.density )) } If we were okay with duplicating the code, we could copy the text of prop.th.abl and change all of the names from depending on th to depending on the name of the other parameter (e.g. a). Instead, we shall write a function which takes the name of the parameter as a variable and returns the equivalent function. Here we interleave comments for what the person ability version would have been: ## Write a function to return a proposal function for that parameter generic.normal.proposal <- function( param.name ) { return( function( state ) { ## Extract parameters from the state param.old <- state[[param.name]] ## th.old <- state$th

MH.param  <- state$MH[[param.name]] ## MH.th <- state$MH$th ## Draw the proposal to update the state state[[param.name]] <- rnorm( length( param.old ), mean=param.old, sd=MH.param ) ## state$th         <- rnorm(
##                       length( th.old ),
##                       mean=th.old,
##                        sd=MH.th )

## Build a function to return the density of the proposal
log.density <- function( state.to, state.from ) {
dnorm( x    = state.to[[param.name]],
##     state.to$th, mean = state.from[[param.name]], ## state.from$th,

sd   = MH.param,
##     MH.th,

log  = TRUE )
}

## Return the proposal object
return( list( param.name=param.name, state=state,
##      list( param.name='th', state=state,
##
log.density=log.density ))
}
)
}



Then implementing the refactored proposal function is simple:

prop.th.abl.refactor <- generic.normal.proposal('th')

# Testing the refactored code

Here we see if the sampler will produce the same output with the refactored code in place.

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

## Check that we have the version from Post 4.
## 1 function (state)
## 2 {
## 3     th.old <- state$th ## 4 MH.th <- state$MH$th ## 5 state$th <- rnorm(length(th.old), mean = th.old, sd = MH.th)
## 6     log.density <- function(state.to, state.from) {

## Set the seed and run the sampler for ten iterations
set.seed(314159)
run.post4 <- run.chain.2pl(
M.burnin = 0, M.keep  = 10, M.thin = 1,
U.data  = U, hyperpars = hyperpars,
th.init = runif( length(theta.abl),
min=-5, max=5     ),
a.init  = a.disc, b.init  = b.diff,
s2.init = sig2.theta,
MH.th=1, MH.a=NULL, MH.b=NULL)

## Swap in the refactored proposal
prop.th.abl.post4 <- prop.th.abl
prop.th.abl       <- prop.th.abl.refactor

## Set the seed and run the sampler for ten iterations
set.seed(314159)
run.refactor <- run.chain.2pl(
M.burnin = 0, M.keep  = 10, M.thin = 1,
U.data  = U, hyperpars = hyperpars,
th.init = runif( length(theta.abl),
min=-5, max=5     ),
a.init  = a.disc, b.init  = b.diff,
s2.init = sig2.theta,
MH.th=1, MH.a=NULL, MH.b=NULL)

## Check that they return the same values
all.equal( run.post4, run.refactor )
##  TRUE

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