AdMitMH: Independence Chain Metropolis-Hastings Algorithm using an...
In AdMit: Adaptive Mixture of Student-t Distributions

Description Usage Arguments Details Value Note Author(s) References See Also Examples

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

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

`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`.

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 Hxd) containing the vectors of modes (in row) of the mixture components.
Sigma: matrix (of size Hxd*d) 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 (>=1) is the number of components and d (>=1) is the dimension of the first argument in KERNEL. Typically, mit is estimated by the function AdMit.

A list with the following components:

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

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

Further details and examples of the R package AdMit can be found in Ardia, Hoogerheide and van Dijk (2009a,b). See also the package vignette by typing vignette("AdMit").

Further information on the Metropolis-Hastings algorithm can be found in Chib and Greenberg (1995) and Koop (2003).

Please cite the package in publications. Use citation("AdMit").

David Ardia

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. doi: 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. doi: 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.

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.

  ## 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)