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) )
## [1] 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) )
## [1] 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) )
## [1] "X"

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

environment(is.three.less.than)$X
## [1] 3

Similarly, for the other function, we get

environment(is.two.less.than)$X
## [1] 2

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. 
head(prop.th.abl)                                                    
## 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 )
## [1] TRUE

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