## 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))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.