sample_meeting_times: Sample meeting times associated with coupled lagged Gibbs...

In pierrejacob/montecarlodsm: Gibbs sampler for Dempster-Shafer inference in Categorical distributions

Description

Sample meeting times to monitor convergence of the Gibbs sampler, as in the article Estimating Convergence of Markov chains with L-Lag Couplings by Niloy Biswas, Pierre E. Jacob, Paul Vanetti, available at https://arxiv.org/abs/1905.09971.

The coupled chains are generated using a mixture of coupled kernels; with probability omega, a "common random numbers (CRN)" coupling is performed; otherwise a (nearly) maximal coupling is used in a Gibbs sweep.

Usage

sample_meeting_times(
  counts,
  lag = 1,
  omega = 0.9,
  max_iterations = 1e+05,
  removezero = TRUE
)

Arguments

counts	vector of counts; could include zeros.
lag	lag between the chains; defaults to 1.
omega	probability of a coupled "CRN" step as opposed to maximal coupling step; defaults to 0.9.
max_iterations	number of iterations after which to stop, in case meeting hasn't occurred; defaults to 1e5.
removezero	remove zeros from counts before running coupled Gibbs sampler; defaults to TRUE

Value

An integer representing the meeting time; or +Inf if meeting has not occurred before 'max_iterations'. The meeting times can be turned into upper bounds on the total variation (TV) distance between the Markov chain at some iteration and the limiting distribution, using the function tv_upper_bound.

Examples

## Not run: 
nrep <- 100
meeting_times <- sapply(1:nrep, function(irep) sample_meeting_times(counts = c(3,2,0,1)))
tmax <- floor(max(meeting_times)*1.2)
ubounds <- sapply(1:tmax, function(t) tv_upper_bound(meeting_times, 1, t))
plot(x = 1:tmax, y = ubounds, type = 'l', xlab = "iteration", ylab = "TV upper bounds")

## End(Not run)

pierrejacob/montecarlodsm documentation built on June 16, 2021, 1:06 p.m.