Terminates an MCMC run when the Gelman-Rubin statistic drops below the
specified value.
The Gelman-Rubin statistic, also referred to as the potential scale reduction
factor (PSRF), compares the variance of the sample pooled from multiple runs to
the sum of variances calculated from individual runs. Accordingly, it can only
be calculated when two or more independent runs are performed, and its value
tends to unity (1) as the runs converge.
The number of samples to be removed as burnin before calculating the test
statistic is determined using the `burninMethod`. If the `"ESS"` (default) or
`"SEM"` options are chosen, different burnin lengths are tested, increasing
from 0 to 50% (for `"ESS"`) or 100% (for `"SEM"`) of the length of the trace in
increments of 10 samples. The `"ESS"` option calculates effective sample sizes
(ESS) for all monitored parameters after removing the number of samples
corresponding to each candidate burnin length. The best burnin length for a
given parameter is the one that maximizes its ESS value. The `"SEM"` option
instead calculates the standard error of the mean (SEM), and the best burnin
length for a given parameter is the one that minimizes its SEM value. In both
cases, the final burnin length is set to the maximum of the parameter-specific
burnin lengths.
Alternatively, the user may set `burninMethod` to `"fixed"`, which discards a
constant fraction of the samples collected up to that point. This fraction can
be specified using the `burnin` argument, and is set to 0.25 by default. The
`burnin` argument has no effect if the `"ESS"` or `"SEM"` options are chosen,
and a corresponding warning is displayed if the user explicitly sets the
argument without specifying `burninMethod="fixed"`. The `"fixed"` option is
appropriate for analyses with very long parameter traces and large numbers of
monitored variables, for which the automatic burnin determination may be too
computationally demanding.
See also the tutorial on [convergence assessment](https://revbayes.github.io/tutorials/convergence/).
# Binomial example: estimate success probability given 7 successes out of 20 trials
r ~ dnExp(10)
p := Probability(ifelse(r < 1, r, 1))
n <- 20
k ~ dnBinomial(n, p)
k.clamp(7)
mymodel = model(k)
moves = VectorMoves()
moves.append( mvSlide(r, delta=0.1, weight=1) )
paramFile = "parameters.log"
monitors = VectorMonitors()
monitors.append( mnModel(filename=paramFile, printgen=100, p) )
# Stop when the potential scale reduction factor falls below 1.01
stopping_rules[1] = srGelmanRubin(1.01, file = paramFile, freq = 1000)
# Create the MCMC object.
# Set nruns = 2 to ensure the Gelman-Rubin statistic is applicable
mymcmc = mcmc(mymodel, monitors, moves, nruns = 2)
# Begin the MCMC run
mymcmc.run(rules = stopping_rules)