stable.GR {stableGR}

Gelman-Rubin diagnostic using stable variance estimators

Description

This function uses fast and strongly consistent estimators estimators of Monte Carlo variance to calculate the Gelman-Rubin convergence diagnostic for Markov chain Monte Carlo. A univariate ‘potential scale reduction factor’ (PSRF) is calculated for each variable in x. For multivariate chains, a multivariate PSRF is calculated to take into account the interdependence of the chain's components. The PSRFs decrease to 1 as the chain length increases. When the PSRF becomes sufficiently close to 1, the sample collected by the Markov chain has converged to the target distribution.

Usage

stable.GR(
  x,
  multivariate = TRUE,
  mapping = "determinant",
  method = "lug",
  size = NULL,
  autoburnin = FALSE,
  blather = FALSE
)

Arguments

x

a list of matrices, where each matrix represents one Markov chain sample. Each row of the matrices represents one step of the chain. Each column of the matrices represents one variable. A list with a single matrix (chain) is allowed. Optionally, this can be an mcmclist object.

multivariate

a logical flag indicating whether the multivariate potential scale reduction factor should be calculated for multivariate chains.

mapping

the function used to map the covariance matrix to a scalar. This is one of “determinant” (determinant of the covariance matrix, the default) or “maxeigen” (the largest eigenvalue of the covariance matrix).

method

the method used to compute the standard error of the chains. This is one of “lug” (lugsail, the default), “bm” (batch means), “obm” (overlapping batch means), “tukey” (spectral variance method with a Tukey-Hanning window), or “bartlett” (spectral variance method with a Bartlett window).

size

options are NULL (default, which calculates an ideal batch size), character values of sqroot and cuberoot, or any numeric value between 1 and n. Size represents the batch size in “bm” (batch means) and the truncation point in “bartlett” and “tukey”. sqroot means size is floor(n^(1/2) and cuberoot means size is floor(n^(1/3)).

autoburnin

a logical flag indicating whether only the second half of the series should be used in the computation. If set to TRUE and start(x) is less than end(x)/2 then start of series will be adjusted so that only second half of series is used.

blather

a logical flag indicating whether to include additional output.

Value

psrf

A vector containing the point estimates of the PSRF.

mpsrf

A scalar point estimate of the multivariate PSRF.

means

A vector containing the sample means based on the chains provided.

n.eff

A scalar point estimate of the effective sample size.

blather

Either FALSE or a list containing intermediate calculations.

Theory

Gelman and Rubin (1992) and Brooks and Gelman (1998) first constructed the univariate and multivariate potential scale reduction factors (PSRF), respectively, to diagnose Markov chain convergence. The function stable.GR stabilizes the PSRF and improves the PSRF's efficiency by incorporating lugsail estimators for the target variance. The PSRF decreases to 1 as the chain length increases; when the PSRF becomes sufficiently close to 1, the sample collected by the Markov chain has converged to to the target distribution. A PSRF convergence threshold can be calculated using choosepsrf.

References

Vats, D. and Knudson, C. Revisiting the Gelman-Rubin Diagnostic. arXiv:1812.09384.

Vats, D. and Flegal, J. Lugsail lag windows and their application to MCMC. arXiv: 1809.04541.

Flegal, J. M. and Jones, G. L. (2010) Batch means and spectral variance estimators in Markov chain Monte Carlo. The Annals of Statistics, 38, 1034–1070.

Gelman, A and Rubin, DB (1992) Inference from iterative simulation using multiple sequences, Statistical Science, 7, 457-511.

Brooks, SP. and Gelman, A. (1998) General methods for monitoring convergence of iterative simulations. Journal of Computational and Graphical Statistics, 7, 434-455.

Examples

library(stableGR)
set.seed(100)
p <- 2
n <- 100 # n is tiny here purely for demo purposes.
# use n much larger for real problems!

sig.mat = matrix(c(1, .8, .8, 1), ncol = 2, nrow = 2)

# Making 3 chains
chain1 <- mvn.gibbs(N = n, p = p, mu = rep(1,p), sigma = sig.mat)
chain2 <- mvn.gibbs(N = n, p = p, mu = rep(1,p), sigma = sig.mat)
chain3 <- mvn.gibbs(N = n, p = p, mu = rep(1,p), sigma = sig.mat)

# find GR diagnostic using all three chains
x <- list(chain1, chain2, chain3)
stable.GR(x)