Transition Kernels for MCMC

The function kernel_new is a helper function that allows creating fmcmc_kernel objects which are passed to the MCMC() function.

kernel_new(proposal, ..., logratio = NULL, kernel_env = new.env(hash = TRUE))

Arguments

proposal, logratio: Functions. Both receive a single argument, an environment. This functions are called later within MCMC (see details).
...: In the case of kernel_new, further arguments to be stored with the kernel.
kernel_env: Environment. This will be used as the main container of the kernel's components. It is returned as an object of class c("environment", "fmcmc_kernel").

Value

An environment of class fmcmc_kernel which contains the following:

proposal A function that receives a single argument, an environment. This is the proposal function used within MCMC().
logratio A function to compute log ratios of the current vs the proposed step of the chain. Also used within MCMC().
... Further arguments passed to kernel_new.

Details

The objects fmcmc_kernels are environments that in general contain the following objects:

proposal: The function used to propose changes in the chain based on the current state. The function must return a vector of length equal to the number of parameters in the model.
logratio: This function is called after a new state has been proposed, and is used to compute the log of the Hastings ratio.

In the case that the logratio function is not specified, then it is assumed that the transition kernel is symmetric, this is, log-ratio is then implemented as function(env) {env$f1 - env$f0}
...: Further objects that are used within those functions.

Both functions, proposal and logratio, receive a single argument, an environment, which is passed by the MCMC() function during each step using the function environment().

The passed environment is actually the environment in which the MCMC function is running, in particular, this environment contains the following objects:

Object		Description
`i`		Integer. The current iteration.
`theta1`		Numeric vector. The last proposed state.
`theta0`		Numeric vector. The current state
`f`		The log-unnormalized posterior function (a wrapper of `fun` passed to MCMC).
`f1`		The last value of `f(theta1)`
`f0`		The last value of `f(theta0)`
`kernel`		The actual `fmcmc_kernel` object.
`ans`		The matrix of samples defined up to `i - 1`.

These are the core component of the MCMC function. The following block of code is how this is actually implemented in the package:

for (i in 1L:nsteps) {
  # Step 1. Propose
  theta1[] <- kernel$proposal(environment())
  f1       <- f(theta1)
  
  # Checking f(theta1) (it must be a number, can be Inf)
  if (is.nan(f1) | is.na(f1) | is.null(f1)) 
    stop(
      "fun(par) is undefined (", f1, ")",
      "Check either -fun- or the -lb- and -ub- parameters.",
      call. = FALSE
    )
  
  # Step 2. Hastings ratio
  if (R[i] < kernel$logratio(environment())) {
    theta0 <- theta1
    f0     <- f1
  }
  
  # Step 3. Saving the state
  ans[i,] <- theta0
  
}

For an extensive example on how to create new kernel objects see the vignette vignette("user-defined-kernels", "fmcmc").

Behavior

In some cases, calls to the proposal() and logratio() functions in fmcmc_kernels can trigger changes or updates of variables stored within them. A concrete example is with adaptive kernels.

Adaptive Metropolis and Robust Adaptive Metropolis implemented in the functions kernel_adapt() and kernel_ram() both update a covariance matrix used during the proposal stage, and furthermore, have a warmup stage that sets the point at which both will start adapting. Because of this, both kernels have internal counters of the absolute step count which allows them activating, scaling, etc. the proposals correctly.

When running multiple chains, MCMC will create independent copies of a baseline passed fmcmc_kernel object. These are managed together in a fmcmc_kernel_list object.
Even if the chains are run in parallel, if a predefined kernel object is passed it will be updated to reflect the last state of the kernels before the MCMC call returns.

References

Brooks, S., Gelman, A., Jones, G. L., & Meng, X. L. (2011). Handbook of Markov Chain Monte Carlo. Handbook of Markov Chain Monte Carlo.

Examples


# Example creating a multivariate normal kernel using the mvtnorm R package
# for a bivariate normal distribution
library(mvtnorm)
#> Error in library(mvtnorm): there is no package called ‘mvtnorm’

# Define your Sigma
sigma <- matrix(c(1, .2, .2, 1), ncol = 2)

# How does it looks like?
sigma
#>      [,1] [,2]
#> [1,]  1.0  0.2
#> [2,]  0.2  1.0
#      [,1] [,2]
# [1,]  1.0  0.2
# [2,]  0.2  1.0

# Create the kernel
kernel_mvn <- kernel_new(
  proposal = function(env) {
  env$theta0 + as.vector(mvtnorm::rmvnorm(1, mean = 0, sigma = sigma.))
  },
  sigma. = sigma
)

# As you can see, in the previous call we passed sigma as it will be used by
# the proposal function
# The logaratio function was not necesary to be passed since this kernel is
# symmetric.