File: debug.Rnw

package info (click to toggle)
r-cran-mcmc 0.9-7-1
• area: main
• in suites: bullseye, sid
• size: 2,352 kB
• sloc: ansic: 1,298; makefile: 14; sh: 8
 file content (274 lines) | stat: -rw-r--r-- 13,032 bytes parent folder | download | duplicates (6)
 123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274 \documentclass{article} \usepackage{amstext} % \VignetteIndexEntry{Debugging MCMC Code} \begin{document} <>= options(keep.source = TRUE, width = 60) foo <- packageDescription("mcmc") @ \title{Debugging MCMC Code} \author{Charles J. Geyer} \maketitle \section{Introduction} This document discusses debugging Markov chain Monte Carlo code using the R contributed package \texttt{mcmc} (Version \Sexpr{foo$Version}) for examples. It also documents the debugging output of the functions \texttt{mcmc} and \texttt{temper}. Debugging MCMC code if the code is taken as a black box is basically impossible. In interesting examples, the only thing one knows about the equilibrium distribution of an MCMC sampler is what one learns from the samples. This obviously doesn't help with testing. If the sampler is buggy, then the only thing you know about the equilibrium distribution is wrong, but if you don't know it is buggy, then you don't know it is wrong. So you don't know anything. There is no way to tell whether random output has the correct distribution when you don't know anything about the distribution it is supposed to have. The secret to debugging MCMC code lies in two principles: \begin{itemize} \item take the randomness out, and \item expose the innards. \end{itemize} The first slogan means consider the algorithm a deterministic function of the elementary pseudo-random numbers that are trusted (for example, the outputs of the R random number generators, which you aren't responsible for debugging and are also well tested). The second slogan means output, at least for debugging purposes, enough intermediate state so that testing is straightforward. For a Gibbs sampler, this means outputting all of the trusted elementary pseudo-random numbers used, the state before and after each elementary Gibbs update, and which update is being done if a random scan is used. Also one needs to output the initial seeds of the pseudo-random number generator (this is true for all situations and will not be mentioned again). For a Metropolis-Hastings sampler, this means outputting all of the trusted elementary pseudo-random numbers used, the state before and after each elementary Metropolis-Hastings update, the proposal for that update, the Hastings ratio for that update, decision (accept or reject) in that update. For more complicated MCMC samplers, there is more innards'' to expose'' (see the discussion of the \texttt{temper} function below), but you get the idea. You can't output too much debugging information. \section{The Metrop Function} The R function \texttt{metrop} in the \texttt{mcmc} package has an argument \verb@debug = FALSE@ that when \verb@TRUE@ causes extra debugging information to be output. Let \texttt{niter} be the number of iterations \verb@nbatch * blen * nspac@, and let \texttt{d} be the dimension of the state vector. The result of invoking \texttt{metrop} is a list. When \verb@debug = TRUE@ it has the following additional components \begin{itemize} \item \texttt{current}, an \texttt{niter} by \texttt{d} matrix of mode \verb@"numeric"@, the state before iteration \texttt{i} is \verb@current[i, ]@ \item \texttt{proposal}, an \texttt{niter} by \texttt{d} matrix of mode \verb@"numeric"@, the proposal for iteration \texttt{i} is \verb@proposal[i, ]@ \item \texttt{z}, an \texttt{niter} by \texttt{d} matrix of mode \verb@"numeric"@, the vector of standard normal random variates used to generate the proposal for iteration \texttt{i} is \verb@z[i, ]@ \item \texttt{log.green}, a vector of length \texttt{niter} and mode \verb@"numeric"@, the logarithm of the Hastings ratio for each iteration \item \texttt{u}, a vector of length \texttt{niter} and mode \verb@"numeric"@, the$\text{Uniform}(0, 1)$random variate compared to the Hastings ratio for each iteration or \texttt{NA} if none is needed (when the log Hastings ratio is nonnegative) \item \texttt{debug.accept}, a vector of length \texttt{niter} and mode \verb@"logical"@, the decision for each iteration, accept the proposal (\texttt{TRUE}) or reject it (\texttt{FALSE}) \end{itemize} (The components \texttt{z} and \texttt{debug.accept} were added in version 0.7-3 of the \texttt{mcmc} package. Before that only the others were output.) Two components of the list returned by the \texttt{metrop} function always (whether \verb@debug = TRUE@ or \verb@debug = FALSE@) are also necessary for debugging. They are \begin{itemize} \item \texttt{initial.seed} the value of the variable \texttt{.Random.seed} that contains the seeds of the R random number generator system before invocation of the \texttt{metrop} function \item \texttt{final}, a vector of length \texttt{d} and mode \verb@"numeric"@, the state after the last iteration \end{itemize} All of the files in the \texttt{tests} directory of the source for the package (not installed but found in the source tarball on CRAN) test the \texttt{metrop} function except those beginning \texttt{temp}, which test the \texttt{temper} function. Since these tests were written many years ago, are spread out over many files, and are not commented, we will not describe them in detail. Suffice it to say that they check every aspect of the functioning of the \texttt{metrop} function. \section{The Temper Function} The R function \texttt{temper} in the \texttt{mcmc} package has an argument \verb@debug = FALSE@ that when \verb@TRUE@ causes extra debugging information to be output. Let \texttt{niter} be the number of iterations \verb@nbatch * blen * nspac@, let \texttt{d} be the dimension of the state vector, and let \texttt{ncomp} be the number of components of the tempering mixture. The result of invoking \texttt{temper} is a list. When \verb@debug = TRUE@ and \verb@parallel = TRUE@ it has the following additional components % which % unif.which % state % log.hastings % unif.hastings % proposal % acceptd % norm % unif.choose % coproposal \begin{itemize} \item \texttt{which}, a vector of length \texttt{niter} and mode \verb@"logical"@ the type of update for each iteration, within component (\texttt{TRUE}) or swap components (\texttt{FALSE}). \item \texttt{unif.which}, a vector of length \texttt{niter} and mode \verb@"numeric"@, the$\text{Uniform}(0, 1)$random variate used to decide which type of update is done. \item \texttt{state}, an \texttt{niter} by \texttt{ncomp} by \texttt{d} array of mode \verb@"numeric"@, the state before iteration \texttt{i} is \verb@state[i, , ]@ \item \texttt{proposal}, an \texttt{niter} by \verb@d + 1@ matrix of mode \verb@"numeric"@, the proposal for iteration \texttt{i} is \verb@proposal[i, ]@ (explanation below) \item \texttt{coproposal}, an \texttt{niter} by \verb@d + 1@ matrix of mode \verb@"numeric"@, the proposal for iteration \texttt{i} is \verb@coproposal[i, ]@ (explanation below) \item \texttt{log.hastings}, a vector of length \texttt{niter} and mode \verb@"numeric"@, the logarithm of the Hastings ratio for each iteration \item \texttt{unif.hastings}, a vector of length \texttt{niter} and mode \verb@"numeric"@, the$\text{Uniform}(0, 1)$random variate compared to the Hastings ratio for each iteration or \texttt{NA} if none is needed (when the log Hastings ratio is nonnegative) \item \texttt{acceptd}, a vector of length \texttt{niter} and mode \verb@"logical"@, the decision for each iteration, accept the proposal (\texttt{TRUE}) or reject it (\texttt{FALSE}) \item \texttt{norm}, an \texttt{niter} by \texttt{d} matrix of mode \verb@"numeric"@, the vector of standard normal random variates used to generate the proposal for iteration \texttt{i} is \verb@z[i, ]@ unless none are needed (for swap updates) when it is \texttt{NA} \item \texttt{unif.choose}, an \texttt{niter} by 2 matrix of mode \verb@"numeric"@, the vector of$\text{Uniform}(0, 1)$random variates used to choose the components to update in iteration \texttt{i} is \verb@unif.choose[i, ]@; in a swap update two are used; in a within-component update only one is used and the second is \texttt{NA} \end{itemize} In a within-component update, one component say \texttt{j} is chosen for update. The \emph{coproposal} is the current value of the state for this component, which is a vector of length \verb@d + 1@, the first component of which is \texttt{j} and the rest of which is \verb@state[i, j, ]@ if we are in iteration \texttt{i}. The \emph{proposal} is a similar vector, the first component of which is again \texttt{j} and the rest of which is a multivariate normal random vector centered at \verb@state[i, j, ]@. The coproposal is the current state; the proposal is the possible value (if accepted) of the state at the next time. In a swap update, two components say \texttt{j1} and \texttt{j2} are chosen for update. Strictly, speaking the coproposal is the pair of vectors \verb@c(j1, state[i, j1, ])@ and \verb@c(j2, state[i, j2, ])@ and the proposal is these swapped, that is, the pair of vectors \verb@c(j2, state[i, j1, ])@ and \verb@c(j1, state[i, j2, ])@ if we are in iteration \texttt{i}. Since, however, there is a lot of redundant information here, the vector \verb@c(j1, state[i, j1, ])@ is output as \verb@coproposal[i, ]@ and the vector \verb@c(j2, state[i, j2, ])@ is output as \verb@proposal[i, ]@. When \verb@debug = TRUE@ and \verb@parallel = FALSE@ the result of invoking \texttt{temper} is a list having the following additional components % which % unif.which % state % log.hastings % unif.hastings % proposal % acceptd % norm % unif.choose \begin{itemize} \item \texttt{which}, a vector of length \texttt{niter} and mode \verb@"logical"@ the type of update for each iteration, within component (\texttt{TRUE}) or jump from one component to another (\texttt{FALSE}). \item \texttt{unif.which}, a vector of length \texttt{niter} and mode \verb@"numeric"@, the$\text{Uniform}(0, 1)$random variate used to decide which type of update is done. \item \texttt{state}, an \texttt{niter} by \verb@d + 1@ matrix of mode \verb@"numeric"@, the state before iteration \texttt{i} is \verb@state[i, ]@ \item \texttt{proposal}, an \texttt{niter} by \verb@d + 1@ matrix of mode \verb@"numeric"@, the proposal for iteration \texttt{i} is \verb@proposal[i, ]@ \item \texttt{log.hastings}, a vector of length \texttt{niter} and mode \verb@"numeric"@, the logarithm of the Hastings ratio for each iteration \item \texttt{unif.hastings}, a vector of length \texttt{niter} and mode \verb@"numeric"@, the$\text{Uniform}(0, 1)$random variate compared to the Hastings ratio for each iteration or \texttt{NA} if none is needed (when the log Hastings ratio is nonnegative) \item \texttt{acceptd}, a vector of length \texttt{niter} and mode \verb@"logical"@, the decision for each iteration, accept the proposal (\texttt{TRUE}) or reject it (\texttt{FALSE}) \item \texttt{norm}, an \texttt{niter} by \texttt{d} matrix of mode \verb@"numeric"@, the vector of standard normal random variates used to generate the proposal for iteration \texttt{i} is \verb@z[i, ]@ unless none are needed (for jump updates) when it is \texttt{NA} \item \texttt{unif.choose}, a vector of length \texttt{niter} and mode \verb@"numeric"@, the$\text{Uniform}(0, 1)\$ random variates used to choose the component to update in iteration \texttt{i} is \verb@unif.choose[i, ]@; in a jump update one is used; in a within-component update none is used and \texttt{NA} is output \end{itemize} All of the files in the \texttt{tests} directory of the source for the package (not installed but found in the source tarball on CRAN) beginning \texttt{temp} test the \texttt{temper} function. They check every aspect of the functioning of the \texttt{temper} function. In the file \texttt{temp-par.R} in the \texttt{tests} directory, the following checks are made according to the comments in that file \begin{enumerate} \item check decision about within-component or jump/swap \item check proposal and coproposal are actually current state or part thereof \item check hastings ratio calculated correctly \item check hastings rejection decided correctly \item check acceptance carried out or not (according to decision) correctly \item check within-component proposal \item check swap proposal \item check standard normal and uniform random numbers are as purported \item check batch means \item check acceptance rates \item check scale vector \item check scale matrix \item check scale list \item check outfun \end{enumerate} In the file \texttt{temp-ser.R} in the \texttt{tests} directory, the all of the same checks are made according to the comments in that file except for check number 2 above, which would make no sense because there is no \texttt{coproposal} component in the serial (\verb@parallel = FALSE@) case. \end{document}