R/AdMitMH.R

In AdMit: Adaptive Mixture of Student-t Distributions

## Function which performs (independence) Metropolis-Hastings (MH) sampling
## using an adaptive mixture of Student-t densities as the candidate density
## __input__
## N        : [integer>0] length of the MCMC output (default: 1e5)
## KERNEL   : [function] which computes the kernel !! must be vectorized !!
## mit      : [list] with mixture information (default: univariate Cauchy)
## ...      : additional parameters used by 'KERNEL'
## __output__
## [list] with the following components:
## $theta   : [Nxk matrix] sample generated by the (independence) MH algorithm
## $accept  : [double] acceptance rate in the (independence) MH algorithm
## __20080429__
'AdMitMH' <- function(N=1e5, KERNEL, mit=list(), ...)
  {
    if (N<2)
      stop ("'N' should be at least larger than 2")
    if (missing(KERNEL))
      stop ("'KERNEL' is missing in 'AdMitMH'")
    KERNEL <- match.fun(KERNEL)
    if (!any(names(formals(KERNEL))=="log"))
      stop ("'KERNEL' MUST have the logical argument 'log' which returns the (natural) logarithm of 'KERNEL' if 'log=TRUE'")      
    theta <- rMit(N, mit)
    lnw <- fn.w(theta, KERNEL=KERNEL, mit=mit, log=TRUE, ...)
    k <- ncol(mit$mu)
    
    r <- .C('fnMH_C',
            theta = as.double(as.vector(t(theta))),
            N = as.integer(N),
            k = as.integer(k),
            lnw = as.double(lnw),
            u = as.double(stats::runif(N)),
            draws = vector('double',N*k),
            ns = as.integer(0),
            PACKAGE = 'AdMit',
            NAOK = TRUE)

    draws <- matrix(r$draws, N, k, byrow=TRUE)
    dimnames(draws) <- list(1:N, paste("k", 1:k, sep=""))

    list(draws=draws,
         accept=as.numeric(r$ns/N))
  }