[r-cran-mcmcpack] 09/90: Imported Upstream version 0.6-6
Andreas Tille
tille at debian.org
Fri Dec 16 09:07:33 UTC 2016
This is an automated email from the git hooks/post-receive script.
tille pushed a commit to branch master
in repository r-cran-mcmcpack.
commit 26a6f65ad227f871ea64d7d516247250819ae234
Author: Andreas Tille <tille at debian.org>
Date: Fri Dec 16 08:07:03 2016 +0100
Imported Upstream version 0.6-6
---
DESCRIPTION | 6 +-
HISTORY | 14 +++
INDEX | 11 ++-
R/MCMClogit.R | 102 +++++++++++++++++++---
R/MCMCmetrop1R.R | 45 ++++++----
R/hidden.R | 2 +
README | 9 +-
man/MCMCdynamicEI.Rd | 4 +-
man/MCMCfactanal.Rd | 21 ++---
man/MCMChierEI.Rd | 4 +-
man/MCMCirt1d.Rd | 13 +--
man/MCMCirtKd.Rd | 4 +-
man/MCMClogit.Rd | 128 ++++++++++++++++++++-------
man/MCMCmetrop1R.Rd | 98 +++++++++------------
man/MCMCmixfactanal.Rd | 15 ++--
man/MCMCmnl.Rd | 22 ++---
man/MCMCoprobit.Rd | 18 ++--
man/MCMCordfactanal.Rd | 12 +--
man/MCMCpanel.Rd | 10 +--
man/MCMCpoisson.Rd | 21 ++---
man/MCMCprobit.Rd | 12 +--
man/MCMCregress.Rd | 12 +--
man/MCMCtobit.Rd | 12 +--
src/MCMCdynamicEI.cc | 2 +-
src/MCMChierEI.cc | 2 +-
src/MCMClogituserprior.cc | 215 ++++++++++++++++++++++++++++++++++++++++++++++
src/MCMCmnlslice.cc | 2 +-
src/la.cc | 43 +++++++---
src/stat.cc | 8 +-
29 files changed, 627 insertions(+), 240 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index 989de0e..d1f2f5b 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
Package: MCMCpack
-Version: 0.6-5
-Date: 2005-6-29
+Version: 0.6-6
+Date: 2005-11-14
Title: Markov chain Monte Carlo (MCMC) Package
Author: Andrew D. Martin <admartin at wustl.edu>, and
Kevin M. Quinn <kevin_quinn at harvard.edu>
@@ -17,4 +17,4 @@ Description: This package contains functions to perform Bayesian
sampling algorithm, and tools for visualization.
License: GPL version 2 or newer
URL: http://mcmcpack.wustl.edu
-Packaged: Wed Jun 29 17:57:40 2005; adm
+Packaged: Mon Nov 14 09:56:45 2005; adm
diff --git a/HISTORY b/HISTORY
index 271a231..698f2d6 100644
--- a/HISTORY
+++ b/HISTORY
@@ -2,6 +2,20 @@
// Changes and Bug Fixes
//
+0.6-5 to 0.6-6
+ * fixed the std::accumulate problem pointed out by Kurt Hornik. Thanks
+ to Dan Pemstein for tracking the problem down and making the fix.
+ * fixed up how the force.samp option works with MCMCmetrop1R-- previously
+ if a diagonal element of the Hessian was 0 the correction wouldn't work.
+ * changed how additional arguments are passed to user-defined function
+ in MCMCmetrop1R. This breaks old code but provides a more standard
+ interface.
+ * when logfun==FALSE in MCMCmetrop1R the initial call to optim() now
+ maximizes the log of the user-defined function
+ * cleaned up the *.Rd files for the model fitting functions a bit
+ (more still needs to be done)
+ * MCMClogit now optionally accepts a user-defined prior density
+
0.6-4 to 0.6-5
* added a check so one can run MCMCirt1d without passing constraints
* fixed the Senate dataset so there are no duplicate names in the
diff --git a/INDEX b/INDEX
index a6eb58c..f99695c 100644
--- a/INDEX
+++ b/INDEX
@@ -9,8 +9,8 @@ MCMChierEI Markov Chain Monte Carlo for Wakefield's
Hierarchial Ecological Inference Model
MCMCirt1d Markov Chain Monte Carlo for One Dimensional
Item Response Theory Model
-MCMCirtKd Markov Chain Monte Carlo for K-Dimensional
- Item Response Theory Model
+MCMCirtKd Markov Chain Monte Carlo for K-Dimensional Item
+ Response Theory Model
MCMClogit Markov Chain Monte Carlo for Logistic
Regression
MCMCmetrop1R Metropolis Sampling from User-Written R
@@ -23,10 +23,9 @@ MCMCoprobit Markov Chain Monte Carlo for Ordered Probit
Regression
MCMCordfactanal Markov Chain Monte Carlo for Ordinal Data
Factor Analysis Model
-MCMCpanel Markov Chain Monte Carlo for the General
- Linear Panel Model
-MCMCpoisson Markov Chain Monte Carlo for Poisson
- Regression
+MCMCpanel Markov Chain Monte Carlo for the General Linear
+ Panel Model
+MCMCpoisson Markov Chain Monte Carlo for Poisson Regression
MCMCprobit Markov Chain Monte Carlo for Probit Regression
MCMCregress Markov Chain Monte Carlo for Gaussian Linear
Regression
diff --git a/R/MCMClogit.R b/R/MCMClogit.R
index 8aabcd9..ca2fca2 100644
--- a/R/MCMClogit.R
+++ b/R/MCMClogit.R
@@ -6,11 +6,13 @@
## Modified to meet new developer specification 7/15/2004 KQ
## Modified for new Scythe and rngs 7/25/2004 KQ
## note: B0 is now a precision
+## Modified to allow user-specified prior density 8/17/2005 KQ
+##
"MCMClogit" <-
function(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
thin=1, tune=1.1, verbose = 0, seed = NA, beta.start = NA,
- b0 = 0, B0 = 0, ...) {
+ b0 = 0, B0 = 0, user.prior.density=NULL, logfun=TRUE, ...) {
## checks
check.offset(list(...))
@@ -35,6 +37,50 @@
b0 <- mvn.prior[[1]]
B0 <- mvn.prior[[2]]
+ ## setup the environment so that fun can see the things passed as ...
+ userfun <- function(ttt) user.prior.density(ttt, ...)
+ my.env <- environment(fun = userfun)
+
+
+
+ ## check to make sure user.prior.density returns a numeric scalar and
+ ## starting values have positive prior mass
+ if (is.function(user.prior.density)){
+ funval <- userfun(beta.start)
+ if (length(funval) != 1){
+ cat("user.prior.density does not return a scalar.\n")
+ stop("Respecify and call MCMClogit() again. \n")
+ }
+ if (!is.numeric(funval)){
+ cat("user.prior.density does not return a numeric value.\n")
+ stop("Respecify and call MCMClogit() again. \n")
+ }
+
+ if (identical(funval, Inf)){
+ cat("user.prior.density(beta.start) == Inf.\n")
+ stop("Respecify and call MCMClogit() again. \n")
+ }
+
+ if (logfun){
+ if (identical(funval, -Inf)){
+ cat("user.prior.density(beta.start) == -Inf.\n")
+ stop("Respecify and call MCMClogit() again. \n")
+ }
+ }
+ else{
+ if (funval <= 0){
+ cat("user.prior.density(beta.start) <= 0.\n")
+ stop("Respecify and call MCMClogit() again. \n")
+ }
+ }
+ }
+ else if (!is.null(user.prior.density)){
+ cat("user.prior.density is neither a NULL nor a function.\n")
+ stop("Respecify and call MCMClogit() again. \n")
+ }
+
+
+
## form the tuning parameter
tune <- vector.tune(tune, K)
V <- vcov(glm(formula=formula, data=data, family=binomial))
@@ -46,22 +92,50 @@
}
- ## define holder for posterior density sample
- sample <- matrix(data=0, mcmc/thin, dim(X)[2] )
+ propvar <- tune %*% V %*% tune
+
+
## call C++ code to draw sample
- auto.Scythe.call(output.object="posterior", cc.fun.name="MCMClogit",
- sample=sample, Y=Y, X=X, burnin=as.integer(burnin),
- mcmc=as.integer(mcmc), thin=as.integer(thin),
- tune=tune, lecuyer=as.integer(lecuyer),
- seedarray=as.integer(seed.array),
- lecuyerstream=as.integer(lecuyer.stream),
- verbose=as.integer(verbose), betastart=beta.start,
- b0=b0, B0=B0, V=V)
+ if (is.null(user.prior.density)){
+ ## define holder for posterior density sample
+ sample <- matrix(data=0, mcmc/thin, dim(X)[2] )
+
+ auto.Scythe.call(output.object="posterior", cc.fun.name="MCMClogit",
+ sample=sample, Y=Y, X=X, burnin=as.integer(burnin),
+ mcmc=as.integer(mcmc), thin=as.integer(thin),
+ tune=tune, lecuyer=as.integer(lecuyer),
+ seedarray=as.integer(seed.array),
+ lecuyerstream=as.integer(lecuyer.stream),
+ verbose=as.integer(verbose), betastart=beta.start,
+ b0=b0, B0=B0, V=V)
+
+ ## put together matrix and build MCMC object to return
+ output <- form.mcmc.object(posterior, names=xnames,
+ title="MCMClogit Posterior Density Sample")
+
+
+ }
+ else {
+ sample <- .Call("MCMClogituserprior_cc",
+ userfun, as.integer(Y), as.matrix(X),
+ as.double(beta.start),
+ my.env, as.integer(burnin), as.integer(mcmc),
+ as.integer(thin),
+ as.integer(verbose),
+ lecuyer=as.integer(lecuyer),
+ seedarray=as.integer(seed.array),
+ lecuyerstream=as.integer(lecuyer.stream),
+ as.logical(logfun),
+ as.matrix(propvar),
+ PACKAGE="MCMCpack")
+
+ output <- mcmc(data=sample, start=burnin+1,
+ end=burnin+mcmc, thin=thin)
+ varnames(output) <- as.list(xnames)
+ attr(output, "title") <- "MCMClogit Posterior Density Sample"
+ }
- ## put together matrix and build MCMC object to return
- output <- form.mcmc.object(posterior, names=xnames,
- title="MCMClogit Posterior Density Sample")
return(output)
}
diff --git a/R/MCMCmetrop1R.R b/R/MCMCmetrop1R.R
index 9380859..11d1dda 100644
--- a/R/MCMCmetrop1R.R
+++ b/R/MCMCmetrop1R.R
@@ -2,8 +2,12 @@
## random walk Metropolis algorithm
##
## KQ 6/24/2004
+##
## modified to work with non-invertible Hessian KQ 6/28/2005
##
+## changed the method used to pass additional arguments to the user-defined
+## function KQ 8/15/2005
+##
"MCMCmetrop1R" <- function(fun, theta.init,
burnin=500, mcmc=20000, thin=1,
@@ -26,24 +30,28 @@
## setup the environment so that fun can see the things passed as ...
- ## it should be the case that users can specify arguments with the same
- ## names as variables defined in MCMCmetrop1R without causing problems.
- my.env <- new.env()
- environment(fun) <- my.env
- dots <- list(...)
- dotnames <- names(dots)
- ndots <- length(dots)
- if (ndots >= 1){
- for (i in 1:ndots){
- assign(x=dotnames[[i]], value=dots[[i]], inherits=FALSE, envir=my.env)
- }
+ userfun <- function(ttt) fun(ttt, ...)
+ my.env <- environment(fun = userfun)
+
+
+ ## setup function for maximization based on value of logfun
+ if (logfun){
+ maxfun <- fun
+ }
+ else if (logfun==FALSE){
+ maxfun <- function(ttt, ...) log(fun(ttt, ...))
+ }
+ else{
+ cat("logfun not a logical value.\n")
+ stop("Respecifiy and call MCMCmetrop1R() again. \n",
+ call.=FALSE)
}
## find approx mode and Hessian using optim()
- opt.out <- optim(theta.init, fun,
+ opt.out <- optim(theta.init, maxfun,
control=list(fnscale=-1, trace=optim.trace,
REPORT=optim.REPORT, maxit=optim.maxit),
- method="BFGS", hessian=TRUE)
+ method="BFGS", hessian=TRUE, ...)
if(opt.out$convergence!=0){
warning("Mode and Hessian were not found with call to optim().\nSampling proceeded anyway. \n")
}
@@ -54,6 +62,13 @@
hess.new <- opt.out$hessian
hess.flag <- 0
if (force.samp==TRUE){
+ if (max(diag(optim.out$hessian)==0)){
+ for (i in 1:nrow(hess.new)){
+ if (hess.new[i,i] == 0){
+ hess.new[i,i] <- -1e-6
+ }
+ }
+ }
while (is.null(CC)){
hess.flag <- 1
hess.new <- hess.new - diag(diag(0.01 * abs(opt.out$hessian)))
@@ -76,12 +91,10 @@
}
propvar <- tune %*% solve(-1*hess.new) %*% tune
- ## the old way that only worked when Hessian was neg. def.
- ##propvar <- tune %*% solve(-1*opt.out$hessian) %*% tune
## Call the C++ function to do the MCMC sampling
- sample <- .Call("MCMCmetrop1R_cc", fun, as.double(theta.init),
+ sample <- .Call("MCMCmetrop1R_cc", userfun, as.double(theta.init),
my.env, as.integer(burnin), as.integer(mcmc),
as.integer(thin),
as.integer(verbose),
diff --git a/R/hidden.R b/R/hidden.R
index 647507d..486a6ad 100644
--- a/R/hidden.R
+++ b/R/hidden.R
@@ -9,6 +9,8 @@
## create an agreement score matrix from a vote matrix
## subjects initially on rows and items on cols of X
+## note: treats missing votes as category for agreement / might be
+## more principled to treat them in another fashion
"agree.mat" <- function(X){
X <- t(X) # put subjects on columns
n <- ncol(X)
diff --git a/README b/README
index eca6e6c..2589dc0 100644
--- a/README
+++ b/README
@@ -1,8 +1,8 @@
//
-// MCMCpack, Version 0.6-1
+// MCMCpack, Version 0.6-6
//
-Release Date: February 18, 2005
+Release Date: November 14, 2005
This package contains functions for posterior simulation for a number of
statistical models. All simulation is done in compiled C++ written in
@@ -10,8 +10,7 @@ the Scythe Statistical Library Version 1.0. All models return coda mcmc
objects that can then be summarized using coda functions or the coda
menu interface. The package also contains some useful utility
functions, including some additional density functions and pseudo-random
-number generators for statistical distributions, a general purpose Metropolis
-sampling algorithm, and tools for visualization.
+number generators for statistical distributions, a general purpose Metropolis sampling algorithm, and tools for visualization.
//
// Authors
@@ -26,7 +25,7 @@ Kevin M. Quinn <kevin_quinn at harvard.edu>
This package (along with Scythe) uses C++ and the Standard Template
Library (STL). It requires use of the GCC compiler 3.0 or greater. It
-has been tested on GCC 3.1, 3.2, 3.3, and 3.4 on Linux and MacOS X.
+has been tested on GCC 3.1, 3.2, 3.3, 3.4, and 4.0 on Linux and MacOS X.
We have also successfully compiled this package for Windows using
a MinGW cross-compiler using binutils 2.13, gcc 3.2, mingw
diff --git a/man/MCMCdynamicEI.Rd b/man/MCMCdynamicEI.Rd
index 1c35988..d549a41 100644
--- a/man/MCMCdynamicEI.Rd
+++ b/man/MCMCdynamicEI.Rd
@@ -74,7 +74,7 @@ MCMCdynamicEI(r0, r1, c0, c1, burnin=5000, mcmc=50000, thin=1,
}
\value{
- An mcmc object that contains the posterior density sample.
+ An mcmc object that contains the sample from the posterior distribution.
This object can be summarized by functions provided by the coda package.
}
@@ -129,7 +129,7 @@ MCMCdynamicEI(r0, r1, c0, c1, burnin=5000, mcmc=50000, thin=1,
Inference centers on \eqn{p_0}{p0}, \eqn{p_1}{p1},
\eqn{\sigma^2_0}{sigma^2_0}, and \eqn{\sigma^2_1}{sigma^2_1}.
Univariate slice sampling (Neal, 2003) together with Gibbs sampling
- is used to sample from the posterior density.
+ is used to sample from the posterior distribution.
}
\references{
diff --git a/man/MCMCfactanal.Rd b/man/MCMCfactanal.Rd
index 3a00c4b..fa866f7 100644
--- a/man/MCMCfactanal.Rd
+++ b/man/MCMCfactanal.Rd
@@ -2,12 +2,12 @@
\alias{MCMCfactanal}
\title{Markov Chain Monte Carlo for Normal Theory Factor Analysis Model}
\description{
- This function generates a posterior density sample from Normal theory
- factor analysis model. Normal priors are assumed on the factor
- loadings and factor scores while inverse Gamma priors are assumed for
- the uniquenesses. The user supplies data and parameters for the prior
- distributions, and a sample from the posterior density is returned as
- an mcmc object, which can be subsequently analyzed with
+ This function generates a sample from the posterior distribution of
+ a normal theory factor analysis model. Normal priors are assumed on
+ the factor loadings and factor scores while inverse Gamma priors are
+ assumed for the uniquenesses. The user supplies data and parameters
+ for the prior distributions, and a sample from the posterior distribution
+ is returned as an mcmc object, which can be subsequently analyzed with
functions provided in the coda package.
}
@@ -114,8 +114,9 @@ MCMCfactanal(x, factors, lambda.constraints=list(),
}
\value{
- An mcmc object that contains the posterior density sample. This
- object can be summarized by functions provided by the coda package.
+ An mcmc object that contains the sample from the posterior
+ distribution. This object can be summarized by functions provided by
+ the coda package.
}
\details{The model takes the following form:
@@ -145,11 +146,11 @@ MCMCfactanal(x, factors, lambda.constraints=list(),
\deqn{\Psi_{ii} \sim \mathcal{IG}(a_{0_i}/2, b_{0_i}/2),
i=1,\ldots,k}{Psi_ii ~ IG(a0_i/2, b0_i/2), i=1,...,k}
- \code{MCMCfactanal} simulates from the posterior density using
+ \code{MCMCfactanal} simulates from the posterior distribution using
standard Gibbs sampling. The simulation proper is done in
compiled C++ code to maximize efficiency. Please consult the
coda documentation for a comprehensive list of functions that
- can be used to analyze the posterior density sample.
+ can be used to analyze the posterior sample.
As is the case with all measurement models, make sure that you have plenty
of free memory, especially when storing the scores.
diff --git a/man/MCMChierEI.Rd b/man/MCMChierEI.Rd
index 47918b2..34a8102 100644
--- a/man/MCMChierEI.Rd
+++ b/man/MCMChierEI.Rd
@@ -73,7 +73,7 @@ MCMChierEI(r0, r1, c0, c1, burnin=5000, mcmc=50000, thin=1,
}
\value{
- An mcmc object that contains the posterior density sample.
+ An mcmc object that contains the sample from the posterior distribution.
This object can be summarized by functions provided by the coda package.
}
@@ -129,7 +129,7 @@ MCMChierEI(r0, r1, c0, c1, burnin=5000, mcmc=50000, thin=1,
\eqn{\mu_1}{mu1}, \eqn{\sigma^2_0}{sigma^2_0}, and
\eqn{\sigma^2_1}{sigma^2_1}.
Univariate slice sampling (Neal, 2003) along with Gibbs sampling is
- used to sample from the posterior density.
+ used to sample from the posterior distribution.
See Section 5.4 of Wakefield (2003) for discussion of the priors used
here. \code{MCMChierEI} departs from the Wakefield model in that the
diff --git a/man/MCMCirt1d.Rd b/man/MCMCirt1d.Rd
index 3f1db55..9077a77 100644
--- a/man/MCMCirt1d.Rd
+++ b/man/MCMCirt1d.Rd
@@ -3,12 +3,12 @@
\title{Markov Chain Monte Carlo for One Dimensional Item Response Theory
Model}
\description{
- This function generates a posterior density sample from a one
+ This function generates a sample from the posterior distribution of a one
dimentional item response theory (IRT) model, with Normal
priors on the subject abilities (ideal points), and
multivariate Normal priors on the item parameters. The user
supplies data and priors, and a sample from the posterior
- density is returned as an mcmc object, which can be
+ distribution is returned as an mcmc object, which can be
subsequently analyzed with functions provided in the coda
package.
@@ -121,19 +121,20 @@ MCMCirt1d(datamatrix, theta.constraints=list(), burnin = 1000,
}
\value{
- An mcmc object that contains the posterior density sample. This
- object can be summarized by functions provided by the coda package.
+ An mcmc object that contains the sample from the posterior
+ distribution. This object can be summarized by functions provided by
+ the coda package.
}
\details{
- \code{MCMCirt1d} simulates from the posterior density using
+ \code{MCMCirt1d} simulates from the posterior distribution using
standard Gibbs sampling using data augmentation (a Normal draw
for the subject abilities, a multivariate Normal
draw for the item parameters, and a truncated Normal draw for
the latent utilities). The simulation proper is done in
compiled C++ code to maximize efficiency. Please consult the
coda documentation for a comprehensive list of functions that
- can be used to analyze the posterior density sample.
+ can be used to analyze the posterior sample.
The model takes the following form. We assume that each
subject has an subject ability (ideal point) denoted
diff --git a/man/MCMCirtKd.Rd b/man/MCMCirtKd.Rd
index 20de18c..3ff5294 100644
--- a/man/MCMCirtKd.Rd
+++ b/man/MCMCirtKd.Rd
@@ -3,7 +3,7 @@
\title{Markov Chain Monte Carlo for K-Dimensional Item Response Theory
Model}
\description{
- This function generates a posterior density sample from a
+ This function generates a sample from the posterior distribution of a
K-dimensional item response theory (IRT) model, with standard
normal priors on the subject abilities (ideal points), and
normal priors on the item parameters. The user
@@ -128,7 +128,7 @@ MCMCirtKd(datamatrix, dimensions, item.constraints=list(),
the latent utilities). The simulation proper is done in
compiled C++ code to maximize efficiency. Please consult the
coda documentation for a comprehensive list of functions that
- can be used to analyze the posterior density sample.
+ can be used to analyze the posterior sample.
The default number of burnin and mcmc iterations is much
smaller than the typical default values in MCMCpack. This is
diff --git a/man/MCMClogit.Rd b/man/MCMClogit.Rd
index dda7d63..c0ab626 100644
--- a/man/MCMClogit.Rd
+++ b/man/MCMClogit.Rd
@@ -2,10 +2,10 @@
\alias{MCMClogit}
\title{Markov Chain Monte Carlo for Logistic Regression}
\description{
- This function generates a posterior density sample
- from a logistic regression model using a random walk Metropolis
+ This function generates a sample from the posterior distribution of
+ a logistic regression model using a random walk Metropolis
algorithm. The user supplies data and priors,
- and a sample from the posterior density is returned as an mcmc
+ and a sample from the posterior distribution is returned as an mcmc
object, which can be subsequently analyzed with functions
provided in the coda package.
}
@@ -13,7 +13,7 @@
\usage{
MCMClogit(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
thin=1, tune=1.1, verbose = 0, seed = NA, beta.start = NA,
- b0 = 0, B0 = 0, ...) }
+ b0 = 0, B0 = 0, user.prior.density=NULL, logfun=TRUE, ...) }
\arguments{
\item{formula}{Model formula.}
@@ -31,7 +31,7 @@ MCMClogit(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
scalar or a \eqn{k}{k}-vector, where \eqn{k}{k} is the length of
\eqn{\beta}{beta}.Make sure that the
acceptance rate is satisfactory (typically between 0.20 and 0.5)
- before using the posterior density sample for inference.}
+ before using the posterior sample for inference.}
\item{verbose}{A switch which determines whether or not the progress of
the sampler is printed to the screen. If \code{verbose} is greater
@@ -57,48 +57,75 @@ MCMClogit(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
use the maximum likelihood estimate of \eqn{\beta}{beta} as the starting
value.}
- \item{b0}{The prior mean of \eqn{\beta}{beta}. This can either be a
- scalar or a column
- vector with dimension equal to the number of betas. If this takes a scalar
- value, then that value will serve as the prior mean for all of the
- betas.}
+ \item{b0}{If \code{user.prior.density==NULL} \code{b0} is the prior
+ mean of \eqn{\beta}{beta} under a multivariate normal prior. This
+ can either be a scalar or a column vector with dimension equal to
+ the number of betas. If this takes a scalar value, then that value
+ will serve as the prior mean for all of the betas.}
- \item{B0}{The prior precision of \eqn{\beta}{beta}. This can either be a
- scalar
- or a square matrix with dimensions equal to the number of betas. If this
- takes a scalar value, then that value times an identity matrix serves
- as the prior precision of \eqn{\beta}{beta}. Default value of 0 is
- equivalent to an improper uniform prior for beta.}
-
+ \item{B0}{If \code{user.prior.density==NULL} \code{B0} is the prior
+ precision of \eqn{\beta}{beta} under a multivariate normal prior.
+ This can either be a scalar or a square matrix with dimensions
+ equal to the number of betas. If this takes a scalar value, then
+ that value times an identity matrix serves as the prior precision
+ of \eqn{\beta}{beta}. Default value of 0 is equivalent to an
+ improper uniform prior for beta.}
+
+ \item{user.prior.density}{If non-NULL, the prior (log)density up to
+ a constant of proportionality. 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 prior (log)density is to be evaluated. Additional
+ arguments can be passed to \code{user.prior.density()} by
+ inserting them in the call to \code{MCMClogit()}. See the Details
+ section and the examples below for more information. }
+
+ \item{logfun}{Logical indicating whether \code{use.prior.density()}
+ returns the natural log of a density function (TRUE) or a density
+ (FALSE).}
+
\item{\ldots}{further arguments to be passed}
}
\value{
- An mcmc object that contains the posterior density sample. This
+ An mcmc object that contains the posterior sample. This
object can be summarized by functions provided by the coda package.
}
-\details{\code{MCMClogit} simulates from the posterior density of a logistic
- regression model using a random walk Metropolis algorithm. The simulation
- proper is done in compiled C++ code to maximize efficiency. Please consult
- the coda documentation for a comprehensive list of functions that can be
- used to analyze the posterior density sample.
+\details{\code{MCMClogit} simulates from the posterior distribution of a
+ logistic regression model using a random walk Metropolis
+ algorithm. The simulation proper is done in compiled C++ code to
+ maximize efficiency. Please consult the coda documentation for a
+ comprehensive list of functions that can be used to analyze the
+ posterior sample.
The model takes the following form:
\deqn{y_i \sim \mathcal{B}ernoulli(\pi_i)}{y_i ~ Bernoulli(pi_i)}
Where the inverse link function:
\deqn{\pi_i = \frac{\exp(x_i'\beta)}{1 + \exp(x_i'\beta)}}{pi_i =
exp(x_i'beta) / [1 + exp(x_i'beta)]}
- We assume a multivariate Normal prior on \eqn{\beta}{beta}:
- \deqn{\beta \sim \mathcal{N}(b_0,B_0^{-1})}{beta ~ N(b0,B0^(-1))}
-
- The Metropollis proposal distribution is centered at the current value of
- \eqn{\theta}{theta} and has variance-covariance \eqn{V = T
- (B_0 + C^{-1})^{-1} T }{V = T (B0 + C^{-1})^{-1} T}, where
- \eqn{T}{T} is a the diagonal positive definite matrix formed from the
- \code{tune}, \eqn{B_0}{B0} is the prior precision, and \eqn{C}{C} is
- the large sample variance-covariance matrix of the MLEs. This last
- calculation is done via an initial call to \code{glm}.
+ By default, we assume a multivariate Normal prior on \eqn{\beta}{beta}:
+ \deqn{\beta \sim \mathcal{N}(b_0,B_0^{-1})}{beta ~ N(b0,B0^(-1))}
+ Additionally, arbitrary user-defined priors can be specified with the
+ \code{user.prior.density} argument.
+
+ If the default multivariate normal prior is used, the Metropollis
+ proposal distribution is centered at the current value of
+ \eqn{\beta}{beta} and has variance-covariance \eqn{V = T (B_0 +
+ C^{-1})^{-1} T }{V = T (B0 + C^{-1})^{-1} T}, where \eqn{T}{T} is a
+ the diagonal positive definite matrix formed from the \code{tune},
+ \eqn{B_0}{B0} is the prior precision, and \eqn{C}{C} is the large
+ sample variance-covariance matrix of the MLEs. This last calculation
+ is done via an initial call to \code{glm}.
+
+ If a user-defined prior is used, the Metropollis
+ proposal distribution is centered at the current value of
+ \eqn{\beta}{beta} and has variance-covariance \eqn{V = T
+ C T }{V = T C T}, where \eqn{T}{T} is a
+ the diagonal positive definite matrix formed from the \code{tune} and
+ \eqn{C}{C} is the large sample variance-covariance matrix of the
+ MLEs. This last calculation is done via an initial call to
+ \code{glm}.
}
\references{
@@ -114,10 +141,45 @@ MCMClogit(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
\examples{
\dontrun{
+ ## default improper uniform prior
data(birthwt)
posterior <- MCMClogit(low~age+as.factor(race)+smoke, data=birthwt)
plot(posterior)
summary(posterior)
+
+
+ ## multivariate normal prior
+ data(birthwt)
+ posterior <- MCMClogit(low~age+as.factor(race)+smoke, b0=0, B0=.001,
+ data=birthwt)
+ plot(posterior)
+ summary(posterior)
+
+
+ ## user-defined independent Cauchy prior
+ logpriorfun <- function(beta){
+ sum(dcauchy(beta, log=TRUE))
+ }
+
+ posterior <- MCMClogit(low~age+as.factor(race)+smoke,
+ data=birthwt, user.prior.density=logpriorfun,
+ logfun=TRUE)
+ plot(posterior)
+ summary(posterior)
+
+
+ ## user-defined independent Cauchy prior with additional args
+ logpriorfun <- function(beta, location, scale){
+ sum(dcauchy(beta, location, scale, log=TRUE))
+ }
+
+ posterior <- MCMClogit(low~age+as.factor(race)+smoke,
+ data=birthwt, user.prior.density=logpriorfun,
+ logfun=TRUE, location=0, scale=10)
+ plot(posterior)
+ summary(posterior)
+
+
}
}
diff --git a/man/MCMCmetrop1R.Rd b/man/MCMCmetrop1R.Rd
index d1ceadd..6e88a8c 100644
--- a/man/MCMCmetrop1R.Rd
+++ b/man/MCMCmetrop1R.Rd
@@ -3,8 +3,7 @@
\title{Metropolis Sampling from User-Written R function}
\description{
This function allows a user to construct a sample from a user-defined
- R function using a random walk Metropolis algorithm. It assumes the
- parameters to be sampled are in a single block.
+ continuous distribution using a random walk Metropolis algorithm.
}
\usage{
MCMCmetrop1R(fun, theta.init, burnin = 500, mcmc = 20000, thin = 1,
@@ -13,14 +12,20 @@ MCMCmetrop1R(fun, theta.init, burnin = 500, mcmc = 20000, thin = 1,
optim.maxit = 500, ...)
}
\arguments{
- \item{fun}{The (log)density from which to take a sample. This should
- be a function defined in R and it should take a single
- argument. Additional arguments can be passed implicitly by either
- putting them in the global environment or by passing them as
- additional arguments to \code{MCMCmetrop1R()}. See the Details section
- and the examples below. }
+ \item{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 \code{fun()} by inserting them in the call to
+ \code{MCMCmetrop1R()}. See the Details section and the examples
+ below for more information.}
+
\item{theta.init}{Starting values for the sampling. Must be of the
- appropriate dimension.}
+ appropriate dimension. It must also be the case that
+ \code{fun(theta.init, ...)} is greater than \code{-Inf} if
+ \code{fun()} is a logdensity or greater than 0 if \code{fun()} is a
+ density.}
\item{burnin}{The number of burn-in iterations for the sampler.}
@@ -29,9 +34,9 @@ MCMCmetrop1R(fun, theta.init, burnin = 500, mcmc = 20000, thin = 1,
\item{thin}{The thinning interval used in the simulation. The number of
MCMC iterations must be divisible by this value.}
- \item{tune}{Can be either a positive scalar or a
- \eqn{k}{k}-vector, where \eqn{k}{k} is the length of
- \eqn{\theta}{theta}.}
+ \item{tune}{The tuning parameter for the Metropolis sampling.
+ Can be either a positive scalar or a \eqn{k}{k}-vector, where
+ \eqn{k}{k} is the length of \eqn{\theta}{theta}.}
\item{verbose}{A switch which determines whether or not the progress of
the sampler is printed to the screen. If \code{verbose} is greater
@@ -54,7 +59,7 @@ MCMCmetrop1R(fun, theta.init, burnin = 500, mcmc = 20000, thin = 1,
\item{force.samp}{Logical indicating whether the sampling should
proceed if the Hessian calculated from the initial call to
- \code{optim} routine to maximize the (log)posterior
+ \code{optim} routine to maximize the (log)density
is not negative definite. If \code{force.samp==TRUE}
and the Hessian from \code{optim} is non-negative definite, the
Hessian is rescaled by subtracting small values from it's main
@@ -64,9 +69,9 @@ MCMCmetrop1R(fun, theta.init, burnin = 500, mcmc = 20000, thin = 1,
\code{optim} is non-negative definite, an error message is
printed and the call to \code{MCMCmetrop1R} is terminated.
- \emph{Please note that a non-negative Hessian at the mode is a very
- good indication that the model being fitted is not
- identified. Thus, \code{force.samp} should only be set equal to
+ \emph{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, \code{force.samp} should only be set equal to
\code{TRUE} with great caution.}
}
@@ -82,7 +87,7 @@ MCMCmetrop1R(fun, theta.init, burnin = 500, mcmc = 20000, thin = 1,
\item{\dots}{Additional arguments.}
}
\details{
-MCMCmetrop1R produces a sample from a user-defined (log)density using a
+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.
@@ -94,15 +99,6 @@ diagonal positive definite matrix formed from the \code{tune} and
\eqn{H}{H} is the approximate Hessian of \code{fun} evaluated at it's
mode. This last calculation is done via an initial call to
\code{optim}.
-
-Passing data into a (log)density function does not conform to standard R
-practice. The (log)density function should have only one argument -- the
-parameter vector. The data used by the (log)density need to be either in the
-global environment or defined in the call to \code{MCMCmetrop1R()}, which will
-create an environment that contains the data and then call the (log)density function in that environment. The names of the data objects
-are not checked, so be careful that these match the names used with the
-(log)density function. \emph{NOTE:} This interface will likely be changed in
-a future implementation of MCMCpack.
}
\value{
An mcmc object that contains the posterior density sample. This
@@ -114,7 +110,7 @@ a future implementation of MCMCpack.
## logistic regression with an improper uniform prior
## X and y are passed as args to MCMCmetrop1R
- logitfun <- function(beta){
+ logitfun <- function(beta, y, X){
eta <- X \%*\% beta
p <- 1.0/(1.0+exp(-eta))
sum( y * log(p) + (1-y)*log(1-p) )
@@ -137,36 +133,10 @@ a future implementation of MCMCpack.
summary(post.samp)
## ##################################################
-
- ## logistic regression with an improper uniform prior
- ## X and y are now in the global environment
-
- logitfun <- function(beta){
- 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)
- X <- cbind(1,x1,x2)
- p <- exp(.5 - x1 + x2)/(1+exp(.5 - x1 + x2))
- y <- rbinom(1000, 1, p)
-
- post.samp <- MCMCmetrop1R(logitfun, theta.init=c(0,0,0),
- thin=1, mcmc=40000, burnin=500,
- tune=c(1.5, 1.5, 1.5),
- verbose=500, logfun=TRUE, optim.maxit=100)
-
- 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){
+ negbinfun <- function(theta, y, X){
k <- length(theta)
beta <- theta[1:(k-1)]
alpha <- exp(theta[k])
@@ -193,7 +163,22 @@ a future implementation of MCMCpack.
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)
+
}
}
\references{
@@ -212,6 +197,7 @@ a future implementation of MCMCpack.
Statistical Methods}. 2nd Edition. New York: Springer.
}
\seealso{\code{\link[coda]{plot.mcmc}},
- \code{\link[coda]{summary.mcmc}}, \code{\link[stats]{optim}}}
+ \code{\link[coda]{summary.mcmc}}, \code{\link[stats]{optim}},
+ \code{\link[mcmc]{metrop}}}
\keyword{models}
diff --git a/man/MCMCmixfactanal.Rd b/man/MCMCmixfactanal.Rd
index 7434d38..ed8aab3 100644
--- a/man/MCMCmixfactanal.Rd
+++ b/man/MCMCmixfactanal.Rd
@@ -2,13 +2,13 @@
\alias{MCMCmixfactanal}
\title{Markov Chain Monte Carlo for Mixed Data Factor Analysis Model}
\description{
- This function generates a posterior density sample from a mixed data
- (both continuous and ordinal) factor analysis model. Normal priors are
- assumed on the factor loadings and factor scores, improper
+ This function generates a sample from the posterior distribution of a
+ mixed data (both continuous and ordinal) factor analysis model. Normal
+ priors are assumed on the factor loadings and factor scores, improper
uniform priors are assumed on the cutpoints, and inverse gamma priors
are assumed for the error variances (uniquenesses). The user supplies
data and parameters for the prior distributions, and a sample from the
- posterior density is returned as an mcmc object, which can be
+ posterior distribution is returned as an mcmc object, which can be
subsequently analyzed with functions provided in the coda package.
}
@@ -143,7 +143,7 @@ MCMCmixfactanal(x, factors, lambda.constraints=list(),
}
\value{
- An mcmc object that contains the posterior density sample. This
+ An mcmc object that contains the posterior sample. This
object can be summarized by functions provided by the coda package.
}
@@ -198,7 +198,7 @@ MCMCmixfactanal(x, factors, lambda.constraints=list(),
i=1,\dots,n}{phi_i(2:d) ~ N(0, I),
i=1,...,n}
- \code{MCMCmixfactanal} simulates from the posterior density using
+ \code{MCMCmixfactanal} simulates from the posterior distribution using
a Metropolis-Hastings within Gibbs sampling algorithm. The algorithm
employed is based on work by Cowles (1996). Note that
the first element of \eqn{\phi_i}{phi_i} is a 1. As a result, the
@@ -208,8 +208,7 @@ MCMCmixfactanal(x, factors, lambda.constraints=list(),
returned in the mcmc object.
The simulation proper is done in compiled C++ code to maximize
efficiency. Please consult the coda documentation for a comprehensive
- list of functions that can be used to analyze the posterior density
- sample.
+ list of functions that can be used to analyze the posterior sample.
As is the case with all measurement models, make sure that you have plenty
of free memory, especially when storing the scores.
diff --git a/man/MCMCmnl.Rd b/man/MCMCmnl.Rd
index f9be68d..7562f8e 100644
--- a/man/MCMCmnl.Rd
+++ b/man/MCMCmnl.Rd
@@ -3,10 +3,10 @@
\title{Markov Chain Monte Carlo for Multinomial Logistic Regression}
\description{
- This function generates a posterior density sample
- from a multinomial logistic regression model using either a random
+ This function generates a sample from the posterior distribution of
+ a multinomial logistic regression model using either a random
walk Metropolis algorithm or a slice sampler. The user supplies data
- and priors, and a sample from the posterior density is returned as an
+ and priors, and a sample from the posterior distribution is returned as an
mcmc object, which can be subsequently analyzed with functions provided
in the coda package.
}
@@ -74,7 +74,7 @@ MCMCmnl(formula, baseline = NULL, data = parent.frame(),
scalar or a \eqn{k}{k}-vector, where \eqn{k}{k} is the length of
\eqn{\beta}{beta}. Make sure that the
acceptance rate is satisfactory (typically between 0.20 and 0.5)
- before using the posterior density sample for inference. }
+ before using the posterior sample for inference. }
\item{verbose}{A switch which determines whether or not the progress of
@@ -115,16 +115,16 @@ MCMCmnl(formula, baseline = NULL, data = parent.frame(),
\item{\dots}{Further arguments to be passed. }
}
\value{
- An mcmc object that contains the posterior density sample. This
+ An mcmc object that contains the posterior sample. This
object can be summarized by functions provided by the coda package.
}
\details{
-\code{MCMCmnl} simulates from the posterior density of a multinomial logistic
- regression model using either a random walk Metropolis algorithm or a
- univariate slice sampler. The simulation
- proper is done in compiled C++ code to maximize efficiency. Please consult
- the coda documentation for a comprehensive list of functions that can be
- used to analyze the posterior density sample.
+\code{MCMCmnl} simulates from the posterior distribution of a
+multinomial logistic regression model using either a random walk
+Metropolis algorithm or a univariate slice sampler. The simulation
+proper is done in compiled C++ code to maximize efficiency. Please
+consult the coda documentation for a comprehensive list of functions
+that can be used to analyze the posterior sample.
The model takes the following form: \deqn{y_i \sim
\mathcal{M}ultinomial(\pi_i)}{y_i ~ Multinomial(pi_i)}
diff --git a/man/MCMCoprobit.Rd b/man/MCMCoprobit.Rd
index b1a4a04..813c575 100644
--- a/man/MCMCoprobit.Rd
+++ b/man/MCMCoprobit.Rd
@@ -2,10 +2,10 @@
\alias{MCMCoprobit}
\title{Markov Chain Monte Carlo for Ordered Probit Regression}
\description{
- This function generates a posterior density sample
- from an ordered probit regression model using the data augmentation
+ This function generates a sample from the posterior distribution of
+ an ordered probit regression model using the data augmentation
approach of Cowles (1996). The user supplies data and priors,
- and a sample from the posterior density is returned as an mcmc
+ and a sample from the posterior distribution is returned as an mcmc
object, which can be subsequently analyzed with functions
provided in the coda package.
}
@@ -70,16 +70,16 @@ MCMCoprobit(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
}
\value{
- An mcmc object that contains the posterior density sample. This
+ An mcmc object that contains the posterior sample. This
object can be summarized by functions provided by the coda package.
}
\details{
-\code{MCMCoprobit} simulates from the posterior density of a ordered probit
- regression model using data augmentation. The simulation proper is
- done in compiled C++ code to maximize efficiency. Please consult the
- coda documentation for a comprehensive list of functions that can be
- used to analyze the posterior density sample.
+\code{MCMCoprobit} simulates from the posterior distribution of a
+ordered probit regression model using data augmentation. The simulation
+proper is done in compiled C++ code to maximize efficiency. Please
+consult the coda documentation for a comprehensive list of functions
+that can be used to analyze the posterior sample.
The observed variable \eqn{y_i}{y_i} is ordinal with a total of \eqn{C}{C}
categories, with distribution
diff --git a/man/MCMCordfactanal.Rd b/man/MCMCordfactanal.Rd
index a4f3b9f..f8f7272 100644
--- a/man/MCMCordfactanal.Rd
+++ b/man/MCMCordfactanal.Rd
@@ -2,11 +2,11 @@
\alias{MCMCordfactanal}
\title{Markov Chain Monte Carlo for Ordinal Data Factor Analysis Model}
\description{
- This function generates a posterior density sample from an ordinal data
- factor analysis model. Normal priors are assumed on the factor
+ This function generates a sample from the posterior distribution of an
+ ordinal data factor analysis model. Normal priors are assumed on the factor
loadings and factor scores while improper uniform priors are assumed
on the cutpoints. The user supplies data and parameters for the prior
- distributions, and a sample from the posterior density is returned as
+ distributions, and a sample from the posterior distribution is returned as
an mcmc object, which can be subsequently analyzed with
functions provided in the coda package.
}
@@ -110,7 +110,7 @@ MCMCordfactanal(x, factors, lambda.constraints=list(),
}
\value{
- An mcmc object that contains the posterior density sample. This
+ An mcmc object that contains the posterior sample. This
object can be summarized by functions provided by the coda package.
}
@@ -164,7 +164,7 @@ MCMCordfactanal(x, factors, lambda.constraints=list(),
The standard two-parameter item response theory model with probit
link is a special case of the model sketched above.
- \code{MCMCordfactanal} simulates from the posterior density using
+ \code{MCMCordfactanal} simulates from the posterior distribution using
a Metropolis-Hastings within Gibbs sampling algorithm. The algorithm
employed is based on work by Cowles (1996). Note that
the first element of \eqn{\phi_i}{phi_i} is a 1. As a result, the
@@ -174,7 +174,7 @@ MCMCordfactanal(x, factors, lambda.constraints=list(),
returned in the mcmc object.
The simulation proper is done in compiled C++ code to maximize
efficiency. Please consult the coda documentation for a comprehensive
- list of functions that can be used to analyze the posterior density
+ list of functions that can be used to analyze the posterior
sample.
As is the case with all measurement models, make sure that you have plenty
diff --git a/man/MCMCpanel.Rd b/man/MCMCpanel.Rd
index 11ad638..ae8331f 100644
--- a/man/MCMCpanel.Rd
+++ b/man/MCMCpanel.Rd
@@ -2,13 +2,13 @@
\alias{MCMCpanel}
\title{Markov Chain Monte Carlo for the General Linear Panel Model}
\description{
- MCMCpanel generates a posterior density sample from a General
+ MCMCpanel generates a sample from the posterior distribution of a General
Linear Panel Model using Algorithm 2 of Chib and Carlin (1999).
This model uses a multivariate Normal prior for the fixed
effects parameters, a Wishart prior on the random effects
precision matrix, and a Gamma prior on the conditional error
precision. The user supplies data and priors, and a sample from
- the posterior density is returned as an mcmc object,
+ the posterior distribution is returned as an mcmc object,
which can be subsequently analyzed with functions provided in
the coda package.
}
@@ -95,17 +95,17 @@ MCMCpanel(obs, Y, X, W, burnin = 1000, mcmc = 10000, thin = 5,
}
\value{
- An mcmc object that contains the posterior density sample. This
+ An mcmc object that contains the posterior sample. This
object can be summarized by functions provided by the coda package.
}
\details{
- \code{MCMCpanel} simulates from the posterior density sample using
+ \code{MCMCpanel} simulates from the posterior distribution sample using
the blocked Gibbs sampler of Chib and Carlin (1999), Algorithm 2.
The simulation proper
is done in compiled C++ code to maximize efficiency. Please consult
the coda documentation for a comprehensive list of functions that can be
- used to analyze the posterior density sample.
+ used to analyze the posterior sample.
The model takes the following form:
\deqn{y_i = X_i \beta + W_i b_i + \varepsilon_i}{y_i = X_i * beta + W_i *
diff --git a/man/MCMCpoisson.Rd b/man/MCMCpoisson.Rd
index 2b7574d..d2ef0aa 100644
--- a/man/MCMCpoisson.Rd
+++ b/man/MCMCpoisson.Rd
@@ -2,10 +2,10 @@
\alias{MCMCpoisson}
\title{Markov Chain Monte Carlo for Poisson Regression}
\description{
- This function generates a posterior density sample
- from a Poisson regression model using a random walk Metropolis
+ This function generates a sample from the posterior distribution
+ of a Poisson regression model using a random walk Metropolis
algorithm. The user supplies data and priors,
- and a sample from the posterior density is returned as an mcmc
+ and a sample from the posterior distribution is returned as an mcmc
object, which can be subsequently analyzed with functions
provided in the coda package.
}
@@ -31,7 +31,7 @@ MCMCpoisson(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
scalar or a \eqn{k}{k}-vector, where \eqn{k}{k} is the length of
\eqn{\beta}{beta}.Make sure that the
acceptance rate is satisfactory (typically between 0.20 and 0.5)
- before using the posterior density sample for inference.}
+ before using the posterior sample for inference.}
\item{verbose}{A switch which determines whether or not the progress of
@@ -76,15 +76,16 @@ MCMCpoisson(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
}
\value{
- An mcmc object that contains the posterior density sample. This
+ An mcmc object that contains the posterior sample. This
object can be summarized by functions provided by the coda package.
}
-\details{\code{MCMCpoisson} simulates from the posterior density of a Poisson
- regression model using a random walk Metropolis algorithm. The simulation
- proper is done in compiled C++ code to maximize efficiency. Please consult
- the coda documentation for a comprehensive list of functions that can be
- used to analyze the posterior density sample.
+\details{\code{MCMCpoisson} simulates from the posterior distribution of
+ a Poisson regression model using a random walk Metropolis
+ algorithm. The simulation proper is done in compiled C++ code to
+ maximize efficiency. Please consult the coda documentation for a
+ comprehensive list of functions that can be used to analyze the
+ posterior sample.
The model takes the following form:
\deqn{y_i \sim \mathcal{P}oisson(\mu_i)}{y_i ~ Poisson(mu_i)}
diff --git a/man/MCMCprobit.Rd b/man/MCMCprobit.Rd
index 4c64805..3aa6b8a 100644
--- a/man/MCMCprobit.Rd
+++ b/man/MCMCprobit.Rd
@@ -2,10 +2,10 @@
\alias{MCMCprobit}
\title{Markov Chain Monte Carlo for Probit Regression}
\description{
- This function generates a posterior density sample
- from a probit regression model using the data augmentation
+ This function generates a sample from the posterior distribution of
+ a probit regression model using the data augmentation
approach of Albert and Chib (1993). The user supplies data and priors,
- and a sample from the posterior density is returned as an mcmc
+ and a sample from the posterior distribution is returned as an mcmc
object, which can be subsequently analyzed with functions
provided in the coda package.
}
@@ -76,16 +76,16 @@ MCMCprobit(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
}
\value{
- An mcmc object that contains the posterior density sample. This
+ An mcmc object that contains the posterior sample. This
object can be summarized by functions provided by the coda package.
}
\details{
-\code{MCMCprobit} simulates from the posterior density of a probit
+\code{MCMCprobit} simulates from the posterior distribution of a probit
regression model using data augmentation. The simulation
proper is done in compiled C++ code to maximize efficiency. Please consult
the coda documentation for a comprehensive list of functions that can be
- used to analyze the posterior density sample.
+ used to analyze the posterior sample.
The model takes the following form:
\deqn{y_i \sim \mathcal{B}ernoulli(\pi_i)}{y_i ~ Bernoulli(pi_i)}
diff --git a/man/MCMCregress.Rd b/man/MCMCregress.Rd
index ba9307a..1ccb195 100644
--- a/man/MCMCregress.Rd
+++ b/man/MCMCregress.Rd
@@ -2,12 +2,12 @@
\alias{MCMCregress}
\title{Markov Chain Monte Carlo for Gaussian Linear Regression}
\description{
- This function generates a posterior density sample
- from a linear regression model with Gaussian errors using
+ This function generates a sample from the posterior distribution
+ of a linear regression model with Gaussian errors using
Gibbs sampling (with a multivariate Gaussian prior on the
beta vector, and an inverse Gamma prior on the conditional
error variance). The user supplies data and priors, and
- a sample from the posterior density is returned as an mcmc
+ a sample from the posterior distribution is returned as an mcmc
object, which can be subsequently analyzed with functions
provided in the coda package.
}
@@ -81,17 +81,17 @@ MCMCregress(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
}
\value{
- An mcmc object that contains the posterior density sample. This
+ An mcmc object that contains the posterior sample. This
object can be summarized by functions provided by the coda package.
}
\details{
- \code{MCMCregress} simulates from the posterior density using
+ \code{MCMCregress} simulates from the posterior distribution using
standard Gibbs sampling (a multivariate Normal draw for the betas, and an
inverse Gamma draw for the conditional error variance). The simulation
proper is done in compiled C++ code to maximize efficiency. Please consult
the coda documentation for a comprehensive list of functions that can be
- used to analyze the posterior density sample.
+ used to analyze the posterior sample.
The model takes the following form:
\deqn{y_i = x_i ' \beta + \varepsilon_{i}}{y_i = x_i'beta + epsilon_i}
diff --git a/man/MCMCtobit.Rd b/man/MCMCtobit.Rd
index 80872fa..907f301 100644
--- a/man/MCMCtobit.Rd
+++ b/man/MCMCtobit.Rd
@@ -2,13 +2,13 @@
\alias{MCMCtobit}
\title{Markov Chain Monte Carlo for Gaussian Linear Regression with a Censored Dependent Variable}
\description{
- This function generates a posterior density sample
- from a linear regression model with Gaussian errors using
+ This function generates a sample from the posterior distribution
+ of a linear regression model with Gaussian errors using
Gibbs sampling (with a multivariate Gaussian prior on the
beta vector, and an inverse Gamma prior on the conditional
error variance). The dependent variable may be censored
from below, from above, or both. The user supplies data
- and priors, and a sample from the posterior density is
+ and priors, and a sample from the posterior distribution is
returned as an mcmc object, which can be subsequently
analyzed with functions provided in the coda package.
}
@@ -87,14 +87,14 @@ MCMCtobit(formula, data = parent.frame(), below = 0, above = Inf,
}
\details{
- \code{MCMCtobit} simulates from the posterior density using standard
+ \code{MCMCtobit} simulates from the posterior distribution using standard
Gibbs sampling (a multivariate Normal draw for the betas, and an
inverse Gamma draw for the conditional error variance). \code{MCMCtobit}
differs from \code{MCMCregress} in that the dependent variable may be
censored from below, from above, or both. The simulation
proper is done in compiled C++ code to maximize efficiency. Please consult
the coda documentation for a comprehensive list of functions that can be
- used to analyze the posterior density sample.
+ used to analyze the posterior sample.
The model takes the following form:
\deqn{y_i = x_i ' \beta + \varepsilon_{i},}{y_i = x_i'beta + epsilon_i,}
@@ -119,7 +119,7 @@ MCMCtobit(formula, data = parent.frame(), below = 0, above = Inf,
}
\value{
- An mcmc object that contains the posterior density sample. This
+ An mcmc object that contains the posterior sample. This
object can be summarized by functions provided by the coda package.
}
diff --git a/src/MCMCdynamicEI.cc b/src/MCMCdynamicEI.cc
index 25e7bd1..c85fac5 100644
--- a/src/MCMCdynamicEI.cc
+++ b/src/MCMCdynamicEI.cc
@@ -68,7 +68,7 @@ static void doubling(double (*logfun)(double[], const double&, const double&,
theta_R[index] = R;
int K = p;
while (K > 0 &&
- (z < logfun(theta_L, r0, r1, c0, mu0, mu1, sigma0, sigma1) |
+ (z < logfun(theta_L, r0, r1, c0, mu0, mu1, sigma0, sigma1) ||
z < logfun(theta_R, r0, r1, c0, mu0, mu1, sigma0, sigma1))){
double V = stream->runif();
if (V < 0.5){
diff --git a/src/MCMChierEI.cc b/src/MCMChierEI.cc
index f77d381..305ecd5 100644
--- a/src/MCMChierEI.cc
+++ b/src/MCMChierEI.cc
@@ -66,7 +66,7 @@ static void doubling(double (*logfun)(double[], const double&, const double&,
theta_R[index] = R;
int K = p;
while (K > 0 &&
- (z < logfun(theta_L, r0, r1, c0, mu0, mu1, sigma0, sigma1) |
+ (z < logfun(theta_L, r0, r1, c0, mu0, mu1, sigma0, sigma1) ||
z < logfun(theta_R, r0, r1, c0, mu0, mu1, sigma0, sigma1))){
double V = stream->runif();
if (V < 0.5){
diff --git a/src/MCMClogituserprior.cc b/src/MCMClogituserprior.cc
new file mode 100644
index 0000000..e5042e6
--- /dev/null
+++ b/src/MCMClogituserprior.cc
@@ -0,0 +1,215 @@
+// MCMClogituserprior.cc samples from the posterior distribution of a
+// logistic regression model with user-written prior coded in R using a
+// random walk Metropolis algorithm
+//
+// Andrew D. Martin
+// Dept. of Political Science
+// Washington University in St. Louis
+// admartin at wustl.edu
+//
+// Kevin M. Quinn
+// Dept. of Government
+// Harvard University
+// kevin_quinn at harvard.edu
+//
+// This software is distributed under the terms of the GNU GENERAL
+// PUBLIC LICENSE Version 2, June 1991. See the package LICENSE
+// file for more information.
+//
+// Copyright (C) 2004 Andrew D. Martin and Kevin M. Quinn
+//
+// KQ 8/17/2005 (based on current version of MCMCmetrop1R.cc)
+//
+//
+
+#include "matrix.h"
+#include "distributions.h"
+#include "stat.h"
+#include "la.h"
+#include "ide.h"
+#include "smath.h"
+#include "MCMCrng.h"
+#include "MCMCfcds.h"
+
+#include <R.h> // needed to use Rprintf()
+#include <R_ext/Utils.h> // needed to allow user interrupts
+
+using namespace SCYTHE;
+using namespace std;
+
+extern "C" {
+
+#include <Rdefines.h>
+#include <Rinternals.h>
+
+
+ // function that evaluatees the user supplied R function
+ static double user_fun_eval(SEXP fun, SEXP theta, SEXP myframe){
+
+ SEXP R_fcall;
+ if(!isFunction(fun)) error("`fun' must be a function");
+ if(!isEnvironment(myframe)) error("myframe must be an environment");
+ PROTECT(R_fcall = lang2(fun, R_NilValue));
+ SETCADR(R_fcall, theta);
+ SEXP funval = eval(R_fcall, myframe);
+ if (!isReal(funval)) error("`fun' must return a double");
+ double fv = REAL(funval)[0];
+ UNPROTECT(1);
+ return fv;
+ }
+
+ static double logit_loglike(const Matrix<double>& Y,
+ const Matrix<double>& X,
+ const Matrix<double>& beta){
+
+ // likelihood
+ const Matrix<double> eta = X * beta;
+ const Matrix<double> p = 1.0/(1.0 + exp(-eta));
+ double loglike = 0.0;
+ for (int i=0; i<Y.rows(); ++i)
+ loglike += Y[i]*::log(p[i]) + (1-Y[i])*::log(1-p[i]);
+
+ return loglike;
+ }
+
+
+ // the function that actually does the sampling and returns a value to R
+ SEXP MCMClogituserprior_cc(SEXP fun, SEXP Y_R, SEXP X_R,
+ SEXP theta, SEXP myframe, SEXP burnin_R,
+ SEXP mcmc_R, SEXP thin_R,
+ SEXP verbose, SEXP lecuyer_R, SEXP seedarray_R,
+ SEXP lecuyerstream_R, SEXP logfun,
+ SEXP propvar_R){
+
+ // put burnin_R, mcmc_R, and thin_R into ints
+ const int burnin = INTEGER(burnin_R)[0];
+ const int mcmc = INTEGER(mcmc_R)[0];
+ const int thin = INTEGER(thin_R)[0];
+
+ // put rng stuff together and initiate stream
+ int seedarray[6];
+ for(int i=0; i<6; ++i) seedarray[i] = INTEGER(seedarray_R)[i];
+ rng *stream = MCMCpack_get_rng(INTEGER(lecuyer_R)[0],
+ seedarray, INTEGER(lecuyerstream_R)[0]);
+
+ // put propvar_R into a Matrix
+ double* propvar_data = REAL(propvar_R);
+ const int propvar_nr = nrows(propvar_R);
+ const int propvar_nc = ncols(propvar_R);
+ Matrix <double> propvar (propvar_nc, propvar_nr, propvar_data);
+ propvar = t(propvar);
+
+ // put Y_R into a Scythe Matrix
+ int* Y_data = INTEGER(Y_R);
+ const int Y_nr = length(Y_R);
+ const int Y_nc = 1;
+ Matrix <int> Y_M (Y_nc, Y_nr, Y_data);
+ Y_M = t(Y_M);
+
+ // put X_R into a Scythe Matrix
+ double* X_data = REAL(X_R);
+ const int X_nr = nrows(X_R);
+ const int X_nc = ncols(X_R);
+ Matrix <double> X_M (X_nc, X_nr, X_data);
+ X_M = t(X_M);
+
+ // define constants
+ const int npar = length(theta);
+ const int tot_iter = burnin + mcmc;
+ const int nsamp = mcmc/thin;
+ const Matrix <double> propc = cholesky(propvar);
+
+ // define matrix to hold the sample
+ Matrix <double> sample (nsamp, npar);
+
+ // put theta into a Scythe Matrix
+ double* theta_data = REAL(theta);
+ const int theta_nr = length(theta);
+ const int theta_nc = 1;
+ Matrix <double> theta_M (theta_nc, theta_nr, theta_data);
+ theta_M = t(theta_M);
+
+ // evaluate userfun at starting value
+ double loglike_val = logit_loglike(Y_M, X_M, theta_M);
+ double userfun_cur = user_fun_eval(fun, theta, myframe);
+ if (INTEGER(logfun)[0]==0) userfun_cur = ::log(userfun_cur);
+ userfun_cur += loglike_val;
+
+
+ // THE METROPOLIS SAMPLING
+ int count = 0;
+ int accepts = 0;
+ for (int iter=0; iter<tot_iter; ++iter){
+
+ // generate candidate value of theta
+ Matrix <double> theta_can_M = theta_M + propc * stream->rnorm(npar,1);
+
+ // put theta_can_M into a SEXP
+ SEXP theta_can;
+ theta_can = PROTECT(allocVector(REALSXP, npar));
+ for (int i=0; i<npar; ++i){
+ REAL(theta_can)[i] = theta_can_M[i];
+ }
+
+ // evaluate user function fun at candidate theta
+ loglike_val = logit_loglike(Y_M, X_M, theta_can_M);
+ double userfun_can = user_fun_eval(fun, theta_can, myframe);
+ if (INTEGER(logfun)[0]==0) userfun_can = ::log(userfun_can);
+ userfun_can += loglike_val;
+
+ const double ratio = ::exp(userfun_can - userfun_cur);
+
+ if (stream->runif() < ratio){
+ theta = theta_can;
+ theta_M = theta_can_M;
+ userfun_cur = userfun_can;
+ ++accepts;
+ }
+
+ // store values in matrices
+ if ((iter%thin)==0 && iter >= burnin){
+ for (int j = 0; j < npar; j++)
+ sample(count, j) = REAL(theta)[j];
+ ++count;
+ }
+
+ if (INTEGER(verbose)[0] > 0 && iter % INTEGER(verbose)[0] == 0) {
+ Rprintf("MCMClogit iteration %i of %i \n", (iter+1), tot_iter);
+ Rprintf("beta = \n");
+ for (int i=0; i<npar; ++i)
+ Rprintf("%10.5f\n", REAL(theta)[i]);
+ Rprintf("function value = %10.5f\n", userfun_cur);
+ Rprintf("Metropolis acceptance rate = %3.5f\n\n",
+ static_cast<double>(accepts) /
+ static_cast<double>(iter+1));
+ }
+
+
+ UNPROTECT(1);
+ R_CheckUserInterrupt(); // allow user interrupts
+ }
+
+ // put the sample into a SEXP and return it
+ SEXP sample_SEXP;
+ sample_SEXP = PROTECT(allocMatrix(REALSXP, nsamp, npar));
+ for (int i=0; i<nsamp; ++i){
+ for (int j=0; j<npar; ++j){
+ REAL(sample_SEXP)[i + nsamp*j] = sample(i,j);
+ }
+ }
+ UNPROTECT(1);
+
+ delete stream; // clean up random number stream
+
+ // print the the acceptance rate to the console in a way that
+ // everyone (even Windows users) can see
+ Rprintf("\n\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n");
+ Rprintf("The Metropolis acceptance rate was %3.5f",
+ static_cast<double>(accepts) / static_cast<double>(tot_iter));
+ Rprintf("\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n");
+
+ // return the sample
+ return sample_SEXP;
+
+ }
+}
diff --git a/src/MCMCmnlslice.cc b/src/MCMCmnlslice.cc
index 3a34cdb..4f1fe5e 100644
--- a/src/MCMCmnlslice.cc
+++ b/src/MCMCmnlslice.cc
@@ -104,7 +104,7 @@ static double mnl_logpost(const Matrix<double>& Y, const Matrix<double>& X,
beta_R[index] = R;
int K = p;
while (K > 0 &&
- (z < logfun(Y, X, beta_L, beta_prior_mean, beta_prior_prec) |
+ (z < logfun(Y, X, beta_L, beta_prior_mean, beta_prior_prec) ||
z < logfun(Y, X, beta_R, beta_prior_mean, beta_prior_prec))){
double V = stream->runif();
if (V < 0.5){
diff --git a/src/la.cc b/src/la.cc
index 72b6c65..998a3f1 100644
--- a/src/la.cc
+++ b/src/la.cc
@@ -23,6 +23,7 @@
#include <cmath>
#include <algorithm>
+#include <numeric>
#include <set>
#ifdef SCYTHE_COMPILE_DIRECT
@@ -198,7 +199,7 @@ namespace SCYTHE {
}
// See how many rows are true
- int N = accumulate(e.begin(), e.end(), (int) 0);
+ int N = std::accumulate(e.begin(), e.end(), (int) 0);
// declare and form output Matrix
Matrix<T> temp(N, A.cols(), false);
@@ -349,17 +350,37 @@ namespace SCYTHE {
Matrix<T>
crossprod (const Matrix<T> &A)
{
- Matrix<T> temp(A.cols(), A.cols(), false);
-
- for (int i = 0; i < A.cols(); ++i) {
- for (int j = i; j < A.cols(); ++j) {
- temp(i,j) = T (0);
- for (int k = 0; k < A.rows(); ++k)
- temp(j,i) = temp(i,j) += A(k,i) * A(k,j);
- }
- }
+ int rows = A.rows();
+ int cols = A.cols();
+ Matrix<T> result(cols, cols, true);
+ T tmp;
+
+ if (rows == 1) {
+ for (int k = 0; k < rows; ++k) {
+ for (int i = 0; i < cols; ++i) {
+ tmp = A[k * cols + i];
+ for (int j = i; j < cols; ++j) {
+ result[j * cols +i] =
+ result[i * cols + j] += tmp * A[k * cols + j];
+ }
+ }
+ }
+ } else {
+ for (int k = 0; k < rows; ++k) {
+ for (int i = 0; i < cols; ++i) {
+ tmp = A[k * cols + i];
+ for (int j = i; j < cols; ++j) {
+ result[i * cols + j] += tmp * A[k * cols + j];
+ }
+ }
+ }
+
+ for (int i = 0; i < cols; ++i)
+ for (int j = i + 1; j < cols; ++j)
+ result[j * cols + i] = result[i * cols + j];
+ }
- return temp;
+ return result;
}
} // end namespace SCYTHE
diff --git a/src/stat.cc b/src/stat.cc
index 147857c..ea9ee9e 100644
--- a/src/stat.cc
+++ b/src/stat.cc
@@ -42,7 +42,7 @@ namespace SCYTHE {
T
sum (const Matrix<T> &A)
{
- return (accumulate(A.begin(), A.end(), (T) 0));
+ return (std::accumulate(A.begin(), A.end(), (T) 0));
}
/* Calculate the sum of each column in a Matrix */
@@ -53,7 +53,7 @@ namespace SCYTHE {
Matrix<T> temp(1, A.cols(), false);
for (int j = 0; j < A.cols(); ++j)
- temp[j] = accumulate(A.vecc(j), A.vecc(j + 1), (T) 0);
+ temp[j] = std::accumulate(A.vecc(j), A.vecc(j + 1), (T) 0);
return temp;
}
@@ -92,7 +92,7 @@ namespace SCYTHE {
T
mean (const Matrix<T> &A)
{
- return (accumulate(A.begin(), A.end(), (T) 0) / A.size());
+ return (std::accumulate(A.begin(), A.end(), (T) 0) / A.size());
}
/* Calculate the mean of each column of a Matrix */
@@ -103,7 +103,7 @@ namespace SCYTHE {
Matrix<T> temp(1, A.cols(), false);
for (int j = 0; j < A.cols(); ++j)
- temp[j] = accumulate(A.vecc(j), A.vecc(j + 1), (T) 0) / A.rows();
+ temp[j] = std::accumulate(A.vecc(j), A.vecc(j + 1), (T) 0) / A.rows();
return temp;
}
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/r-cran-mcmcpack.git
More information about the debian-science-commits
mailing list