MCMCmetrop1R  R Documentation 
This function allows a user to construct a sample from a userdefined 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 burnin 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 
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 nonnegative Hessian at the mode is often an
indication that the function of interest is not a proper density. Thus,

V 
The variancecovariance 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 userdefined 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 variancecovariance 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
.
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 MetropolisHastings Algorithm." The American Statistician, 49, 327335.
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): 121. \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.lsa.umich.edu.
Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2006. “Output Analysis and Diagnostics for MCMC (CODA)”, R News. 6(1): 711. https://CRAN.Rproject.org/doc/Rnews/Rnews_20061.pdf.
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) + (1y)*log(1p) )
}
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:(k1)]
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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.