**Markov Chain Monte Carlo in Item Response Models**, and kindly contributed to R-bloggers)

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

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

**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 on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...