Built-in set of functions to be used in companion with the argument
conv_checker in MCMC. These functions are not intended to be used
in a context other than the MCMC function.
The object LAST_CONV_CHECK is an environment that holds
information regarding the convergence checker used. This information can be
updated every time that the conv_checker function is called by MCMC using
the functions convergence_data_set and convergence_msg_set. The function
convergence_data_get is just a wrapper of get().
The msg member of LAST_CONV_CHECK is resetted before conv_checker is
called.
LAST_CONV_CHECK
convergence_data_set(x)
convergence_data_get(x)
convergence_msg_set(msg = NA_character_)
convergence_msg_get()
convergence_gelman(freq = 1000L, threshold = 1.1, check_invariant = TRUE, ...)
convergence_geweke(
freq = 1000L,
threshold = 0.025,
check_invariant = TRUE,
...
)
convergence_heildel(freq = 1000L, ..., check_invariant = TRUE)
convergence_auto(freq = 1000L)An object of class fmcmc_output_conv_check (inherits from environment) of length 1.
In the case of convergence_data_set, a named list. For
convergence_data_get, a character vector.
Character scalar. Personalized message to print.
Integer scalar. Frequency of checking.
Numeric value. A Gelman statistic below the threshold
will return TRUE.
Logical. When TRUE the function only computes
the Gelman diagnostic using variables with greater than 1e-10 variance.
Further arguments passed to the method.
A function passed to MCMC to check automatic convergence.
convergence_gelman is a wrapper of coda::gelman.diag().
In the case of convergence_geweke, threshold sets the p-value
for the null \(H_0: Z = 0\), i.e. equal means between the first and last
chunks of the chain. See coda::geweke.diag. This implies that the higher
the threshold, the lower the probability of stopping the chain.
In the case that the chain has more than one parameter, the algorithm will return true if and only if the test fails to reject the null for all the parameters.
For the convergence_heildel, see coda::heidel.diag for details.
The convergence_auto function is the default and is just a wrapper
for convergence_gelman and convergence_geweke. This function returns a
convergence checker that will be either of the other two depending on whether
nchains in MCMC is greater than one--in which case it will use the Gelman
test--or not--in which case it will use the Geweke test.
Convergence checkers are simply a function that receives as argument a matrix
(or list of them) with sampled values, and returns a logical scalar with the
value TRUE if the chain converged. An example of a personalized convergence
checker is provided below. The frequency with which the check is performed is
retrieved from the attribute "freq" from the convergence checker function,
i.e., attr(..., "freq"). If missing, convergence will be checked halfway
the number of steps in the chain, i.e., floor(nsteps/2).
# Example 1: Presonalized conv checker --------------------------------------
# Dummy rule, if acceptance rate is near between .2 and .3.
convergence_example <- function(x) {
arate <- 1 - coda::rejectionRate(x)
all(
abs(arate - .25) < .05
)
}
# Tell fmcmc::MCMC what is the frequency
attr(convergence_example, "freq") <- 2e3
set.seed(223)
x <- rnorm(1000)
y <- x * 2 + rnorm(1000)
logpost <- function(p) {
sum(dnorm(y, mean = x * p, log = TRUE))
}
ans <- MCMC(
initial = 0, fun = logpost, nsteps = 5e4,
kernel= kernel_ram(),
conv_checker = convergence_example
)
#> No convergence yet (steps count: 2000). Trying with the next bulk.
#> No convergence yet (steps count: 4000). Trying with the next bulk.
#> No convergence yet (steps count: 6000). Trying with the next bulk.
#> No convergence yet (steps count: 8000). Trying with the next bulk.
#> No convergence yet (steps count: 10000). Trying with the next bulk.
#> Convergence has been reached with 12000 steps. (12000 final count of samples).
# Example 2: Adding information ---------------------------------------------
# Here we do two things: Save a value and set a message for the user
convergence_example_with_info <- structure(function(x) {
arate <- 1 - coda::rejectionRate(x)
# Saving a value
if (!exists("arates", envir = LAST_CONV_CHECK, inherits = FALSE)) {
convergence_data_set(list(arates = arate))
} else {
convergence_data_set(list(
arates = rbind(convergence_data_get("arates"), arate)
))
}
# Setting up the message
convergence_msg_set(
sprintf("Current Avg. Accept. Rate: %.2f", mean(arate))
)
all(
abs(arate - .25) < .05
)
}, freq = 2000)
ans <- MCMC(
initial = 0, fun = logpost, nsteps = 5e4,
kernel= kernel_ram(),
conv_checker = convergence_example_with_info,
seed = 112,
progress = FALSE
)
#> No convergence yet (steps count: 2000). Current Avg. Accept. Rate: 0.56 Trying with the next bulk.
#> No convergence yet (steps count: 4000). Current Avg. Accept. Rate: 0.40 Trying with the next bulk.
#> No convergence yet (steps count: 6000). Current Avg. Accept. Rate: 0.34 Trying with the next bulk.
#> No convergence yet (steps count: 8000). Current Avg. Accept. Rate: 0.31 Trying with the next bulk.
#> Convergence has been reached with 10000 steps. Current Avg. Accept. Rate: 0.30 (10000 final count of samples).