In the chapter, we argue that a useful way to develop a Metropolis-Hastings (MH) within Gibbs sampler is to split the code into two levels. The top level is the “shell” of the MH within Gibbs algorithm, which sets up the data structures, implements burn-in and thinning of the chain, and calls a lower-level function to do the actual sampling from the complete conditionals. In this way, we can test the shell of the algorithm separately from testing the complete conditional samplers.
This post quotes the code presented in the chapter, discusses a data structure choice, gives an example of how to run the “empty” shell, and gives an example of how to visualize output from the shell.
MH within Gibbs 2PL shell
The code that we present in the chapter for the shell is:
## Define a complete run of the MH within Gibbs sampler run.chain.2pl
blocked.mcmc.update is defined with respect to the complete conditional samplers
## Define a single iteration of the MH within Gibbs sampler blocked.mcmc.update
which are, for now, defined as dummy samplers
## Define the dummy complete conditional sampler dummy.cc.sampler
Please see the chapter for details on the setup.
A technical aside about
Users of programming languages like C or C++ may be confused as to why we chose to pass around the omnibus object
cur, instead of passing the output matrix by reference to the complete conditional samplers and having them update it directly.
For technical and historical reasons, R implements a hybrid between pass by reference and pass by value. In R, an object is passed by reference if the function does not change the object, and passed value if it does. Therefore we use a small, omnibus object (
cur) to be copied about as the chain moves through a single iteration of
blocked.mcmc.update. This is much faster than passing the whole output matrix.
We could improve speed by mimicking pass by reference behavior with a package like R.oo, but at a large sacrifice to code readability. Therefore we define the
cur object to contain only the current state of the chain and the minimal additional information to calculate updates and monitor their success.
Running the shell
In Post 1 we discussed specific choices for the hyper priors. We can encode them in R as follows:
## Define the hyper parameters to values from Post 1 hyperpars
In Post 2 we generated fake data. We can run the sampler on the fake data as follows:
## Load the necessary code from this and previous posts source("http://mcmcinirt.stat.cmu.edu/setup/post-3.R") ## Set the seed to keep results reproducible set.seed(314159) ## Run the sampler at the true, simulated values ## for 1000 iterations run.A
Exploring the chain
run.A is a large matrix, printing it to the screen is not helpful. We can explore it as follows:
## run.A is a large matrix dim(run.A) ##  2061 1000 ## We can select a single parameter, by selecting a ## row of the matrix. For example: run.A['theta.abl 1',] ## selects the row for the first person ability parameter ## However, since there are 1000 columns, even one ## row is too much to print to the screen. Here ## we select columns 1 through 5 run.A['theta.abl 1', 1:5] ##  -1.682354 -1.682354 -1.682354 -1.682354 -1.682354 ## ## Which is more manageable. Note, that as expected, they ## are all equal. ## We can also define a vector of names first.names
Specifying a different vector for each set of outputs can be tedious. We define a function
get.2pl.params to build the selection vector:
## get.2pl.params takes ranges of parameters and appends the ## appropriate variable names. get.2pl.params
We test that the
get.2pl.params function is working by comparing its output to the
first.names variable we constructed earlier.
## N.B. get.2pl.params can give the same values as first.names all.equal( get.2pl.params( 1, 1, 1, 1 ), first.names) ##  TRUE ## So now we can look at the first five samples of the first ## three item discrimination parameters run.A[ get.2pl.params(NULL, 1:3, NULL, NULL), 1:5] ## [,1] [,2] [,3] [,4] [,5] ## a.disc 1 0.7120615 0.7120615 0.7120615 0.7120615 0.7120615 ## a.disc 2 1.1440699 1.1440699 1.1440699 1.1440699 1.1440699 ## a.disc 3 0.7326655 0.7326655 0.7326655 0.7326655 0.7326655
Visualizing sampler output
We can use coda to analyze and visualize the output of the sampler. Because we chose a matrix for the output, converting
coda‘s format is trivial:
require(coda) ## If the above command fails, you must install the ## coda library. To do so, un-comment the install.packages ## line below, run it, and follow the directions it gives. ## ## install.packages('coda') ## N.B. t(...) takes the transpose of a matrix run.A.mcmc
coda mcmc object behaves like a matrix. We can use the
get.2pl.params function to extract out the variables we want. Note that for the
mcmc object, the variables are in the columns.
run.A.mcmc[1:5, get.2pl.params(1,1,1,1)] ## theta.abl 1 a.disc 1 b.diff 1 sig2.theta ## [1,] -1.682354 0.7120615 -3 1.5625 ## [2,] -1.682354 0.7120615 -3 1.5625 ## [3,] -1.682354 0.7120615 -3 1.5625 ## [4,] -1.682354 0.7120615 -3 1.5625 ## [5,] -1.682354 0.7120615 -3 1.5625
plot( run.A.mcmc[, get.2pl.params(1,1,1,1)], density=FALSE )
for the help page, which gives details on how to customize the graph.
Note that the graph is not particularly impressive because the sampler is not doing anything yet.