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