FUNCTIONS FOR CONTROLLING AN MCMC RUN The functions mcmc and more_mcmc are the top-level functions for doing an MCMC run. They return the states from the simulation, along with other information, such as acceptance probabilities. The distribution to sample from is specified by a function for the log probability density provided by the user (or that is one supplied with this package). The MCMC updates to do each iteration are specified by updates functions (either supplied with this package or written by the user), which may be general-purpose or specialized to the particular distribution being sampled. FUNCTION TO DO AN MCMC SIMULATION mcmc (lpr, initial, iterations, ..., updates=list(), rep=1, fix=NULL, initial.p=NULL, last="lpr", rec="acc", ave="apr", lpr.initial=NULL, initial.arg=initial, seed=1, saved.rand=NULL) Does the specified number of iterations of an MCMC simulation, with each iteration consisting of a specified sequence of updates. Most of the arguments to mcmc are also returned as part of its result, so that these arguments will be available without re-specification if more iterations are done with the more_mcmc function (see below), and so that a saved result will contain documentation on what was done. Arguments are as follows: lpr Function for evaluating the log probability density of a state (up to an arbitary additive constant), and possibly its gradient. May have attributes specifiying bounds on variables. initial Initial state of the chain. May be a single vector, or a list of vectors. The total length of these vectors is the dimension of the state (excluding the momentum part if that exists). The vectors may have named elements, which will be used as names for the result. iterations The number of iterations to simulate. ... Zero or more arguments describing the updates to be done in an iteration. Each argument is either an update function, or a list whose first element is an update function and the remaining elements arguments to it, or a vector of strings, which restricts the following updates to just the parts of the state named, or a vector of integer indexes, which further restricts the part of the state updated. (See below for details.) updates A list of additional updates, done after those specified by the ... arguments. Specifying updates this way may sometimes be more convenient. (Optional) rep Number of times to repeat the list of updates for each iteration (default 1). fix A list of vectors specifying parts of the state that are fixed during the simulation. These parts will be in the list passed to the lpr function, but are not part of the Markov chain state, and not present in initial. (Optional) initial.p Initial values for the momentum variables, either a vector or a list with elements that are a subset of those in initial. Specification of initial.p also implies that these momentum variables are part of the state. (Optional) last A vector of names of quantities returned by an update function, whose values should be recorded for the last update in the last repetition of the sequence of update functions. (Default is "lpr" only.) rec A vector of names of quantities returned by update functions, whose values should be recorded after all updates in the last repetition of the sequence of update functions. (Default is "acc" only.) ave A vector of names of quantities returned by update functions, whose values after each update should be averaged over the repetitions of the sequence of update functions. (Default is "apr" only.) lpr.initial The value of lpr for the initial state. Must, of course, be correct, or wrong answers will result! (Optional) initial.arg The initial argument for the sampling run, typically the same as initial, but provided for use by more_mcmc (see below). (Default is initial.) seed The random number seed for the run. This should be set by the user to different values if more than one MCMC run is done for the same problem. Ignored if saved.rand is set. (Default is 1.) saved.rand If not NULL, a saved state of the random number generator (wrapped in an environment, to suppress a long display, and allow error checking). This state is restored instead of the state being set from the seed. (Optional) The result of mcmc is a list with the following elements: final Final state of variables (except momentum) after all iterations, with names as in initial. Does not include variables that are fixed, only those in "initial". final.p Final state of momentum variables after all iterations (absent if momentum is not saved as part of the state). Includes only variables that were specified by initial.p. q Matrix or list of matrices holding the state (except momentum) after each iteration, with one row for each iteration, and column names taken from the "initial" argument. p Matrix or list of matrices holding the momentum part of the state after each iteration, with column names taken from the corresponding parts of "initial". (Absent if momentum is not saved as part of the state.) last A matrix with one row for each iteration and column names from the argument last, except columns for quantities that aren't returned by the last update function called are omitted. No last element is present if the last argument was NULL. rec A matrix with one row for each iteration and column names from the rec argument, repeated for each update function, omitting values that aren't returned. These values are from the last repetition of the updates. No rec element is present if the rec argument is NULL. ave A matrix with one row for each iteration and column names from the ave argument, repeated for each update function, omitting values that aren't returned. These values are averages over all repetitions. No ave element is present if the ave argument is NULL. lpr.final The value of lpr for the final state. saved.rand State of the random number generator after the last iteration. lpr.arg The lpr argument passed. initial.arg The initial.arg argument passed, usually the same as initial. rep.arg The rep argument passed (may be the default of 1). seed.arg The seed argument passed (may be the default of 1). fix.arg The fix argument passed, or absent. last.arg The last argument passed, or absent. rec.arg The rec argument passed, or absent ave.arg The ave argument passed, or absent updates The MCMC updates done (ie, a list of what was passed as ... concatenated with the updates argument). As described above, the updates specified via the ... arguments or via the updates argument may include vectors of strings or integers that restrict subsequent updates to a portion of the state. When a string vector is encountered, it is concatenated with any immediately preceding string vectors and is then used to restrict updates to the subset of the state consisting only of list elements with names in this vector of strings. An integer vector is concatenated with any immediately preceding integer vectors and is then used to restrict updates to a subset of the set of values that would otherwise have been updated according to the usual R indexing convention (eg, 1:10 specifies only the first ten values, and -10 specifies all except the tenth value). Integer restrictions are allowed only when the state is a vector, or the state is a list with one element, or the current string restriction specifies just one element of the state. A entry of TRUE in the list of updates will cancel all restrictions. A string restriction cancels any previous string restriction (that doesn't immediately precede it) and any integer restriction. An integer restriction cancels any previous restriction specified by an integer vector (that doesn't immediately precede it), but it does not cancel the current string restriction. Updates from ... or in the updates argument are specified by either an update function or a list with the first element being an update function and the remaining elements being arguments for this update function. For special purpose update functions (with a "special" attribute that is TRUE), these arguments are passed on unchanged. Unnamed arguments of general-purpose update functions are also passed on unchanged. A named argument of a general-purpose update function may be o A numeric scalar. This will be passed unchanged to the update function, which might interpret it as applying to every value in the state, or (like "rep") it might have a meaning that is not related to individual values in the state. This scalar may have attributes, which may be of any type, and whose meaning (if any) will depend on the update function. o A list, if the state is a list, in which case it must have elements for all parts of the state being updated, that are numeric scalars (applying to that whole part) or numeric vectors of the same length as that part of the state. This list will be converted to a vector according to the current string and integer restrictions before it is passed to the update function. o A numeric vector, if the state is a vector, in which case it must either be scalar or of the same length as the state. If not a scalar, it will be subsetted with the current integer restriction (if any) before being passed. o A function, in which case this function will be called each time the update is done to provide the value to be passed to the update function. o Something else (eg, a logical vector or a string), which will be passed unchanged to the update function. If the state is a vector, a named argument that is a function will be called before each update with a single argument that is a copy of the state in which positions that are being updated (according to to the integer restriction) are set to NA, except if there is no integer restriction, no argument will be passed. It should return either a scalar or a vector of values for the positions being updated. If the state is a list, a named argument that is a function will be called with a single argument that is a list containing all elements of the state not being updated (according to the string restriction), plus fixed elements. If there is an integer restriction as well, the element being updated will appear as an element in this argument, with NA in positions being updated. It should return a scalar, or a vector of the same length as the total number of values being updated, or a list with elements corresponding to the parts being updated, in which each element is either a scalar or a vector of length equal to the size of that part, or the size of the subset being updated (when there is an integer restriction). If the state is a list and there is no current restriction, the argument passed will be a list of only the fixed elements (which might be an empty list). When mcmc creates a wrapper for the lpr function it is passed, the wrapper function will have a grad argument if the original lpr function did. If the original lpr function had a ch.value argument, any wrapper passed to a specialized update function will have the same set of ch.value, ch.elem, ch.pos, and lpr.value arguments as the original lpr function. A wrapper for lpr passed to a general-purpose update function will not have ch.value or related arguments if the state is a list and is not currently restricted to just one element, even if the original lpr function has these arguments. If the state is a vector, the wrapper will have ch.value, ch.pos, and lpr.value arguments if the original lpr function did. If the state is a list and only one element is currently being updated (possibly with an integer restriction), the wrapped lpr function will have ch.value, ch.pos, and lpr.value arguments if the original lpr function has ch.value, ch.elem, ch.pos, and lpr.value arguments. It is permissible for the initial.p argument to have only a subset of the elements in initial, in which case momentum variables are retained in the state for only that subset. However, at a point where an update function with an initial.p argument is used, it must be that either all of the elements in the current string restriction are present in initial.p, or none of them are present (in which case no initial.p argument will be passed to the update function). DO MORE MCMC ITERATIONS OF AN MCMC SIMULATION more_mcmc (previous, iterations) Extends an MCMC run by doing more iterations. Arguments are as follows: previous The result of a previous mcmc or more_mcmc call iterations The number of additional iterations to do The result returned is of the same form as is returned by the mcmc function. Doing a run of n iterations with mcmc, and then extending by m more iterations with more_mcmc should produced the same result as doing n+m iterations with mcmc would have. Note that the state of the random number generator at the end of the first run is saved, and restored for the additional iterations, and that the initial state and the seed that were passed to mcmc are passed on to more_mcmc for inclusion in its result (as initial.arg and seed.arg elements).