expect_mcmc_reversible {mcunit}

Test of MCMC chain assuming reversibility of the chain

Description

Test of MCMC steps having the correct stationary distribution assuming reversibility of the chain. Uses ideas from Besag and Clifford (1989) as extended to possible ties in Gandy and Scott (2020).

Usage

expect_mcmc_reversible(
  object,
  control = NULL,
  thinning = NULL,
  nsteps = 10,
  p = 1,
  tolerance = 1e-08
)

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.

nsteps

Number of steps of the chain to use. Default: 10.

p

The probability with which the MCMC updates the parameter as opposed to the data in a given step. If less than 1.0, then the MCMC is a random scan on both data and parameters. Default: 1.0.

tolerance

Absolute error where value of the samplers are treated as equal. Default: 1e-8.

Value

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

References

Besag J, Clifford P (1989). “Generalized Monte Carlo significance tests.” Biometrika, 76(4), 633–642. doi: 10.1093/biomet/76.4.633, https://doi.org/10.1093/biomet/76.4.633.

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

See Also

expect_mcmc

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
               },
               test=function(x) x
               )
expect_mcmc_reversible(object)