AdMit (version 2.1.9)

AdMitMH: Independence Chain Metropolis-Hastings Algorithm using an Adaptive Mixture of Student-t Distributions as the Candidate Density

Description

Performs independence chain Metropolis-Hastings (M-H) sampling using an adaptive mixture of Student-t distributions as the candidate density

Usage

AdMitMH(N = 1e5, KERNEL, mit = list(), ...)

Arguments

N

number of draws generated by the independence chain M-H algorithm (positive integer number). Default: N = 1e5.

KERNEL

kernel function of the target density on which the adaptive mixture is fitted. This function should be vectorized for speed purposes (i.e., its first argument should be a matrix and its output a vector). Moreover, the function must contain the logical argument log. If log = TRUE, the function returns (natural) logarithm values of kernel function. NA and NaN values are not allowed. (See the function AdMit for examples of KERNEL implementation.)

mit

list containing information on the mixture approximation (see *Details*).

…

further arguments to be passed to KERNEL.

Value

A list with the following components:

draws: matrix (of size N\(\times d\)) of draws generated by the independence chain M-H algorithm, where N \((\geq 1)\) is the number of draws and \(d (\geq 1)\) is the dimension of the first argument in KERNEL.

accept: acceptance rate of the independence chain M-H algorithm.

Details

The argument mit is a list containing information on the adaptive mixture of Student-t distributions. The following components must be provided:

p

vector (of length \(H\)) of mixing probabilities.

mu

matrix (of size \(H \times d\)) containing the vectors of modes (in row) of the mixture components.

Sigma

matrix (of size \(H \times d^2\)) containing the scale matrices (in row) of the mixture components.

df

degrees of freedom parameter of the Student-t components (real number not smaller than one).

where \(H (\geq 1)\) is the number of components and \(d (\geq 1)\) is the dimension of the first argument in KERNEL. Typically, mit is estimated by the function AdMit.

References

Ardia, D., Hoogerheide, L.F., van Dijk, H.K. (2009a). AdMit: Adaptive Mixture of Student-t Distributions. The R Journal 1(1), pp.25-30. 10.32614/RJ-2009-003

Ardia, D., Hoogerheide, L.F., van Dijk, H.K. (2009b). Adaptive Mixture of Student-t Distributions as a Flexible Candidate Distribution for Efficient Simulation: The R Package AdMit. Journal of Statistical Software 29(3), pp.1-32. 10.18637/jss.v029.i03

Chib, S., Greenberg, E. (1995). Understanding the Metropolis-Hasting Algorithm. The American Statistician 49(4), pp.327-335.

Koop, G. (2003). Bayesian Econometrics. Wiley-Interscience (London, UK). ISBN: 0470845678.

See Also

AdMitIS for importance sampling using an adaptive mixture of Student-t distributions as the importance density, AdMit for fitting an adaptive mixture of Student-t distributions to a target density through its KERNEL function; the package coda for MCMC output analysis.

Examples

# NOT RUN {
<!-- % -->
# }
# NOT RUN {
  ## NB : Low number of draws for speedup. Consider using more draws!
  ## Gelman and Meng (1991) kernel function
  GelmanMeng <- function(x, A = 1, B = 0, C1 = 3, C2 = 3, log = TRUE)
  {
    if (is.vector(x))
      x <- matrix(x, nrow = 1)
      r <- -.5 * (A * x[,1]^2 * x[,2]^2 + x[,1]^2 + x[,2]^2
                - 2 * B * x[,1] * x[,2] - 2 * C1 * x[,1] - 2 * C2 * x[,2])
      if (!log)
        r <- exp(r)
      as.vector(r)
    }

  ## Run the AdMit function to fit the mixture approximation
  set.seed(1234)
  outAdMit <- AdMit(KERNEL = GelmanMeng, 
                    mu0 = c(0.0, 0.1), control = list(Ns = 1e4))

  ## Run M-H using the mixture approximation as the candidate density
  outAdMitMH <- AdMitMH(N = 1e4, KERNEL = GelmanMeng, mit = outAdMit$mit)
  options(digits = 4, max.print = 40)
  print(outAdMitMH)

  ## Use functions provided by the package coda to obtain
  ## quantities of interest for the density whose kernel is 'GelmanMeng'
  library("coda")
  draws <- as.mcmc(outAdMitMH$draws)
  draws <- window(draws, start = 1001)
  colnames(draws) <- c("X1", "X2")
  summary(draws)
  summary(draws)$stat[,3]^2 / summary(draws)$stat[,4]^2 ## RNE
  plot(draws)
# }