MCMCmetrop1R: Metropolis Sampling from User-Written R function

In MCMCpack: Markov Chain Monte Carlo (MCMC) Package

Metropolis Sampling from User-Written R function

Description

This function allows a user to construct a sample from a user-defined continuous distribution using a random walk Metropolis algorithm.

Usage

MCMCmetrop1R(
  fun,
  theta.init,
  burnin = 500,
  mcmc = 20000,
  thin = 1,
  tune = 1,
  verbose = 0,
  seed = NA,
  logfun = TRUE,
  force.samp = FALSE,
  V = NULL,
  optim.method = "BFGS",
  optim.lower = -Inf,
  optim.upper = Inf,
  optim.control = list(fnscale = -1, trace = 0, REPORT = 10, maxit = 500),
  ...
)

Arguments

fun	The unnormalized (log)density of the distribution from which to take a sample. This must be a function defined in R whose first argument is a continuous (possibly vector) variable. This first argument is the point in the state space at which the (log)density is to be evaluated. Additional arguments can be passed to fun() by inserting them in the call to MCMCmetrop1R(). See the Details section and the examples below for more information.
theta.init	Starting values for the sampling. Must be of the appropriate dimension. It must also be the case that fun(theta.init, ...) is greater than -Inf if fun() is a logdensity or greater than 0 if fun() is a density.
burnin	The number of burn-in iterations for the sampler.
mcmc	The number of MCMC iterations after burnin.
thin	The thinning interval used in the simulation. The number of MCMC iterations must be divisible by this value.
tune	The tuning parameter for the Metropolis sampling. Can be either a positive scalar or a k-vector, where k is the length of \theta.
verbose	A switch which determines whether or not the progress of the sampler is printed to the screen. If verbose is greater than 0 the iteration number, the \theta vector, the function value, and the Metropolis acceptance rate are sent to the screen every verboseth iteration.
seed	The seed for the random number generator. If NA, the Mersenne Twister generator is used with default seed 12345; if an integer is passed it is used to seed the Mersenne twister. The user can also pass a list of length two to use the L'Ecuyer random number generator, which is suitable for parallel computation. The first element of the list is the L'Ecuyer seed, which is a vector of length six or NA (if NA a default seed of rep(12345,6) is used). The second element of list is a positive substream number. See the MCMCpack specification for more details.
logfun	Logical indicating whether fun returns the natural log of a density function (TRUE) or a density (FALSE).
force.samp	Logical indicating whether the sampling should proceed if the Hessian calculated from the initial call to optim routine to maximize the (log)density is not negative definite. If force.samp==TRUE and the Hessian from optim is non-negative definite, the Hessian is rescaled by subtracting small values from it's main diagonal until it is negative definite. Sampling proceeds using this rescaled Hessian in place of the original Hessian from optim. By default, if force.samp==FALSE and the Hessian from optim is non-negative definite, an error message is printed and the call to MCMCmetrop1R is terminated. Please note that a non-negative Hessian at the mode is often an indication that the function of interest is not a proper density. Thus, force.samp should only be set equal to TRUE with great caution.
V	The variance-covariance matrix for the Gaussian proposal distribution. Must be a square matrix or NULL. If a square matrix, V must have dimension equal to the length of theta.init. If NULL, V is calculated from tune and an initial call to optim. See the Details section below for more information. Unless the log-posterior is expensive to compute it will typically be best to use the default V = NULL.
optim.method	The value of the method parameter sent to optim during an initial maximization of fun. See optim for more details.
optim.lower	The value of the lower parameter sent to optim during an initial maximization of fun. See optim for more details.
optim.upper	The value of the upper parameter sent to optim during an initial maximization of fun. See optim for more details.
optim.control	The value of the control parameter sent to optim during an initial maximization of fun. See optim for more details.
...	Additional arguments.

Details

MCMCmetrop1R produces a sample from a user-defined distribution using a random walk Metropolis algorithm with multivariate normal proposal distribution. See Gelman et al. (2003) and Robert & Casella (2004) for details of the random walk Metropolis algorithm.

The proposal distribution is centered at the current value of \theta and has variance-covariance V. If V is specified by the user to be NULL then V is calculated as: V = T (-1\cdot H)^{-1} T , where T is a the diagonal positive definite matrix formed from the tune and H is the approximate Hessian of fun evaluated at its mode. This last calculation is done via an initial call to optim.

Value

An mcmc object that contains the posterior sample. This object can be summarized by functions provided by the coda package.

References

Siddhartha Chib; Edward Greenberg. 1995. “Understanding the Metropolis-Hastings Algorithm." The American Statistician, 49, 327-335.

Andrew Gelman, John B. Carlin, Hal S. Stern, and Donald B. Rubin. 2003. Bayesian Data Analysis. 2nd Edition. Boca Raton: Chapman & Hall/CRC.

Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. “MCMCpack: Markov Chain Monte Carlo in R.”, Journal of Statistical Software. 42(9): 1-21. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.18637/jss.v042.i09")}.

Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. Scythe Statistical Library 1.0. http://scythe.wustl.edu.s3-website-us-east-1.amazonaws.com/.

Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2006. “Output Analysis and Diagnostics for MCMC (CODA)”, R News. 6(1): 7-11. https://CRAN.R-project.org/doc/Rnews/Rnews_2006-1.pdf.

Christian P. Robert and George Casella. 2004. Monte Carlo Statistical Methods. 2nd Edition. New York: Springer.

See Also

plot.mcmc, summary.mcmc, optim, metrop

Examples

  ## Not run: 

    ## logistic regression with an improper uniform prior
    ## X and y are passed as args to MCMCmetrop1R

    logitfun <- function(beta, y, X){
      eta <- X %*% beta
      p <- 1.0/(1.0+exp(-eta))
      sum( y * log(p) + (1-y)*log(1-p) )
    }

    x1 <- rnorm(1000)
    x2 <- rnorm(1000)
    Xdata <- cbind(1,x1,x2)
    p <- exp(.5 - x1 + x2)/(1+exp(.5 - x1 + x2))
    yvector <- rbinom(1000, 1, p)

    post.samp <- MCMCmetrop1R(logitfun, theta.init=c(0,0,0),
                              X=Xdata, y=yvector,
                              thin=1, mcmc=40000, burnin=500,
                              tune=c(1.5, 1.5, 1.5),
                              verbose=500, logfun=TRUE)

    raftery.diag(post.samp)
    plot(post.samp)
    summary(post.samp)
    ## ##################################################


    ##  negative binomial regression with an improper unform prior
    ## X and y are passed as args to MCMCmetrop1R
    negbinfun <- function(theta, y, X){
      k <- length(theta)
      beta <- theta[1:(k-1)]
      alpha <- exp(theta[k])
      mu <- exp(X %*% beta)
      log.like <- sum(
                      lgamma(y+alpha) - lfactorial(y) - lgamma(alpha) +
                      alpha * log(alpha/(alpha+mu)) +
                      y * log(mu/(alpha+mu))
                     )
    }

    n <- 1000
    x1 <- rnorm(n)
    x2 <- rnorm(n)
    XX <- cbind(1,x1,x2)
    mu <- exp(1.5+x1+2*x2)*rgamma(n,1)
    yy <- rpois(n, mu)

    post.samp <- MCMCmetrop1R(negbinfun, theta.init=c(0,0,0,0), y=yy, X=XX,
                              thin=1, mcmc=35000, burnin=1000,
                              tune=1.5, verbose=500, logfun=TRUE,
                              seed=list(NA,1))
    raftery.diag(post.samp)
    plot(post.samp)
    summary(post.samp)
    ## ##################################################


    ## sample from a univariate normal distribution with
    ## mean 5 and standard deviation 0.1
    ##
    ## (MCMC obviously not necessary here and this should
    ##  really be done with the logdensity for better
    ##  numerical accuracy-- this is just an illustration of how
    ##  MCMCmetrop1R works with a density rather than logdensity)

    post.samp <- MCMCmetrop1R(dnorm, theta.init=5.3, mean=5, sd=0.1,
                          thin=1, mcmc=50000, burnin=500,
                          tune=2.0, verbose=5000, logfun=FALSE)

    summary(post.samp)

  
## End(Not run)

MCMCpack documentation built on Sept. 11, 2024, 8:13 p.m.