MCMCmetrop1R {MCMCpack} | R Documentation |
This function allows a user to construct a sample from a user-defined continuous distribution using a random walk Metropolis algorithm.
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), ...)
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 |
theta.init |
Starting values for the sampling. Must be of the
appropriate dimension. It must also be the case that
|
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 |
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 |
logfun |
Logical indicating whether |
force.samp |
Logical indicating whether the sampling should
proceed if the Hessian calculated from the initial call to
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, |
V |
The variance-covariance matrix for the Gaussian proposal
distribution. Must be a square matrix or |
optim.method |
The value of the |
optim.lower |
The value of the |
optim.upper |
The value of the |
optim.control |
The value of the |
... |
Additional arguments. |
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*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
.
An mcmc object that contains the posterior sample. This object can be summarized by functions provided by the coda package.
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. http://www.jstatsoft.org/v42/i09/.
Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. Scythe Statistical Library 1.0. http://scythe.wustl.edu.
Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2002. Output Analysis and Diagnostics for MCMC (CODA). http://www-fis.iarc.fr/coda/.
Christian P. Robert and George Casella. 2004. Monte Carlo Statistical Methods. 2nd Edition. New York: Springer.
plot.mcmc
,
summary.mcmc
, optim
,
metrop
## 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)