expect_mcmc {mcunit}

Test of MCMC chain

Description

Test of MCMC steps having the correct stationary distribution without assuming reversibility of the chain. Details of this are in Gandy and Scott (2020); it uses ideas described in the appendix of Gandy and Veraart (2017).

Usage

expect_mcmc(object, control = NULL, thinning = NULL, joint = FALSE)

Arguments

object

A list describing the MCMC sampler with the following elements:

• genprior: A function with no arguments that generates a sample of the prior distribution. No default value.

• gendata: A function that takes as input the parameter value (of the type generated by genprior) and returns the observed data as an arbitrary R object. No default value.

• stepMCMC: A function that takes three arguments:

  • theta: the current position of the chain (of the same type as produced by the prior),

  • dat: the observed data (of the same type as produced by gendat)

  • thinning: the number of steps the chain should take. 1 corresponds to one step.

• test: Function that takes either one or two arguments, and returns a vector with components which will be used for checking the MCMC sampler. The first argument is interpreted as a parameter value, and if a second argument exists, it is interpreted as a data value. An example is the identity function: function(x) x. Alternatively, if you have access to the model's likelihood function, you could use: function(x,y) likelihood(x,y).

control

a list controlling the algorithm

• n number of samples to be taken in the first step. Default: 1e3

• maxseqsteps: Number of sequential attempts to use. Default: 7.

• incn: Factor by which to multiply n from the second sequential attempt onward. Default: 4.

• level: bound on the type I error, ie the probability of wrongly rejecting a sampler with the correct distribution. Default: 1e-5.

• debug: If positive then debug information will be printed via 'message()'. Default: 0.

thinning

Amount of thinning for the MCMC chain. 1 corresponds to no thinning. Default: automatically computed to ensure an autocorrelation of at most 0.5 at lag 1.

joint

If TRUE, then the MCMC uses systematic scan of both data and parameters, rather than just updating parameters with the sampler to be tested. Default: FALSE.

Value

The first argument, invisibly, to allow chaining of expectations.

References

Gandy A, Scott J (2020). “Unit Testing for MCMC and other Monte Carlo Methods.” https://arxiv.org/abs/2001.06465.

Gandy A, Veraart LAM (2017). “A Bayesian Methodology for Systemic Risk Assessment in Financial Networks.” Management Science, 63(12), 4428–4446. doi: 10.1287/mnsc.2016.2546, https://doi.org/10.1287/mnsc.2016.2546.

See Also

expect_mcmc_reversible

Examples

object <- list(genprior=function() rnorm(1),
               gendata=function(theta) rnorm(5,theta),
               stepMCMC=function(theta,data,thinning){
                  f <- function(x) prod(dnorm(data,x))*dnorm(x)  
                  for (i in 1:thinning){
                     thetanew = rnorm(1,mean=theta,sd=1)
                     if (runif(1)<f(thetanew)/f(theta))
                     theta <- thetanew
                  }
                  theta
               }
               )
expect_mcmc(object)