[r-cran-mcmcpack] 11/90: Imported Upstream version 0.7-1
Andreas Tille
tille at debian.org
Fri Dec 16 09:07:34 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 c5b5ece50af46d4c1fb9a0c7dfe97b7fe019c0c9
Author: Andreas Tille <tille at debian.org>
Date: Fri Dec 16 08:07:04 2016 +0100
Imported Upstream version 0.7-1
---
DESCRIPTION | 10 +-
HISTORY | 32 ++
INDEX | 15 +
NAMESPACE | 14 +-
R/BayesFactors.R | 278 ++++++++++++
R/MCMCSVDreg.R | 174 ++++++++
R/MCMCdynamicEI.R | 2 +-
R/MCMCfactanal.R | 6 +-
R/MCMChierEI.R | 4 +-
R/MCMCirt1d.R | 35 +-
R/MCMCirtKdRob.R | 413 +++++++++++++++++
R/MCMClogit.R | 70 ++-
R/MCMCmixfactanal.R | 4 +-
R/MCMCmnl.R | 11 +-
R/MCMCoprobit.R | 6 +-
R/MCMCordfactanal.R | 6 +-
R/MCMCpanel.R | 2 +-
R/MCMCpoisson.R | 47 +-
R/MCMCprobit.R | 50 ++-
R/MCMCregress.R | 86 +++-
R/MCMCtobit.R | 6 +-
R/MCmodels.R | 109 +++++
R/automate.R | 73 ++-
R/distn.R | 8 +-
R/hidden.R | 82 +++-
R/procrust.R | 64 +++
R/utility.R | 4 +-
man/BayesFactor.Rd | 74 +++
man/MCMCSVDreg.Rd | 161 +++++++
man/MCMCirt1d.Rd | 20 +-
man/MCMCirtKdRob.Rd | 358 +++++++++++++++
man/MCMClogit.Rd | 8 +-
man/MCMCpoisson.Rd | 9 +-
man/MCMCprobit.Rd | 10 +-
man/MCMCregress.Rd | 20 +-
man/MCMCtobit.Rd | 4 +-
man/MCbinomialbeta.Rd | 60 +++
man/MCmultinomdirichlet.Rd | 56 +++
man/MCnormalnormal.Rd | 62 +++
man/MCpoissongamma.Rd | 60 +++
man/PostProbMod.Rd | 66 +++
man/procrust.Rd | 51 +++
man/xpnd.Rd | 11 +-
src/MCMCSVDreg.cc | 235 ++++++++++
src/MCMCfcds.cc | 6 +-
src/MCMCirt1d.cc | 37 +-
src/MCMCirtKdRob.cc | 1062 ++++++++++++++++++++++++++++++++++++++++++++
src/MCMCregress.cc | 78 +++-
src/Makevars | 1 -
src/Makevars.in | 1 -
50 files changed, 3903 insertions(+), 158 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index d1f2f5b..7b33dbc 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,11 +1,11 @@
Package: MCMCpack
-Version: 0.6-6
-Date: 2005-11-14
+Version: 0.7-1
+Date: 2006-1-28
Title: Markov chain Monte Carlo (MCMC) Package
Author: Andrew D. Martin <admartin at wustl.edu>, and
Kevin M. Quinn <kevin_quinn at harvard.edu>
Maintainer: Andrew D. Martin <admartin at wustl.edu>
-Depends: R (>= 2.0.1), coda (>= 0.9-2), MASS
+Depends: R (>= 2.2.0), coda (>= 0.10-3), MASS
Description: This package contains functions to perform Bayesian
inference using posterior simulation for a number of statistical
models. All simulation is done in compiled C++ written in the Scythe
@@ -15,6 +15,6 @@ Description: This package contains functions to perform Bayesian
additional density functions and pseudo-random number generators
for statistical distributions, a general purpose Metropolis
sampling algorithm, and tools for visualization.
-License: GPL version 2 or newer
+License: GPL version 2
URL: http://mcmcpack.wustl.edu
-Packaged: Mon Nov 14 09:56:45 2005; adm
+Packaged: Sat Jan 28 09:00:50 2006; adm
diff --git a/HISTORY b/HISTORY
index 698f2d6..be2eaa7 100644
--- a/HISTORY
+++ b/HISTORY
@@ -2,6 +2,38 @@
// Changes and Bug Fixes
//
+0.6-6 to 0.7-1
+ * Added robust k-dimensional IRT model MCMCirtKdRob().
+ * Added SVD regression MCMCSVDreg().
+ * Updated auto.Scythe.call() to have any number of non-constants args; this
+ is necessary to return other things from the C++ code besides the
+ posterior density sample (including log-marginal likelihoods,
+ acceptance rates, etc.).
+ * Updated form.mcmc.object() to add additional attributes to an
+ mcmc object to hold other quantities of interest (such as
+ log-marginal likelihoods, acceptance rates, data, etc.).
+ * For linear regression model, added option to compute log-marginal
+ likelihood via the method of Chib (1995) or via the Laplace Approximation.
+ * For logistic regression model, added option to compute log-marginal
+ likelihood via the Laplace Approximation.
+ * For probit regression model, added option to compute log-marginal
+ likelihood via the Laplace Approximation.
+ * For Poisson regression model, added option to compute log-marginal
+ likelihood via the Laplace Approximation.
+ * Added methods to calculate Bayes factors and posterior probabilities
+ of models and handle log-marginal likelihoods.
+ * Added teaching model MCbinomialbeta().
+ * Added teaching model MCpoissongamma().
+ * Added teaching model MCnormalnormal().
+ * Added teaching model MCmultinomdirichlet().
+ * patched xpnd() function to provide more functionality. Thanks to
+ Gregor Gorjan for the patches.
+ * added the function procrustes() that performs a Procrustes
+ transformation of a matrix.
+ * made some minor changes to MCMCirt1d() to help conserve memory when
+ dealing with large datasets (MORE OPTIMIZATION NEEDS TO BE DONE IN
+ TERMS OF BOTH SPEED AND MEMORY)
+
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.
diff --git a/INDEX b/INDEX
index f99695c..6f3e753 100644
--- a/INDEX
+++ b/INDEX
@@ -1,6 +1,9 @@
+BayesFactor Create an object of class BayesFactor from
+ MCMCpack output
Dirichlet The Dirichlet Distribution
InvGamma The Inverse Gamma Distribution
InvWishart The Inverse Wishart Distribution
+MCMCSVDreg Markov Chain Monte Carlo for SVD Regression
MCMCdynamicEI Markov Chain Monte Carlo for Quinn's Dynamic
Ecological Inference Model
MCMCfactanal Markov Chain Monte Carlo for Normal Theory
@@ -11,6 +14,8 @@ MCMCirt1d Markov Chain Monte Carlo for One Dimensional
Item Response Theory Model
MCMCirtKd Markov Chain Monte Carlo for K-Dimensional Item
Response Theory Model
+MCMCirtKdRob Markov Chain Monte Carlo for Robust
+ K-Dimensional Item Response Theory Model
MCMClogit Markov Chain Monte Carlo for Logistic
Regression
MCMCmetrop1R Metropolis Sampling from User-Written R
@@ -31,16 +36,26 @@ MCMCregress Markov Chain Monte Carlo for Gaussian Linear
Regression
MCMCtobit Markov Chain Monte Carlo for Gaussian Linear
Regression with a Censored Dependent Variable
+MCbinomialbeta Monte Carlo Simulation from a Binomial
+ Likelihood with a Beta Prior
+MCmultinomdirichlet Monte Carlo Simulation from a Multinomial
+ Likelihood with a Dirichlet Prior
+MCnormalnormal Monte Carlo Simulation from a Normal Likelihood
+ (with known variance) with a Normal Prior
+MCpoissongamma Monte Carlo Simulation from a Poisson
+ Likelihood with a Gamma Prior
Nethvote Dutch Voting Behavior in 1989
NoncenHypergeom The Noncentral Hypergeometric Distribution
PErisk Political Economic Risk Data from 62 Countries
in 1987
+PostProbMod Calculate Posterior Probability of Model
Senate 106th U.S. Senate Roll Call Vote Matrix
SupremeCourt U.S. Supreme Court Vote Matrix
Wishart The Wishart Distribution
choicevar Handle Choice-Specific Covariates in
Multinomial Choice Models
dtomogplot Dynamic Tomography Plot
+procrustes Procrustes Transformation
read.Scythe Read a Matrix from a File written by Scythe
tomogplot Tomography Plot
vech Extract Lower Triangular Elements from a
diff --git a/NAMESPACE b/NAMESPACE
index ac1ab79..33acf6c 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -2,7 +2,8 @@ useDynLib(MCMCpack)
import(coda)
import(MASS)
-export(
+export(
+ BayesFactor,
choicevar,
ddirichlet,
dinvgamma,
@@ -10,11 +11,16 @@ export(
dnoncenhypergeom,
dtomogplot,
dwish,
+ MCbinomialbeta,
+ MCpoissongamma,
+ MCnormalnormal,
+ MCmultinomdirichlet,
MCMCdynamicEI,
MCMCfactanal,
MCMChierEI,
MCMCirt1d,
MCMCirtKd,
+ MCMCirtKdRob,
MCMClogit,
MCMCmetrop1R,
MCMCmixfactanal,
@@ -25,7 +31,9 @@ export(
MCMCpoisson,
MCMCprobit,
MCMCregress,
+ MCMCSVDreg,
MCMCtobit,
+ PostProbMod,
rdirichlet,
read.Scythe,
rinvgamma,
@@ -37,3 +45,7 @@ export(
write.Scythe,
xpnd
)
+
+
+S3method(print, BayesFactor)
+S3method(summary, BayesFactor)
diff --git a/R/BayesFactors.R b/R/BayesFactors.R
new file mode 100644
index 0000000..7dc5b20
--- /dev/null
+++ b/R/BayesFactors.R
@@ -0,0 +1,278 @@
+## BayesFactor.R contains functions useful for calculating and comparing
+## marginal likelihoods
+##
+## Originally written by KQ 1/26/2006
+##
+
+## log densities
+"logdinvgamma" <- function(sigma2, a, b){
+ logf <- a * log(b) - lgamma(a) + -(a+1) * log(sigma2) + -b/sigma2
+ return(logf)
+}
+
+"logdmvnorm" <- function(theta, mu, Sigma){
+ d <- length(theta)
+ logf <- -0.5*d * log(2*pi) - 0.5*log(det(Sigma)) -
+ 0.5 * t(theta - mu) %*% solve(Sigma) %*% (theta - mu)
+ return(logf)
+}
+
+
+
+
+
+
+## log posterior densities
+"logpost.regress" <- function(theta, y, X, b0, B0, c0, d0){
+ n <- length(y)
+ k <- ncol(X)
+ beta <- theta[1:k]
+ sigma2 <- exp(theta[k+1])
+
+ Sigma <- solve(B0)
+
+ loglike <- sum(dnorm(y, X%*%beta, sqrt(sigma2), log=TRUE))
+
+ ## note the change to the prior for sigma2 b/c of the transformation
+ logprior <- logdinvgamma(sigma2, c0/2, d0/2) + theta[k+1] +
+ logdmvnorm(beta, b0, Sigma)
+
+ return (loglike + logprior)
+}
+
+
+"logpost.probit" <- function(theta, y, X, b0, B0){
+ n <- length(y)
+ k <- ncol(X)
+ beta <- theta
+
+ p <- pnorm(X %*% beta)
+
+ Sigma <- solve(B0)
+
+ loglike <- sum( y * log(p) + (1-y)*log(1-p) )
+ logprior <- logdmvnorm(beta, b0, Sigma)
+
+ return (loglike + logprior)
+}
+
+"logpost.logit" <- function(theta, y, X, b0, B0){
+ n <- length(y)
+ k <- ncol(X)
+ beta <- theta
+
+ eta <- X %*% beta
+ p <- 1.0/(1.0+exp(-eta))
+
+ Sigma <- solve(B0)
+
+ loglike <- sum( y * log(p) + (1-y)*log(1-p) )
+ logprior <- logdmvnorm(beta, b0, Sigma)
+
+ return (loglike + logprior)
+}
+
+"logpost.logit.userprior" <- function(theta, y, X, userfun, logfun,
+ my.env){
+ n <- length(y)
+ k <- ncol(X)
+ beta <- theta
+
+ eta <- X %*% beta
+ p <- 1.0/(1.0+exp(-eta))
+
+ loglike <- sum( y * log(p) + (1-y)*log(1-p) )
+ if (logfun){
+ logprior <- eval(userfun(theta), envir=my.env)
+ }
+ else{
+ logprior <- log(eval(userfun(theta), envir=my.env))
+ }
+
+ return (loglike + logprior)
+}
+
+
+"logpost.poisson" <- function(theta, y, X, b0, B0){
+ n <- length(y)
+ k <- ncol(X)
+ beta <- theta
+
+ eta <- X %*% beta
+ lambda <- exp(eta)
+
+ Sigma <- solve(B0)
+
+ loglike <- sum(dpois(y, lambda, log=TRUE))
+ logprior <- logdmvnorm(beta, b0, Sigma)
+
+ return (loglike + logprior)
+}
+
+
+
+
+## functions for working with BayesFactor objects
+"BayesFactor" <- function(...){
+ model.list <- list(...)
+ M <- length(model.list)
+ #model.names <- paste("model", 1:M, sep="")
+ this.call <- match.call()
+ this.call.string <- deparse(this.call)
+ this.call.string <- strsplit(this.call.string, "BayesFactor\\(")
+ this.call.string <- this.call.string[[1]][length(this.call.string[[1]])]
+ this.call.string <- strsplit(this.call.string, ",")
+
+ model.names <- NULL
+ for (i in 1:M){
+ model.names <- c(model.names, this.call.string[[1]][i])
+ }
+ model.names <- gsub(")", "", model.names)
+ model.names <- gsub(" ", "", model.names)
+
+ for (i in 1:M){
+ if (!is.mcmc(model.list[[i]])){
+ stop("argument not of class mcmc\n")
+ }
+ }
+
+ BF.mat <- matrix(NA, M, M)
+ BF.log.mat <- matrix(NA, M, M)
+ rownames(BF.mat) <- colnames(BF.mat) <-
+ rownames(BF.log.mat) <- colnames(BF.log.mat) <- model.names
+ BF.logmarglike <- array(NA, M, dimnames=model.names)
+ BF.call <- NULL
+
+ for (i in 1:M){
+ BF.logmarglike[i] <- attr(model.list[[i]], "logmarglike")
+ BF.call <- c(BF.call, attr(model.list[[i]], "call"))
+ for (j in 1:M){
+ if (identical(attr(model.list[[i]], "y"), attr(model.list[[j]], "y"))){
+ BF.log.mat[i,j] <- attr(model.list[[i]], "logmarglike") -
+ attr(model.list[[j]], "logmarglike")
+ BF.mat[i,j] <- exp(BF.log.mat[i,j])
+ }
+ }
+ }
+
+ return(structure(list(BF.mat=BF.mat, BF.log.mat=BF.log.mat,
+ BF.logmarglike=BF.logmarglike,
+ BF.call=BF.call),
+ class="BayesFactor"))
+}
+
+
+"is.BayesFactor" <- function(BF){
+ return(class(BF) == "BayesFactor")
+}
+
+
+"print.BayesFactor" <- function(x, ...){
+
+ cat("The matrix of Bayes Factors is:\n")
+ print(x$BF.mat, digits=3)
+
+ cat("\nThe matrix of the natural log Bayes Factors is:\n")
+ print(x$BF.log.mat, digits=3)
+
+ M <- length(x$BF.call)
+ for (i in 1:M){
+ cat("\n", rownames(x$BF.mat)[i], ":\n")
+ cat(" call = \n")
+ print(x$BF.call[[i]])
+ cat("\n log marginal likelihood = ", x$BF.logmarglike[i], "\n\n")
+
+ }
+
+}
+
+
+"summary.BayesFactor" <- function(object, ...){
+
+ cat("The matrix of Bayes Factors is:\n")
+ print(object$BF.mat, digits=3)
+
+ cat("\nThe matrix of the natural log Bayes Factors is:\n")
+ print(object$BF.log.mat, digits=3)
+
+ BF.mat.NA <- object$BF.mat
+ diag(BF.mat.NA) <- NA
+ minvec <- apply(BF.mat.NA, 1, min, na.rm=TRUE)
+ best.model <- which.max(minvec)
+ if (minvec[best.model] > 150){
+ cat("\nThere is very strong evidence to support",
+ rownames(object$BF.mat)[best.model],
+ "over\nall other models considered.\n")
+ }
+ else if(minvec[best.model] > 20){
+ cat("\nThere is strong evidence or better to support",
+ rownames(object$BF.mat)[best.model],
+ "over\nall other models considered.\n")
+ }
+ else if(minvec[best.model] > 3){
+ cat("\nThere is positive evidence or better to support",
+ rownames(object$BF.mat)[best.model],
+ "over\nall other models considered.\n")
+ }
+ else {
+ cat("\nThe evidence to support",
+ rownames(object$BF.mat)[best.model],
+ "over all\nother models considered is worth no more\n than a bare mention.\n")
+ }
+
+ cat("\n\nStrength of Evidence Guidelines\n(from Kass and Raftery, 1995, JASA)\n")
+ cat("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n")
+ cat("2log(BF[i,j]) BF[i,j] Evidence Against Model j\n")
+ cat("------------------------------------------------------------\n")
+ cat(" 0 to 2 1 to 3 Not worth more than a\n")
+ cat(" bare mention\n")
+ cat(" 2 to 6 3 to 20 Positive\n")
+ cat(" 6 to 10 20 to 150 Strong\n")
+ cat(" >10 >150 Very Strong\n")
+ cat("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n\n")
+
+ M <- length(object$BF.call)
+ for (i in 1:M){
+ cat("\n", rownames(object$BF.mat)[i], ":\n")
+ cat(" call = \n")
+ print(object$BF.call[[i]])
+ cat("\n log marginal likelihood = ", object$BF.logmarglike[i], "\n\n")
+
+ }
+
+}
+
+
+"PostProbMod" <- function(BF, prior.probs=1){
+ if (!is.BayesFactor(BF)){
+ stop("BF is not of class BayesFactor\n")
+ }
+
+ M <- length(BF$BF.call)
+
+ if (min(prior.probs) <= 0){
+ stop("An element of prior.probs is non-positive\n")
+ }
+
+ prior.probs <- rep(prior.probs, M)[1:M]
+ prior.probs <- prior.probs / sum(prior.probs)
+
+ lognumer <- BF$BF.logmarglike + log(prior.probs)
+ maxlognumer <- max(lognumer)
+
+ logpostprobs <- array(NA, M)
+ denom <- 0
+ for (i in 1:M){
+ denom <- denom + exp(lognumer[i]-maxlognumer)
+ }
+ logdenom <- log(denom)
+
+ for (i in 1:M){
+ logpostprobs[i] <- (lognumer[i] - maxlognumer) - logdenom
+ }
+ postprobs <- exp(logpostprobs)
+
+ names(postprobs) <- rownames(BF$BF.mat)
+
+ return(postprobs)
+}
diff --git a/R/MCMCSVDreg.R b/R/MCMCSVDreg.R
new file mode 100644
index 0000000..61bfab2
--- /dev/null
+++ b/R/MCMCSVDreg.R
@@ -0,0 +1,174 @@
+## MCMCSVDreg.R samples from the posterior distribution of a Gaussian
+## linear regression model in which the X matrix has been decomposed
+## with an SVD. Useful for prediction when number of columns of X
+## is (possibly much) greater than the number of rows of X.
+##
+## See West, Mike. 2000. "Bayesian Regression in the 'Large p, Small n'
+## Paradigm." Duke ISDS Discussion Paper 2000-22.
+##
+## KQ 9/9/2005
+##
+
+parse.formula.SVDreg <- function(formula, data, intercept){
+
+ ## extract Y, X, and variable names for model formula and frame
+ mt <- terms(formula, data=data)
+ if(missing(data)) data <- sys.frame(sys.parent())
+ mf <- match.call(expand.dots = FALSE)
+ mf$intercept <- NULL
+ mf$drop.unused.levels <- TRUE
+ mf$na.action <- na.pass ## for specialty handling of missing data
+
+ mf[[1]] <- as.name("model.frame")
+
+ mf <- eval(mf, sys.frame(sys.parent()))
+ if (!intercept){
+ attributes(mt)$intercept <- 0
+ }
+
+ ## null model support
+ X <- if (!is.empty.model(mt)) model.matrix(mt, mf, contrasts)
+ X <- as.matrix(X) # X matrix
+ Y <- as.matrix(model.response(mf, "numeric")) # Y matrix
+
+
+ ## delete obs that are missing in X but potentially keep obs that are
+ ## missing Y
+ ## These are used to for the full SVD of X'
+ keep.indic <- apply(is.na(X), 1, sum) == 0
+ Y.full <- as.matrix(Y[keep.indic,])
+ X.full <- X[keep.indic,]
+ xvars.full <- dimnames(X.full)[[2]] # X variable names
+ xobs.full <- dimnames(X.full)[[1]] # X observation names
+
+ return(list(Y.full, X.full, xvars.full, xobs.full))
+
+}
+
+
+"MCMCSVDreg" <-
+ function(formula, data=parent.frame(), burnin = 1000, mcmc = 10000,
+ thin=1, verbose = 0, seed = NA, tau2.start = 1,
+ g0 = 0, a0 = 0.001, b0 = 0.001, c0=2, d0=2, w0=1,
+ beta.samp=FALSE, intercept=TRUE, ...) {
+
+ # checks
+ check.offset(list(...))
+ check.mcmc.parameters(burnin, mcmc, thin)
+
+ # seeds
+ seeds <- form.seeds(seed)
+ lecuyer <- seeds[[1]]
+ seed.array <- seeds[[2]]
+ lecuyer.stream <- seeds[[3]]
+
+ ## form response and model matrices
+ holder <- parse.formula.SVDreg(formula, data, intercept)
+ Y <- holder[[1]]
+ X <- holder[[2]]
+ xnames <- holder[[3]]
+ obsnames <- holder[[4]]
+ K <- ncol(X) # number of covariates in unidentified model
+ N <- nrow(X) # number of obs (including obs for which predictions
+ # are required) N is also the length of gamma
+ Y.miss.indic <- as.numeric(is.na(Y))
+ n.miss.Y <- sum(Y.miss.indic)
+ Y[is.na(Y)] <- mean(Y, na.rm=TRUE)
+
+
+ ## create SVD representation of t(X)
+ svd.out <- svd(t(X)) # t(X) = A %*% D %*% F
+ A <- svd.out$u
+ D <- diag(svd.out$d)
+ F <- t(svd.out$v)
+
+ ## starting values and priors
+ if (length(tau2.start) < N){
+ tau2.start <- rep(tau2.start, length.out=N)
+ }
+ tau2.start <- matrix(tau2.start, N, 1)
+ mvn.prior <- form.mvn.prior(g0, 0, N)
+ g0 <- mvn.prior[[1]]
+ c0 <- rep(c0, length.out=N)
+ d0 <- rep(d0, length.out=N)
+ check.ig.prior(a0, b0)
+ for (i in 1:N){
+ check.ig.prior(c0[i], d0[i])
+ }
+ w0 <- rep(w0, length.out=N)
+ if (min(w0) < 0 | max(w0) > 1){
+ cat("Element(s) of w0 not in [0, 1].\n")
+ stop("Please respecify and call MCMCSVDreg again.\n")
+ }
+
+ ## define holder for posterior sample
+ if (beta.samp){
+ sample <- matrix(data=0, mcmc/thin, n.miss.Y + 2*N + 1 + K)
+ }
+ else{
+ sample <- matrix(data=0, mcmc/thin, n.miss.Y + 2*N + 1)
+ }
+
+
+ ## call C++ code to draw sample
+ posterior <- .C("MCMCSVDreg",
+ sampledata = as.double(sample),
+ samplerow = as.integer(nrow(sample)),
+ samplecol = as.integer(ncol(sample)),
+ Y = as.double(Y),
+ Yrow = as.integer(nrow(Y)),
+ Ycol = as.integer(ncol(Y)),
+ Ymiss = as.integer(Y.miss.indic),
+ A = as.double(A),
+ Arow = as.integer(nrow(A)),
+ Acol = as.integer(ncol(A)),
+ D = as.double(D),
+ Drow = as.integer(nrow(D)),
+ Dcol = as.integer(ncol(D)),
+ F = as.double(F),
+ Frow = as.integer(nrow(F)),
+ Fcol = as.integer(ncol(F)),
+ burnin = as.integer(burnin),
+ mcmc = as.integer(mcmc),
+ thin = as.integer(thin),
+ lecuyer = as.integer(lecuyer),
+ seedarray = as.integer(seed.array),
+ lecuyerstream = as.integer(lecuyer.stream),
+ verbose = as.integer(verbose),
+ tau2.start = as.double(tau2.start),
+ tau2row = as.integer(nrow(tau2.start)),
+ tau2col = as.integer(ncol(tau2.start)),
+ g0 = as.double(g0),
+ g0row = as.integer(nrow(g0)),
+ g0col = as.integer(ncol(g0)),
+ a0 = as.double(a0),
+ b0 = as.double(b0),
+ c0 = as.double(c0),
+ d0 = as.double(d0),
+ w0 = as.double(w0),
+ betasamp = as.integer(beta.samp),
+ PACKAGE="MCMCpack"
+ )
+
+ ## pull together matrix and build MCMC object to return
+ Y.miss.names <- NULL
+ if (sum(Y.miss.indic) > 0){
+ Y.miss.names <- paste("y", obsnames[Y.miss.indic==1], sep=".")
+ }
+ gamma.names <- paste("gamma", 1:N, sep=".")
+ tau2.names <- paste("tau^2", 1:N, sep=".")
+ beta.names <- paste("beta", xnames, sep=".")
+ if (beta.samp){
+ output <- form.mcmc.object(posterior,
+ names=c(Y.miss.names, gamma.names, tau2.names,
+ "sigma^2", beta.names),
+ title="MCMCSVDreg Posterior Sample")
+ }
+ else{
+ output <- form.mcmc.object(posterior,
+ names=c(Y.miss.names, gamma.names, tau2.names,
+ "sigma^2"),
+ title="MCMCSVDreg Posterior Sample")
+ }
+ return(output)
+ }
diff --git a/R/MCMCdynamicEI.R b/R/MCMCdynamicEI.R
index ac43d71..6e80575 100644
--- a/R/MCMCdynamicEI.R
+++ b/R/MCMCdynamicEI.R
@@ -119,7 +119,7 @@
p1names <- paste("p1table", 1:ntables, sep="")
varnames(output) <- c(p0names, p1names, "sigma^2_0", "sigma^2_1")
- attr(output, "title") <- "MCMCpack Quinn's Dynamic EI Model Posterior Density Sample"
+ attr(output, "title") <- "MCMCpack Quinn's Dynamic EI Model Posterior Sample"
return(output)
diff --git a/R/MCMCfactanal.R b/R/MCMCfactanal.R
index 670fa78..2214f3b 100644
--- a/R/MCMCfactanal.R
+++ b/R/MCMCfactanal.R
@@ -103,7 +103,7 @@
Psi <- factuniqueness.start(psi.start, X)
- ## define holder for posterior density sample
+ ## define holder for posterior sample
if(store.scores == FALSE) {
sample <- matrix(data=0, mcmc/thin, K*factors+K)
}
@@ -113,7 +113,7 @@
## call C++ code to do the sampling
auto.Scythe.call(output.object="posterior", cc.fun.name="MCMCfactanal",
- sample=sample, X=X, burnin=as.integer(burnin),
+ sample.nonconst=sample, X=X, burnin=as.integer(burnin),
mcmc=as.integer(mcmc), thin=as.integer(thin),
lecuyer=as.integer(lecuyer),
seedarray=as.integer(seed.array),
@@ -136,7 +136,7 @@
rep(1:factors,factors), sep="_")
par.names <- c(par.names, phi.names)
}
- title <- "MCMCpack Factor Analysis Posterior Density Sample"
+ title <- "MCMCpack Factor Analysis Posterior Sample"
output <- form.mcmc.object(posterior, par.names, title)
## get rid of columns for constrained parameters
output.df <- as.data.frame(as.matrix(output))
diff --git a/R/MCMChierEI.R b/R/MCMChierEI.R
index 31d50e1..688d969 100644
--- a/R/MCMChierEI.R
+++ b/R/MCMChierEI.R
@@ -1,4 +1,4 @@
-# sample from the posterior distribution of Wakefield's baseline model
+# sample from the posterior distribution of Wakefield's hierarchical model
# for ecological inference in R using linked C++ code in Scythe
#
# KQ 10/22/2002
@@ -127,7 +127,7 @@
varnames(output) <- c(p0names, p1names, "mu0", "mu1", "sigma^2.0",
"sigma^2.1")
- attr(output, "title") <- "MCMCpack Wakefield's Hierarchical EI Model Posterior Density Sample"
+ attr(output, "title") <- "MCMCpack Wakefield's Hierarchical EI Model Posterior Sample"
return(output)
diff --git a/R/MCMCirt1d.R b/R/MCMCirt1d.R
index 972635e..5efdb22 100644
--- a/R/MCMCirt1d.R
+++ b/R/MCMCirt1d.R
@@ -3,12 +3,13 @@
##
## ADM and KQ 1/23/2003
## updated extensively ADM & KQ 7/28/2004
-
+## store.ability arg added KQ 1/27/2006
+
"MCMCirt1d" <-
function(datamatrix, theta.constraints=list(), burnin = 1000,
mcmc = 20000, thin=1, verbose = 0, seed = NA,
theta.start = NA, alpha.start = NA, beta.start = NA, t0 = 0,
- T0 = 1, ab0=0, AB0=.25, store.item = FALSE,
+ T0 = 1, ab0=0, AB0=.25, store.item = FALSE, store.ability=TRUE,
drop.constant.items=TRUE, ... ) {
## checks
@@ -109,13 +110,21 @@
call.=FALSE)
}
- ## define holder for posterior density sample
- if(store.item == FALSE) {
+ ## define holder for posterior sample
+ if(store.item == FALSE & store.ability == TRUE) {
sample <- matrix(data=0, mcmc/thin, J)
}
- else {
+ else if (store.item == TRUE & store.ability == FALSE){
+ sample <- matrix(data=0, mcmc/thin, 2*K)
+ }
+ else if (store.item == TRUE & store.ability == TRUE){
sample <- matrix(data=0, mcmc/thin, J + 2 * K)
}
+ else{
+ cat("Error: store.item == FALSE & store.ability == FALSE.\n")
+ stop("Please respecify and call MCMCirt1d() again.\n",
+ call.=FALSE)
+ }
## seeds
seeds <- form.seeds(seed)
@@ -161,7 +170,8 @@
thetaineqdata = as.double(theta.ineq.constraints),
thetaineqrow = as.integer(nrow(theta.ineq.constraints)),
thetaineqcol = as.integer(ncol(theta.ineq.constraints)),
- store = as.integer(store.item),
+ storei = as.integer(store.item),
+ storea = as.integer(store.ability),
PACKAGE="MCMCpack"
)
@@ -175,16 +185,17 @@
posterior$samplecol,
byrow=TRUE)
output <- mcmc(data=sample, start=burnin+1, end=burnin+mcmc, thin=thin)
-
- if(store.item == FALSE) {
- names <- theta.names
+
+ names <- NULL
+ if(store.ability == TRUE) {
+ names <- c(names, theta.names)
}
- else {
- names <- c(theta.names, alpha.beta.names)
+ if (store.item == TRUE){
+ names <- c(names, alpha.beta.names)
}
varnames(output) <- names
attr(output,"title") <-
- "MCMCirt1d Posterior Density Sample"
+ "MCMCirt1d Posterior Sample"
return(output)
}
diff --git a/R/MCMCirtKdRob.R b/R/MCMCirtKdRob.R
new file mode 100644
index 0000000..f2fa4f7
--- /dev/null
+++ b/R/MCMCirtKdRob.R
@@ -0,0 +1,413 @@
+##########################################################################
+## sample from a K-dimensional two-parameter item response model with
+## logit link that has been rescaled so that the inverse link is:
+##
+## \delta0 + (1 - \delta0 - \delta1)*\Phi(.)
+##
+## where \delta0 \in (0, k0) and \delta1 \in (0, k1)
+##
+## priors for deltas are rescaled beta with parameters c0, d0, and c1, d1
+##
+##
+## datamatrix is assumed to be nsubjects by nitems
+##
+## Andrew D. Martin
+## Washington University
+##
+## Kevin M. Quinn
+## Harvard University
+##
+## Feb. 17, 2005
+##
+##########################################################################
+
+
+
+"MCMCirtKdRob" <-
+ function(datamatrix, dimensions, item.constraints=list(),
+ ability.constraints=list(),
+ burnin = 500, mcmc = 5000, thin=1, interval.method="step",
+ theta.w=0.5, theta.mp=4, alphabeta.w=1.0, alphabeta.mp=4,
+ delta0.w=NA, delta0.mp=3, delta1.w=NA, delta1.mp=3,
+ verbose = FALSE, seed = NA, theta.start=NA,
+ alphabeta.start = NA, delta0.start = NA,
+ delta1.start = NA, b0=0, B0=0,
+ k0=.1, k1=.1, c0=1, d0=1, c1=1, d1=1,
+ store.item=TRUE, store.ability=FALSE,
+ drop.constant.items=TRUE, ... ) {
+
+ ## set X up
+ if (drop.constant.items==TRUE){
+ x.col.var <- apply(datamatrix, 2, var, na.rm=TRUE)
+ keep.inds <- x.col.var>0
+ keep.inds[is.na(keep.inds)] <- FALSE
+ datamatrix <- datamatrix[,keep.inds]
+ }
+ X <- as.data.frame(datamatrix)
+ xvars <- dimnames(X)[[2]]
+ xobs <- dimnames(X)[[1]]
+ N <- nrow(X) # number of subjects
+ K <- ncol(X) # number of items
+ for (i in 1:K){
+ X[is.na(X[,i]), i] <- -999
+ }
+ if(sum(datamatrix==1 | datamatrix==0 | is.na(datamatrix)) != (N*K)) {
+ cat("Error: Data matrix contains elements other than 0, 1 or NA.\n")
+ stop("Please check data and try MCMCirtKdRob() again.\n",
+ call.=FALSE)
+ }
+ X <- as.matrix(X)
+
+ ## take care of the case where X has no row names
+ if (is.null(xobs)){
+ xobs <- 1:N
+ }
+
+ check.offset(list(...))
+ check.mcmc.parameters(burnin, mcmc, thin)
+
+ ## check slice sampling parameters
+ if (!(interval.method %in% c("step", "doubling"))){
+ cat("Error: interval.method not equal to 'step' or 'doubling'.\n")
+ stop("Please check data and try MCMCirtKdRob() again.\n",
+ call.=FALSE)
+ }
+ method.step <- 0
+ if (interval.method == "step"){
+ method.step <- 1
+ }
+ if (theta.w <= 0 ){
+ cat("Error: theta.w not > 0.\n")
+ stop("Please check data and try MCMCirtKdRob() again.\n",
+ call.=FALSE)
+ }
+ if (theta.mp < 1 ){
+ cat("Error: theta.mp not >= 1.\n")
+ stop("Please check data and try MCMCirtKdRob() again.\n",
+ call.=FALSE)
+ }
+ if (alphabeta.w <= 0 ){
+ cat("Error: alphabeta.w not > 0.\n")
+ stop("Please check data and try MCMCirtKdRob() again.\n",
+ call.=FALSE)
+ }
+ if (alphabeta.mp < 1 ){
+ cat("Error: alphabeta.mp not >= 1.\n")
+ stop("Please check data and try MCMCirtKdRob() again.\n",
+ call.=FALSE)
+ }
+ if (is.na(delta0.w)){
+ delta0.w <- 0.25*k0
+ }
+ if (delta0.w <= 0 ){
+ cat("Error: delta0.w not > 0.\n")
+ stop("Please check data and try MCMCirtKdRob() again.\n",
+ call.=FALSE)
+ }
+ if (delta0.mp < 1 ){
+ cat("Error: delta0.mp not >= 1.\n")
+ stop("Please check data and try MCMCirtKdRob() again.\n",
+ call.=FALSE)
+ }
+ if (is.na(delta1.w)){
+ delta1.w <- 0.25*k1
+ }
+ if (delta1.w <= 0 ){
+ cat("Error: delta1.w not > 0.\n")
+ stop("Please check data and try MCMCirtKdRob() again.\n",
+ call.=FALSE)
+ }
+ if (delta1.mp < 1 ){
+ cat("Error: delta1.mp not >= 1.\n")
+ stop("Please check data and try MCMCirtKdRob() again.\n",
+ call.=FALSE)
+ }
+
+ ## error check the prior parameters for delta
+ if (k0 < 0 | k0 > 0.5){
+ cat("Error: k0 not in (0, 0.5).\n")
+ stop("Please check data and try MCMCirtKdRob() again.\n",
+ call.=FALSE)
+ }
+ if (k1 < 0 | k1 > 0.5){
+ cat("Error: k1 not in (0, 0.5).\n")
+ stop("Please check data and try MCMCirtKdRob() again.\n",
+ call.=FALSE)
+ }
+ if (c0 < 0){
+ cat("Error: c0 < 0.\n")
+ stop("Please check data and try MCMCirtKdRob() again.\n",
+ call.=FALSE)
+ }
+ if (c1 < 0){
+ cat("Error: c1 < 0.\n")
+ stop("Please check data and try MCMCirtKdRob() again.\n",
+ call.=FALSE)
+ }
+ if (d0 < 0){
+ cat("Error: d0 < 0.\n")
+ stop("Please check data and try MCMCirtKdRob() again.\n",
+ call.=FALSE)
+ }
+ if (d1 < 0){
+ cat("Error: d1 < 0.\n")
+ stop("Please check data and try MCMCirtKdRob() again.\n",
+ call.=FALSE)
+ }
+
+ ## setup constraints on Lambda = (alpha, beta)
+ holder <- build.factor.constraints(item.constraints, X, K, dimensions+1)
+ Lambda.eq.constraints <- holder[[1]]
+ Lambda.ineq.constraints <- holder[[2]]
+ X.names <- holder[[3]]
+
+ ## setup constraints on theta
+ holder <- build.factor.constraints(ability.constraints, t(X),
+ N, dimensions)
+ theta.eq.constraints <- holder[[1]]
+ theta.ineq.constraints <- holder[[2]]
+
+ ## setup prior on Lambda
+ holder <- form.factload.norm.prior(b0, B0, K, dimensions+1, X.names)
+ Lambda.prior.mean <- holder[[1]]
+ Lambda.prior.prec <- holder[[2]]
+
+ # seeds
+ seeds <- form.seeds(seed)
+ lecuyer <- seeds[[1]]
+ seed.array <- seeds[[2]]
+ lecuyer.stream <- seeds[[3]]
+
+ ## Starting values for delta0 and delta1
+ if (is.na(delta0.start)){
+ delta0.start <- 0.5 * k0;
+ }
+ if (is.na(delta1.start)){
+ delta1.start <- 0.5 * k1;
+ }
+ if (delta0.start < 0 | delta0.start > k0){
+ cat("Error: delta0 not in (0, k0).\n")
+ stop("Please check data and try MCMCirtKdRob() again.\n",
+ call.=FALSE)
+ }
+ if (delta1.start < 0 | delta1.start > k1){
+ cat("Error: delta1 not in (0, k1).\n")
+ stop("Please check data and try MCMCirtKdRob() again.\n",
+ call.=FALSE)
+ }
+
+
+
+ ## Starting values for Lambda
+ Lambda <- matrix(0, K, dimensions+1)
+ if (is.na(alphabeta.start)){# sets Lambda to equality constraints & 0s
+ for (i in 1:K){
+ for (j in 1:(dimensions+1)){
+ if (Lambda.eq.constraints[i,j]==-999){
+ if(Lambda.ineq.constraints[i,j]==0){
+ if (j==1){
+ probit.out <- glm(as.factor(X[X[,i]!=-999,i])~1,
+ family=binomial(link=logit))
+ probit.beta <- coef(probit.out)
+ Lambda[i,j] <- -1 * probit.beta[1]
+ }
+ }
+ if(Lambda.ineq.constraints[i,j]>0){
+ Lambda[i,j] <- 1.0
+ }
+ if(Lambda.ineq.constraints[i,j]<0){
+ Lambda[i,j] <- -1.0
+ }
+ }
+ else Lambda[i,j] <- Lambda.eq.constraints[i,j]
+ }
+ }
+ }
+ else if (is.matrix(alphabeta.start)){
+ if (nrow(alphabeta.start)==K && ncol(alphabeta.start)==(dimensions+1))
+ Lambda <- alphabeta.start
+ else {
+ cat("Starting values not of correct size for model specification.\n")
+ stop("Please respecify and call MCMCirtKdRob() again\n", call.=FALSE)
+ }
+ }
+ else if (length(alphabeta.start)==1 && is.numeric(alphabeta.start)){
+ Lambda <- matrix(alphabeta.start, K, dimensions+1)
+ for (i in 1:K){
+ for (j in 1:(dimensions+1)){
+ if (Lambda.eq.constraints[i,j] != -999)
+ Lambda[i,j] <- Lambda.eq.constraints[i,j]
+ }
+ }
+ }
+ else {
+ cat("Starting values for alpha & beta neither NA, matrix, nor scalar.\n")
+ stop("Please respecify and call MCMCirtKdRob() again\n", call.=FALSE)
+ }
+ for (i in 1:K){
+ lam.sqdist <- sum(Lambda[i,]^2)
+ while (lam.sqdist > 100){
+ Lambda[i,] <- Lambda[i,] * 0.95
+ lam.sqdist <- sum(Lambda[i,]^2)
+ }
+ }
+
+
+
+ ## Starting values for theta
+ if (is.na(theta.start)){
+ theta <- matrix(0, N, dimensions)
+ }
+ else if(is.null(dim(theta.start)) & is.numeric(theta.start)){
+ theta <- matrix(theta.start, N, dimensions)
+ }
+ else if(nrow(theta.start)==N & ncol(theta.start)==dimensions){
+ theta <- theta.start
+ }
+ else{
+ cat("Starting values for theta neither NA, matrix, nor scalar.\n")
+ stop("Please respecify and call MCMCirtKdRob() again\n", call.=FALSE)
+ }
+ for (i in 1:N){
+ for (j in 1:dimensions){
+ if (theta.eq.constraints[i,j]==-999){
+ if(theta.ineq.constraints[i,j]>0){
+ theta[i,j] <- 0.5
+ }
+ if(theta.ineq.constraints[i,j]<0){
+ theta[i,j] <- -0.5
+ }
+ }
+ else theta[i,j] <- theta.eq.constraints[i,j]
+ }
+ }
+
+
+
+
+
+ ## define holder for posterior sample
+ if (store.ability == FALSE && store.item == FALSE){
+ cat("You need to store either the ability or item parameters.\n")
+ stop("Please respecify and call MCMCirtKdRob() again\n", call.=FALSE)
+ }
+ else if (store.ability == TRUE && store.item == FALSE){
+ sample <- matrix(data=0, mcmc/thin, (dimensions+1)*N+2)
+ }
+ else if(store.ability == FALSE && store.item == TRUE) {
+ sample <- matrix(data=0, mcmc/thin, K*(dimensions+1)+2)
+ }
+ else { # store.ability==TRUE && store.item==TRUE
+ sample <- matrix(data=0, mcmc/thin, K*(dimensions+1)+(dimensions+1)*N+2)
+ }
+
+
+ ## Call the C++ code to do the real work
+ posterior <- .C("irtKdRobpost",
+ samdata = as.double(sample),
+ samrow = as.integer(nrow(sample)),
+ samcol = as.integer(ncol(sample)),
+ X = as.integer(X),
+ Xrow = as.integer(nrow(X)),
+ Xcol = as.integer(ncol(X)),
+ burnin = as.integer(burnin),
+ mcmc = as.integer(mcmc),
+ thin = as.integer(thin),
+ lecuyer = as.integer(lecuyer),
+ seedarray = as.integer(seed.array),
+ lecuyerstream = as.integer(lecuyer.stream),
+ verbose = as.integer(verbose),
+ method.step = as.integer(method.step),
+ theta.w = as.double(theta.w),
+ theta.mp = as.integer(theta.mp),
+ ab.w = as.double(alphabeta.w),
+ ab.mp = as.integer(alphabeta.mp),
+ delta0.w = as.double(delta0.w),
+ delta0.mp = as.integer(delta0.mp),
+ delta1.w = as.double(delta1.w),
+ delta1.mp = as.integer(delta1.mp),
+ delta0 = as.double(delta0.start),
+ delta1 = as.double(delta1.start),
+ Lambda = as.double(Lambda),
+ Lambdarow = as.integer(nrow(Lambda)),
+ Lambdacol = as.integer(ncol(Lambda)),
+ theta = as.double(theta),
+ thetarow = as.integer(nrow(theta)),
+ thetacol = as.integer(ncol(theta)),
+ Lameq = as.double(Lambda.eq.constraints),
+ Lameqrow = as.integer(nrow(Lambda.eq.constraints)),
+ Lameqcol = as.integer(ncol(Lambda.ineq.constraints)),
+ Lamineq = as.double(Lambda.ineq.constraints),
+ Lamineqrow = as.integer(nrow(Lambda.ineq.constraints)),
+ Lamineqcol = as.integer(ncol(Lambda.ineq.constraints)),
+ theteq = as.double(theta.eq.constraints),
+ theteqrow = as.integer(nrow(theta.eq.constraints)),
+ theteqcol = as.integer(ncol(theta.ineq.constraints)),
+ thetineq = as.double(theta.ineq.constraints),
+ thetineqrow = as.integer(nrow(theta.ineq.constraints)),
+ thetineqcol = as.integer(ncol(theta.ineq.constraints)),
+ Lampmean = as.double(Lambda.prior.mean),
+ Lampmeanrow = as.integer(nrow(Lambda.prior.mean)),
+ Lampmeancol = as.integer(ncol(Lambda.prior.prec)),
+ Lampprec = as.double(Lambda.prior.prec),
+ Lampprecrow = as.integer(nrow(Lambda.prior.prec)),
+ Lamppreccol = as.integer(ncol(Lambda.prior.prec)),
+ k0 = as.double(k0),
+ k1 = as.double(k1),
+ c0 = as.double(c0),
+ c1 = as.double(c1),
+ d0 = as.double(d0),
+ d1 = as.double(d1),
+ storeitem = as.integer(store.item),
+ storesability = as.integer(store.ability),
+ PACKAGE="MCMCpack"
+ )
+
+
+
+ ## put together matrix and build MCMC object to return
+ sample <- matrix(posterior$samdata, posterior$samrow, posterior$samcol,
+ byrow=TRUE)
+ output <- mcmc(data=sample,start=1, end=mcmc, thin=thin)
+
+ par.names <- NULL
+ if (store.item==TRUE){
+ alpha.hold <- paste("alpha", X.names, sep=".")
+ beta.hold <- paste("beta", X.names, sep = ".")
+ beta.hold <- rep(beta.hold, dimensions, each=dimensions)
+ beta.hold <- paste(beta.hold, 1:dimensions, sep=".")
+
+ Lambda.names <- t(cbind(matrix(alpha.hold, K, 1),
+ matrix(beta.hold,K,dimensions,byrow=TRUE)))
+ dim(Lambda.names) <- NULL
+ par.names <- c(par.names, Lambda.names)
+ }
+
+ if (store.ability==TRUE){
+ phi.names <- paste(paste("theta",
+ rep(xobs, each=(dimensions+1)), sep="."),
+ rep(0:dimensions,(dimensions+1)), sep=".")
+ par.names <- c(par.names, phi.names)
+ }
+
+ par.names <- c("delta0", "delta1", par.names)
+
+ varnames(output) <- par.names
+
+ ## get rid of columns for constrained parameters
+ output.df <- as.data.frame(as.matrix(output))
+ output.var <- diag(var(output.df))
+ output.df <- output.df[,output.var != 0]
+ output <- mcmc(as.matrix(output.df), start=1, end=mcmc, thin=thin)
+
+ ## add constraint info so this isn't lost
+ attr(output, "constraints") <- item.constraints
+ attr(output, "n.items") <- K
+ attr(output, "n.dimensions") <- dimensions
+ attr(output,"title") <-
+ "MCMCpack Robust K-Dimensional Item Response Theory Model Posterior Sample"
+
+return(output)
+
+
+ }
diff --git a/R/MCMClogit.R b/R/MCMClogit.R
index ca2fca2..a30cb4c 100644
--- a/R/MCMClogit.R
+++ b/R/MCMClogit.R
@@ -7,16 +7,31 @@
## 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
+## Modified to handle marginal likelihood calculation 1/27/2006 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, user.prior.density=NULL, logfun=TRUE, ...) {
+ b0 = 0, B0 = 0, user.prior.density=NULL, logfun=TRUE,
+ marginal.likelihood = c("none", "Laplace"), ...) {
## checks
check.offset(list(...))
check.mcmc.parameters(burnin, mcmc, thin)
+
+ cl <- match.call()
+
+ ## get marginal likelihood argument
+ marginal.likelihood <- match.arg(marginal.likelihood)
+ if (max(B0==0) == 1 & is.null(user.prior.density)){
+ if (marginal.likelihood != "none"){
+ warning("Cannot calculate marginal likelihood with improper prior\n")
+ marginal.likelihood <- "none"
+ }
+ }
+ logmarglike <- NULL
+
## seeds
seeds <- form.seeds(seed)
@@ -102,7 +117,8 @@
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),
+ sample.nonconst=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),
@@ -110,10 +126,29 @@
verbose=as.integer(verbose), betastart=beta.start,
b0=b0, B0=B0, V=V)
+ ## marginal likelihood calculation if Laplace
+ if (marginal.likelihood == "Laplace"){
+ theta.start <- beta.start
+ optim.out <- optim(theta.start, logpost.logit, method="BFGS",
+ control=list(fnscale=-1),
+ hessian=TRUE, y=Y, X=X, b0=b0, B0=B0)
+
+ theta.tilde <- optim.out$par
+ beta.tilde <- theta.tilde[1:K]
+
+ Sigma.tilde <- solve(-1*optim.out$hessian)
+
+ logmarglike <- (length(theta.tilde)/2)*log(2*pi) +
+ log(sqrt(det(Sigma.tilde))) +
+ logpost.logit(theta.tilde, Y, X, b0, B0)
+
+ }
+
+
## put together matrix and build MCMC object to return
output <- form.mcmc.object(posterior, names=xnames,
- title="MCMClogit Posterior Density Sample")
-
+ title="MCMClogit Posterior Sample",
+ y=Y, call=cl, logmarglike=logmarglike)
}
else {
@@ -129,11 +164,36 @@
as.logical(logfun),
as.matrix(propvar),
PACKAGE="MCMCpack")
+
+
+ ## marginal likelihood calculation if Laplace
+ if (marginal.likelihood == "Laplace"){
+ theta.start <- beta.start
+ optim.out <- optim(theta.start, logpost.logit.userprior, method="BFGS",
+ control=list(fnscale=-1),
+ hessian=TRUE, y=Y, X=X, userfun=userfun,
+ logfun=logfun, my.env=my.env)
+
+ theta.tilde <- optim.out$par
+ beta.tilde <- theta.tilde[1:K]
+
+ Sigma.tilde <- solve(-1*optim.out$hessian)
+
+ logmarglike <- (length(theta.tilde)/2)*log(2*pi) +
+ log(sqrt(det(Sigma.tilde))) +
+ logpost.logit.userprior(theta.tilde, Y, X, userfun=userfun,
+ logfun=logfun, my.env=my.env)
+
+ }
+
output <- mcmc(data=sample, start=burnin+1,
end=burnin+mcmc, thin=thin)
varnames(output) <- as.list(xnames)
- attr(output, "title") <- "MCMClogit Posterior Density Sample"
+ attr(output, "title") <- "MCMClogit Posterior Sample"
+ attr(output, "y") <- Y
+ attr(output, "call") <- cl
+ attr(output, "logmarglike") <- logmarglike
}
return(output)
diff --git a/R/MCMCmixfactanal.R b/R/MCMCmixfactanal.R
index 14b017f..0b9a3e6 100644
--- a/R/MCMCmixfactanal.R
+++ b/R/MCMCmixfactanal.R
@@ -243,7 +243,7 @@
}
}
- # define holder for posterior density sample
+ # define holder for posterior sample
if (store.scores == FALSE && store.lambda == FALSE){
sample <- matrix(data=0, mcmc/thin, length(gamma)+K)
}
@@ -367,7 +367,7 @@
attr(output, "n.factors") <- factors
attr(output, "accept.rates") <- t(accepts) / (posterior$burnin+posterior$mcmc)
attr(output,"title") <-
- "MCMCpack Mixed Data Factor Analysis Posterior Density Sample"
+ "MCMCpack Mixed Data Factor Analysis Posterior Sample"
return(output)
diff --git a/R/MCMCmnl.R b/R/MCMCmnl.R
index 3827222..19f476a 100644
--- a/R/MCMCmnl.R
+++ b/R/MCMCmnl.R
@@ -1,4 +1,4 @@
-## MCMCmnl.R samples from the posterior density of a multinomial
+## MCMCmnl.R samples from the posterior distribution of a multinomial
## logit model using Metropolis-Hastings.
##
## KQ 12/22/2004
@@ -247,13 +247,14 @@
b0 <- mvn.prior[[1]]
B0 <- mvn.prior[[2]]
- ## define holder for posterior density sample
+ ## define holder for posterior sample
sample <- matrix(data=0, mcmc/thin, dim(X)[2] )
if (mcmc.method=="MH"){
## call C++ code to draw sample
auto.Scythe.call(output.object="posterior", cc.fun.name="MCMCmnlMH",
- sample=sample, Y=Y, X=X, burnin=as.integer(burnin),
+ sample.nonconst=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),
@@ -265,7 +266,7 @@
## put together matrix and build MCMC object to return
output <- form.mcmc.object(posterior, names=xnames,
- title="MCMCmnl Posterior Density Sample")
+ title="MCMCmnl Posterior Sample")
}
else if (mcmc.method=="slice"){
## call C++ code to draw sample
@@ -280,7 +281,7 @@
## put together matrix and build MCMC object to return
output <- form.mcmc.object(posterior, names=xnames,
- title="MCMCmnl Posterior Density Sample")
+ title="MCMCmnl Posterior Sample")
}
diff --git a/R/MCMCoprobit.R b/R/MCMCoprobit.R
index f072c47..36633dd 100644
--- a/R/MCMCoprobit.R
+++ b/R/MCMCoprobit.R
@@ -110,12 +110,12 @@
gamma[3:ncat] <- (polr.out$zeta[2:(ncat-1)] - polr.out$zeta[1])*.588
gamma[ncat+1] <- 300
- ## posterior density sample
+ ## posterior sample
sample <- matrix(data=0, mcmc/thin, K + ncat + 1)
## call C++ code to draw sample
auto.Scythe.call(output.object="posterior", cc.fun.name="MCMCoprobit",
- sample=sample, Y=as.integer(Y), X=X,
+ sample.nonconst=sample, Y=as.integer(Y), X=X,
burnin=as.integer(burnin),
mcmc=as.integer(mcmc), thin=as.integer(thin),
tune=as.double(tune), lecuyer=as.integer(lecuyer),
@@ -131,7 +131,7 @@
output <- mcmc(data=sample, start=burnin+1, end=burnin+mcmc, thin=thin)
xnames <- c(X.names, paste("gamma", 2:(ncat-1), sep=""))
varnames(output) <- xnames
- attr(output, "title") <- "MCMCoprobit Posterior Density Sample"
+ attr(output, "title") <- "MCMCoprobit Posterior Sample"
return(output)
}
diff --git a/R/MCMCordfactanal.R b/R/MCMCordfactanal.R
index 11a7daf..31a0d90 100644
--- a/R/MCMCordfactanal.R
+++ b/R/MCMCordfactanal.R
@@ -225,7 +225,7 @@
}
}
- ## define holder for posterior density sample
+ ## define holder for posterior sample
if (store.scores == FALSE && store.lambda == FALSE){
sample <- matrix(data=0, mcmc/thin, length(gamma))
}
@@ -362,11 +362,11 @@
attr(output, "accept.rates") <- t(accepts) / (posterior$burnin+posterior$mcmc)
if(case.switch==1) {
attr(output,"title") <-
- "MCMCpack Ordinal Data Factor Analysis Posterior Density Sample"
+ "MCMCpack Ordinal Data Factor Analysis Posterior Sample"
}
if(case.switch==2) {
attr(output,"title") <-
- "MCMCpack K-Dimensional Item Response Theory Model Posterior Density Sample"
+ "MCMCpack K-Dimensional Item Response Theory Model Posterior Sample"
}
return(output)
diff --git a/R/MCMCpanel.R b/R/MCMCpanel.R
index 5b32ada..64ac758 100644
--- a/R/MCMCpanel.R
+++ b/R/MCMCpanel.R
@@ -241,7 +241,7 @@
output <- mcmc(data=sample, start=1, end=mcmc, thin=thin)
varnames(output) <- names
attr(output,"title") <-
- "MCMCpack Linear Panel Model Posterior Density Sample"
+ "MCMCpack Linear Panel Model Posterior Sample"
return(output)
}
diff --git a/R/MCMCpoisson.R b/R/MCMCpoisson.R
index 3eecf13..16253e7 100644
--- a/R/MCMCpoisson.R
+++ b/R/MCMCpoisson.R
@@ -5,15 +5,29 @@
## KQ 3/17/2003 [bug fix]
## Modified to meet new developer specification 7/15/2004 KQ
## Modified for new Scythe and rngs 7/26/2004 KQ
+## Modified to handle marginal likelihood calculation 1/27/2006 KQ
"MCMCpoisson" <-
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,
+ marginal.likelihood = c("none", "Laplace"),...) {
## checks
check.offset(list(...))
check.mcmc.parameters(burnin, mcmc, thin)
+
+ cl <- match.call()
+
+ ## get marginal likelihood argument
+ marginal.likelihood <- match.arg(marginal.likelihood)
+ if (max(B0==0) == 1){
+ if (marginal.likelihood != "none"){
+ warning("Cannot calculate marginal likelihood with improper prior\n")
+ marginal.likelihood <- "none"
+ }
+ }
+ logmarglike <- NULL
## seeds
seeds <- form.seeds(seed)
@@ -44,12 +58,34 @@
stop("\n Check data and call MCMCpoisson() again. \n")
}
- ## define holder for posterior density sample
+ ## define holder for posterior sample
sample <- matrix(data=0, mcmc/thin, dim(X)[2] )
-
+
+
+
+ ## marginal likelihood calculation if Laplace
+ if (marginal.likelihood == "Laplace"){
+ theta.start <- beta.start
+ optim.out <- optim(theta.start, logpost.poisson, method="BFGS",
+ control=list(fnscale=-1),
+ hessian=TRUE, y=Y, X=X, b0=b0, B0=B0)
+
+ theta.tilde <- optim.out$par
+ beta.tilde <- theta.tilde[1:K]
+
+ Sigma.tilde <- solve(-1*optim.out$hessian)
+
+ logmarglike <- (length(theta.tilde)/2)*log(2*pi) +
+ log(sqrt(det(Sigma.tilde))) +
+ logpost.poisson(theta.tilde, Y, X, b0, B0)
+
+ }
+
+
## call C++ code to draw sample
auto.Scythe.call(output.object="posterior", cc.fun.name="MCMCpoisson",
- sample=sample, Y=Y, X=X, burnin=as.integer(burnin),
+ sample.nonconst=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),
@@ -59,7 +95,8 @@
## put together matrix and build MCMC object to return
output <- form.mcmc.object(posterior, names=xnames,
- title="MCMCpoisson Posterior Density Sample")
+ title="MCMCpoisson Posterior Sample",
+ y=Y, call=cl, logmarglike=logmarglike)
return(output)
}
diff --git a/R/MCMCprobit.R b/R/MCMCprobit.R
index 2de612e..dde2e9a 100644
--- a/R/MCMCprobit.R
+++ b/R/MCMCprobit.R
@@ -4,16 +4,30 @@
## ADM and KQ 5/21/2002
## Modified to meet new developer specification 7/26/2004 KQ
## Modified for new Scythe and rngs 7/26/2004 KQ
+## Modified to handle marginal likelihood calculation 1/27/2006 KQ
"MCMCprobit" <-
function(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
thin = 1, verbose = 0, seed = NA, beta.start = NA,
- b0 = 0, B0 = 0, bayes.resid=FALSE, ...) {
+ b0 = 0, B0 = 0, bayes.resid=FALSE,
+ marginal.likelihood = c("none", "Laplace"), ...) {
## checks
check.offset(list(...))
check.mcmc.parameters(burnin, mcmc, thin)
+
+ cl <- match.call()
+
+ ## get marginal likelihood argument
+ marginal.likelihood <- match.arg(marginal.likelihood)
+ if (max(B0==0) == 1){
+ if (marginal.likelihood != "none"){
+ warning("Cannot calculate marginal likelihood with improper prior\n")
+ marginal.likelihood <- "none"
+ }
+ }
+ logmarglike <- NULL
## seeds
seeds <- form.seeds(seed)
@@ -25,7 +39,7 @@
holder <- parse.formula(formula, data)
Y <- holder[[1]]
X <- holder[[2]]
- xnames <- holder[[3]]
+ xnames <- holder[[3]]
K <- ncol(X) # number of covariates
## starting values and priors
@@ -53,14 +67,34 @@
cat("Elements of Y equal to something other than 0 or 1.\n")
stop("Check data and call MCMCprobit() again.\n")
}
-
+
+
+ ## marginal likelihood calculation if Laplace
+ if (marginal.likelihood == "Laplace"){
+ theta.start <- beta.start
+ optim.out <- optim(theta.start, logpost.probit, method="BFGS",
+ control=list(fnscale=-1),
+ hessian=TRUE, y=Y, X=X, b0=b0, B0=B0)
+
+ theta.tilde <- optim.out$par
+ beta.tilde <- theta.tilde[1:K]
+
+ Sigma.tilde <- solve(-1*optim.out$hessian)
+
+ logmarglike <- (length(theta.tilde)/2)*log(2*pi) +
+ log(sqrt(det(Sigma.tilde))) +
+ logpost.probit(theta.tilde, Y, X, b0, B0)
+
+ }
+
if (is.null(resvec)){
## define holder for posterior density sample
sample <- matrix(data=0, mcmc/thin, dim(X)[2] )
## call C++ code to draw sample
auto.Scythe.call(output.object="posterior", cc.fun.name="MCMCprobit",
- sample=sample, Y=Y, X=X, burnin=as.integer(burnin),
+ sample.nonconst=sample, Y=Y, X=X,
+ burnin=as.integer(burnin),
mcmc=as.integer(mcmc), thin=as.integer(thin),
lecuyer=as.integer(lecuyer),
seedarray=as.integer(seed.array),
@@ -70,7 +104,8 @@
## put together matrix and build MCMC object to return
output <- form.mcmc.object(posterior, names=xnames,
- title="MCMCprobit Posterior Density Sample")
+ title="MCMCprobit Posterior Sample",
+ y=Y, call=cl, logmarglike=logmarglike)
}
else{
@@ -79,7 +114,7 @@
## call C++ code to draw sample
auto.Scythe.call(output.object="posterior", cc.fun.name="MCMCprobitres",
- sample=sample, Y=Y, X=X, resvec=resvec,
+ sample.nonconst=sample, Y=Y, X=X, resvec=resvec,
burnin=as.integer(burnin),
mcmc=as.integer(mcmc), thin=as.integer(thin),
lecuyer=as.integer(lecuyer),
@@ -92,7 +127,8 @@
xnames <- c(xnames, paste("epsilonstar", as.character(resvec), sep="") )
output <- form.mcmc.object(posterior, names=xnames,
- title="MCMCprobit Posterior Density Sample")
+ title="MCMCprobit Posterior Sample",
+ y=Y, call=cl, logmarglike=logmarglike)
}
return(output)
diff --git a/R/MCMCregress.R b/R/MCMCregress.R
index d22cb8a..c5a084a 100644
--- a/R/MCMCregress.R
+++ b/R/MCMCregress.R
@@ -1,44 +1,64 @@
-# MCMCregress.R samples from the posterior distribution of a Gaussian
-# linear regression model in R using linked C++ code in Scythe
-#
-# Original written by ADM and KQ 5/21/2002
-# Updated with helper functions ADM 5/28/2004
-# Modified to meet new developer specification 6/18/2004 KQ
-# Modified for new Scythe and rngs 7/22/2004 ADM
+## MCMCregress.R samples from the posterior distribution of a Gaussian
+## linear regression model in R using linked C++ code in Scythe
+##
+## Original written by ADM and KQ 5/21/2002
+## Updated with helper functions ADM 5/28/2004
+## Modified to meet new developer specification 6/18/2004 KQ
+## Modified for new Scythe and rngs 7/22/2004 ADM
+## Modified to handle marginal likelihood calculation 1/26/2006 KQ
"MCMCregress" <-
function(formula, data=parent.frame(), burnin = 1000, mcmc = 10000,
- thin=1, verbose = 0, seed = NA, beta.start = NA,
- b0 = 0, B0 = 0, c0 = 0.001, d0 = 0.001, ...) {
+ thin=1, verbose = 0, seed = NA, beta.start = NA,
+ b0 = 0, B0 = 0, c0 = 0.001, d0 = 0.001,
+ marginal.likelihood = c("none", "Laplace", "Chib95"),
+ ...) {
- # checks
+ ## checks
check.offset(list(...))
check.mcmc.parameters(burnin, mcmc, thin)
+
+ cl <- match.call()
+
+
+ ## get marginal likelihood argument
+ marginal.likelihood <- match.arg(marginal.likelihood)
+ if (max(B0==0) == 1){
+ if (marginal.likelihood != "none"){
+ warning("Cannot calculate marginal likelihood with improper prior\n")
+ marginal.likelihood <- "none"
+ }
+ }
+ logmarglike <- NULL
+ chib <- 0
+ if (marginal.likelihood == "Chib95"){
+ chib <- 1
+ }
- # seeds
+ ## seeds
seeds <- form.seeds(seed)
lecuyer <- seeds[[1]]
seed.array <- seeds[[2]]
lecuyer.stream <- seeds[[3]]
- # form response and model matrices
+ ## form response and model matrices
holder <- parse.formula(formula, data)
Y <- holder[[1]]
X <- holder[[2]]
xnames <- holder[[3]]
K <- ncol(X) # number of covariates
- # starting values and priors
+ ## starting values and priors
beta.start <- coef.start(beta.start, K, formula, family=gaussian, data)
mvn.prior <- form.mvn.prior(b0, B0, K)
b0 <- mvn.prior[[1]]
B0 <- mvn.prior[[2]]
check.ig.prior(c0, d0)
- # define holder for posterior density sample
+ ## define holder for posterior sample
sample <- matrix(data=0, mcmc/thin, K+1)
- # call C++ code to draw sample
+ ## call C++ code to draw sample
auto.Scythe.call(output.object="posterior", cc.fun.name="MCMCregress",
sample=sample, Y=Y, X=X, burnin=as.integer(burnin),
mcmc=as.integer(mcmc), thin=as.integer(thin),
@@ -46,11 +66,39 @@
seedarray=as.integer(seed.array),
lecuyerstream=as.integer(lecuyer.stream),
verbose=as.integer(verbose), betastart=beta.start,
- b0=b0, B0=B0, c0=as.double(c0), d0=as.double(d0))
-
- # pull together matrix and build MCMC object to return
+ b0=b0, B0=B0, c0=as.double(c0), d0=as.double(d0),
+ logmarglikeholder.nonconst=as.double(0.0),
+ chib=as.integer(chib))
+
+ ## get marginal likelihood if Chib95
+ if (marginal.likelihood == "Chib95"){
+ logmarglike <- posterior$logmarglikeholder
+ }
+ ## marginal likelihood calculation if Laplace
+ if (marginal.likelihood == "Laplace"){
+ theta.start <- c(beta.start, log(0.5*var(Y)))
+ optim.out <- optim(theta.start, logpost.regress, method="BFGS",
+ control=list(fnscale=-1),
+ hessian=TRUE, y=Y, X=X, b0=b0, B0=B0, c0=c0, d0=d0)
+
+ theta.tilde <- optim.out$par
+ beta.tilde <- theta.tilde[1:K]
+ sigma2.tilde <- exp(theta.tilde[K+1])
+
+ Sigma.tilde <- solve(-1*optim.out$hessian)
+
+ logmarglike <- (length(theta.tilde)/2)*log(2*pi) +
+ log(sqrt(det(Sigma.tilde))) +
+ logpost.regress(theta.tilde, Y, X, b0, B0, c0, d0)
+
+ }
+
+ ## pull together matrix and build MCMC object to return
output <- form.mcmc.object(posterior,
names=c(xnames, "sigma2"),
- title="MCMCregress Posterior Density Sample")
+ title="MCMCregress Posterior Sample",
+ y=Y, call=cl,
+ logmarglike=logmarglike
+ )
return(output)
}
diff --git a/R/MCMCtobit.R b/R/MCMCtobit.R
index 633f3e4..ba35697 100644
--- a/R/MCMCtobit.R
+++ b/R/MCMCtobit.R
@@ -42,12 +42,12 @@ function(formula, data=parent.frame(), below = 0.0, above = Inf,
B0 <- mvn.prior[[2]]
check.ig.prior(c0, d0)
- # define holder for posterior density sample
+ # define holder for posterior sample
sample <- matrix(data=0, mcmc/thin, K+1)
# call C++ code to draw sample
auto.Scythe.call(output.object="posterior", cc.fun.name="MCMCtobit",
- sample=sample, Y=Y, X=X, below=as.double(below),
+ sample.nonconst=sample, Y=Y, X=X, below=as.double(below),
above=as.double(above),
burnin=as.integer(burnin), mcmc=as.integer(mcmc),
thin=as.integer(thin),
@@ -60,6 +60,6 @@ function(formula, data=parent.frame(), below = 0.0, above = Inf,
# pull together matrix and build MCMC object to return
output <- form.mcmc.object(posterior,
names=c(xnames, "sigma2"),
- title="MCMCtobit Posterior Density Sample")
+ title="MCMCtobit Posterior Sample")
return(output)
}
diff --git a/R/MCmodels.R b/R/MCmodels.R
new file mode 100644
index 0000000..fca6c63
--- /dev/null
+++ b/R/MCmodels.R
@@ -0,0 +1,109 @@
+## Monte Carlo simulation from the likelihood of a
+## binomial distribution with a Beta(alpha, beta) prior
+## ADM 1/25/2006
+MCbinomialbeta <- function(y, n, alpha=1, beta=1, mc=1000, ...) {
+
+ # check data
+ if(y < 0) {
+ cat("Error: Number of successes negative.\n")
+ stop("Please respecify and call function again.")
+ }
+ if(n < 0) {
+ cat("Error: Number of trials negative.\n")
+ stop("Please respecify and call function again.")
+ }
+ if(y > n) {
+ cat("Error: Number of successes greater than number of trials.\n")
+ stop("Please respecify and call function again.")
+ }
+
+ # check other parameters
+ check.beta.prior(alpha, beta)
+ check.mc.parameter(mc)
+
+ # draw sample and return
+ output <- mcmc(matrix(rbeta(mc, alpha+y, beta+n-y),mc,1))
+ varnames(output) <- as.list("pi")
+ attr(output,"title") <- "MCbinomialbeta Posterior Sample"
+
+ return(output)
+}
+
+## Monte Carlo simulation from the likelihood of a
+## Poisson distribution with a Gamma(alpha, beta) prior
+## ADM 1/25/2006
+MCpoissongamma <- function(y, alpha, beta, mc=1000, ...) {
+
+ # check data
+ if(any(y < 0)) {
+ cat("Error: Some counts negative in y.\n")
+ stop("Please respecify and call function again.")
+ }
+ n <- length(y)
+
+ # check other parameters
+ check.gamma.prior(alpha, beta)
+ check.mc.parameter(mc)
+
+ # draw sample and return
+ output <- mcmc(matrix(rgamma(mc, alpha+sum(y), beta+n),mc,1))
+ varnames(output) <- as.list("lambda")
+ attr(output,"title") <- "MCpoissongamma Posterior Sample"
+
+ return(output)
+}
+
+## Monte Carlo simulation from the likelihood of a
+## Normal distribution with a Normal(mu0, tau20) prior
+## the variance sigma2 is known
+## ADM 1/26/2006
+MCnormalnormal <- function(y, sigma2, mu0, tau20, mc=1000, ...) {
+
+ n <- length(y)
+ if(sigma2 <= 0) {
+ cat("Error: Known variance sigma2 is less than or equal to zero.\n")
+ stop("Please respecify and call function again.")
+ }
+
+ # check other parameters
+ check.normal.prior(mu0, tau20)
+ check.mc.parameter(mc)
+
+ # draw sample and return
+ mu1 = (1/tau20 * mu0 + n/sigma2 * mean(y)) / (1/tau20 + n/sigma2)
+ tau21 = 1 / (1/tau20 + n/sigma2)
+ output <- mcmc(matrix(rnorm(mc, mu1, sqrt(tau21)),mc,1))
+ varnames(output) <- as.list("mu")
+ attr(output,"title") <- "MCnormalnormal Posterior Sample"
+
+ return(output)
+}
+
+## Monte Carlo simulation from the likelihood of a
+## multinomal distribution with a Dirichlet(alpha) prior
+MCmultinomdirichlet <- function(y, alpha0, mc=1000, ...) {
+
+ # check data
+ d <- length(y)
+ if(any(y < 0)) {
+ cat("Error: Some counts negative in y.\n")
+ stop("Please respecify and call function again.")
+ }
+
+ # check alpha
+ if(length(alpha0) != d) {
+ cat("Error: Dimension of alpha and y do not match.\n")
+ stop("Please respecify and call function again.")
+ }
+ if(any(alpha0 <= 0)) {
+ cat("Error: At least one alpha in Dirichlet prior less than or equal to zero.\n")
+ stop("Please respecify and call function again.")
+ }
+
+ # draw sample and return
+ output <- mcmc(rdirichlet(mc,y + alpha0))
+ varnames(output) <- paste("pi.", 1:d, sep="")
+ attr(output,"title") <- "MCmultinomdirichlet Posterior Sample"
+
+ return(output)
+}
diff --git a/R/automate.R b/R/automate.R
index 46949b9..633953f 100644
--- a/R/automate.R
+++ b/R/automate.R
@@ -31,10 +31,17 @@
## Matrices will be, for example, Xdata, Xrow, and
## Xcol.
##
+## Any objects that are changed in the C++ code
+## need to have C++ names like sample.nonconst. All
+## *.nonconst have the tag stripped and a pointer
+## is used to change values. This is used in all models
+## for the posterior density sample (sample.nonconst),
+## and other quantities of interest.
+##
## This also build a skeleton C++ template and clean R template
-## for MCMCpack if developer=TRUE. All are arguments constants
-## except the matrix named "sample".
+## for MCMCpack if developer=TRUE.
##
+## Updated by ADM and KQ 1/25/2006 (to allow for multiple nonconsts)
"auto.Scythe.call" <-
function(output.object, cc.fun.name, package="MCMCpack",
@@ -44,9 +51,15 @@
objects <- list(...)
K <- length(objects)
c.names <- names(objects)
- if (max(c.names=="sample") != 1){
- stop("argument named 'sample' must be specified in call to auto.Scythe.call()\n")
- }
+ nonconst.indic <- rep(FALSE, K)
+ test <- grep("\.nonconst$", c.names)
+ if(length(test)==0)
+ stop("something must be declared non-constant in auto.Scythe.call()\n")
+ nonconst.indic[test] <- TRUE
+ c.names <- sub("\.nonconst$", "", c.names)
+ if(length(unique(c.names)) != K)
+ stop("non-unique nonconst names passed in auto.Scythe.call()\n")
+
R.names <- strsplit(toString(match.call()), ",")[[1]]
R.names <- R.names[(length(R.names)-K+1):length(R.names)]
@@ -80,8 +93,6 @@
objects <- c(objects, verbose=as.integer(FALSE))
R.names[length(objects)] <- "as.integer(FALSE)"
}
- K <- length(objects)
- c.names <- names(objects)
## check parameters
check.mcmc.parameters(objects$burnin, objects$mcmc, objects$thin)
@@ -139,8 +150,6 @@
stop("Integers, doubles, or matrices only to auto.Scythe.call().")
}
}
-
-
# clean up and return R call
middle <- paste(middle, sep=" ", collapse=" ")
@@ -159,42 +168,32 @@
for(k in 1:K) {
if(is.double(objects[[k]]) & is.null(dim(objects[[k]]))) {
- holder <- paste("const double *", c.names[[k]], ",", sep="")
+ holder <- paste("double *", c.names[[k]], ",", sep="")
+ if(!nonconst.indic[k]) holder <- paste("const ", holder, sep="")
c.middle <- c(c.middle, holder)
}
else if(is.integer(objects[[k]]) & is.null(dim(objects[[k]]))) {
- holder <- paste("const int *", c.names[[k]], ",", sep="")
+ holder <- paste("int *", c.names[[k]], ",", sep="")
+ if(!nonconst.indic[k]) holder <- paste("const ", holder, sep="")
c.middle <- c(c.middle, holder)
}
else if(is.matrix(objects[[k]])) {
-
- # pull together Scythe call [note sample]
- if(c.names[[k]]=="sample") {
- holder.data <- paste("double *", c.names[[k]], "data,", sep="")
- scythe <- NULL
- together.call <- paste(together.call, scythe, sep="")
- sample.block <- paste(" const int size = *", c.names[[k]],
- "row * *", c.names[[k]],
- "col;\n for(int i = 0; i < size; ++i)\n",
- " ", c.names[[k]], "data[i] = STORAGEMATRIX[i];\n",
- sep="")
- }
- else {
- holder.data <- paste("const double *", c.names[[k]],
- "data,", sep="")
- scythe <- paste(" Matrix <double> ", c.names[[k]],
- " = r2scythe(*",
- c.names[[k]], "row, *", c.names[[k]], "col, ",
- c.names[[k]], "data);\n", sep="")
+ holder.data <- paste("double *", c.names[[k]],
+ "data,", sep="")
+ if(!nonconst.indic[k]) holder.data <-
+ paste("const ", holder.data, sep="")
+ scythe <- paste(" Matrix <double> ", c.names[[k]],
+ " = r2scythe(*",
+ c.names[[k]], "row, *", c.names[[k]], "col, ",
+ c.names[[k]], "data);\n", sep="")
- together.call <- paste(together.call, scythe, sep="")
+ together.call <- paste(together.call, scythe, sep="")
}
- holder.row <- paste("const int *", c.names[[k]], "row,", sep="")
- holder.col <- paste("const int *", c.names[[k]], "col,", sep="")
- c.middle <- c(c.middle, holder.data, holder.row, holder.col)
+ holder.row <- paste("const int *", c.names[[k]], "row,", sep="")
+ holder.col <- paste("const int *", c.names[[k]], "col,", sep="")
+ c.middle <- c(c.middle, holder.data, holder.row, holder.col)
}
- }
-
+
# clean up and print C++ function call
c.middle <- paste(c.middle, sep=" ", collapse=" ")
c.call <- paste(c.start, c.middle, c.end, sep="")
@@ -230,7 +229,7 @@
cat(comment.block, includes.block, main.block, function.call,
together.block, together.call, constants.block, storage.block,
seed.block, startval.block, sample.call,
- sample.block, end.block, sep="", file=cc.file)
+ end.block, sep="", file=cc.file)
if (cc.file != "") {
cat("\nCreated file named '", cc.file, "'.\n", sep="")
cat("Edit the file and move it to the appropriate directory.\n")
diff --git a/R/distn.R b/R/distn.R
index bfa10f6..1d7b592 100644
--- a/R/distn.R
+++ b/R/distn.R
@@ -277,10 +277,10 @@
"rdirichlet" <-
function(n, alpha) {
- l<-length(alpha);
- x<-matrix(rgamma(l*n,alpha),ncol=l,byrow=TRUE);
- sm<-x%*%rep(1,l);
- return(x/as.vector(sm));
+ l <- length(alpha)
+ x <- matrix(rgamma(l*n,alpha),ncol=l,byrow=TRUE)
+ sm <- x%*%rep(1,l)
+ return(x/as.vector(sm))
}
##
diff --git a/R/hidden.R b/R/hidden.R
index 486a6ad..1a262a2 100644
--- a/R/hidden.R
+++ b/R/hidden.R
@@ -5,7 +5,8 @@
## fixes them in all functions simultaneously.
##
## updated by ADM 7/22/04
-## re-organized (alphabetical) by adm 7/28/04
+## re-organized (alphabetical) by ADM 7/28/04
+## added a number of functions for teaching models by ADM 1/25/2006
## create an agreement score matrix from a vote matrix
## subjects initially on rows and items on cols of X
@@ -101,6 +102,66 @@
return(0)
}
+# check inverse Gamma prior
+"check.ig.prior" <-
+ function(c0, d0) {
+
+ if(c0 <= 0) {
+ cat("Error: IG(c0/2,d0/2) prior c0 less than or equal to zero.\n")
+ stop("Please respecify and call ", calling.function(), " again.\n",
+ call.=FALSE)
+ }
+ if(d0 <= 0) {
+ cat("Error: IG(c0/2,d0/2) prior d0 less than or equal to zero.\n")
+ stop("Please respecify and call ", calling.function(), " again.\n",
+ call.=FALSE)
+ }
+ return(0)
+ }
+
+# check Gamma prior
+# ADM 1/25/2006
+"check.gamma.prior" <-
+ function(alpha, beta) {
+
+ if(alpha <= 0) {
+ cat("Error: Gamma(alpha,beta) prior alpha less than or equal to zero.\n")
+ stop("Please respecify and call ", calling.function(), " again.\n",
+ call.=FALSE)
+ }
+ if(alpha <= 0) {
+ cat("Error: Gamma(alpha,beta) prior beta less than or equal to zero.\n")
+ stop("Please respecify and call ", calling.function(), " again.\n",
+ call.=FALSE)
+ }
+ return(0)
+ }
+
+# check Normal prior
+# ADM 1/26/2006
+"check.normal.prior" <-
+ function(mu, sigma2) {
+
+ if(sigma2 <= 0) {
+ cat("Error: Normal(mu0,tau20) prior sigma2 less than or equal to zero.\n")
+ stop("Please respecify and call ", calling.function(), " again.\n",
+ call.=FALSE)
+ }
+ }
+
+# check mc parameter
+# ADM 1/25/2006
+"check.mc.parameter" <-
+ function(mc) {
+
+ if(mc < 0) {
+ cat("Error: Monte Carlo iterations negative.\n")
+ stop("Please respecify and call ", calling.function(), " again.",
+ call.=FALSE)
+ }
+ return(0)
+ }
+
# check mcmc parameters
"check.mcmc.parameters" <-
function(burnin, mcmc, thin) {
@@ -133,7 +194,8 @@
function(args) {
if(sum(names(args)=="offset")==1) {
cat("Error: Offsets are currently not supported in MCMCpack.\n")
- stop("Please respecify and call ", calling.function(), " again.\n")
+ stop("Please respecify and call ", calling.function(), " again.\n",
+ call.=FALSE)
}
return(0)
}
@@ -490,17 +552,29 @@
# pull together the posterior density sample
"form.mcmc.object" <-
- function(posterior.object, names, title) {
+ function(posterior.object, names, title, ...) {
holder <- matrix(posterior.object$sampledata,
posterior.object$samplerow,
posterior.object$samplecol,
byrow=TRUE)
-
+
+
output <- mcmc(data=holder, start=(posterior.object$burnin+1),
end=(posterior.object$burnin+posterior.object$mcmc),
thin=posterior.object$thin)
varnames(output) <- as.list(names)
attr(output,"title") <- title
+
+ attribs <- list(...)
+ K <- length(attribs)
+ attrib.names <- names(attribs)
+
+ if (K>0){
+ for (i in 1:K){
+ attr(output, attrib.names[i]) <- attribs[[i]]
+ }
+ }
+
return(output)
}
diff --git a/R/procrust.R b/R/procrust.R
new file mode 100644
index 0000000..d6bf7d8
--- /dev/null
+++ b/R/procrust.R
@@ -0,0 +1,64 @@
+## function that performs procrustes transformation on X with target Xstar
+##
+## returns the rotation matrix R, translation vector tt,
+## and dilation factor s for which:
+##
+## s X R + 1 tt' \approx Xstar
+##
+## along with X.new = s X R + 1 tt'
+##
+## Based on Borg and Groenen 1997. Modern Multidimensional
+## Scaling. New York: Springer. pp. 340-342.
+##
+## Kevin Quinn
+## Harvard University
+## 6/13/2005
+##
+procrustes <- function(X, Xstar, translation=FALSE, dilation=FALSE){
+ if (nrow(X) != nrow(Xstar)){
+ cat("X and Xstar do not have same number of rows.\n")
+ stop("Check data and call procrustes() again. \n")
+ }
+ if (ncol(X) != ncol(Xstar)){
+ cat("X and Xstar do not have same number of columns.\n")
+ stop("Check data and call procrustes() again. \n")
+ }
+
+ n <- nrow(X)
+ m <- ncol(X)
+
+ if (translation){
+ J <- diag(n) - 1/n * matrix(1, n, n)
+ }
+ else{
+ J <- diag(n)
+ }
+
+ C <- t(Xstar) %*% J %*% X
+ svd.out <- svd(C)
+ R <- svd.out$v %*% t(svd.out$u)
+ s <- 1
+ if (dilation){
+ mat1 <- t(Xstar) %*% J %*% X %*% R
+ mat2 <- t(X) %*% J %*% X
+ s.numer <- 0
+ s.denom <- 0
+ for (i in 1:m){
+ s.numer <- s.numer + mat1[i,i]
+ s.denom <- s.denom + mat2[i,i]
+ }
+ s <- s.numer / s.denom
+ }
+ tt <- matrix(0, m, 1)
+ if (translation){
+ tt <- 1/n * t(Xstar - s*X %*% R) %*% matrix(1, n, 1)
+ }
+
+ X.new <- s * X %*% R + matrix(tt, nrow(X), ncol(X), byrow=TRUE)
+
+ return(list(X.new=X.new, R=R, tt=tt, s=s))
+}
+
+
+
+
diff --git a/R/utility.R b/R/utility.R
index 425f011..b2515c6 100644
--- a/R/utility.R
+++ b/R/utility.R
@@ -21,10 +21,12 @@
#
# ADM 4/18/2003
# ADM 11/13/2003 [bug fix]
+# ADM 1/25/2006 [patch to automatically compute matrix size]
"xpnd" <-
- function (x, nrow) {
+ function (x, nrow = NULL) {
dim(x) <- NULL
+ if(is.null(nrow)) nrow <- (-1 + sqrt(1 + 8 * length(x))) / 2
output <- matrix(0, nrow, nrow)
output[lower.tri(output, diag = TRUE)] <- x
hold <- output
diff --git a/man/BayesFactor.Rd b/man/BayesFactor.Rd
new file mode 100644
index 0000000..d7def3b
--- /dev/null
+++ b/man/BayesFactor.Rd
@@ -0,0 +1,74 @@
+\name{BayesFactor}
+\alias{BayesFactor}
+\alias{is.BayesFactor}
+
+\title{Create an object of class BayesFactor from MCMCpack output}
+\description{
+This function creates an object of class \code{BayesFactor} from MCMCpack
+output.
+
+}
+\usage{
+BayesFactor(...)
+is.BayesFactor(BF)
+}
+
+\arguments{
+ \item{...}{MCMCpack output objects. These have to be of class
+ \code{mcmc} and have a \code{logmarglike} attribute. In what
+ follows, we let \code{M} denote the total number of models to be
+ compared.}
+
+ \item{BF}{An object to be checked for membership in class
+ \code{BayesFactor}.}
+}
+
+\value{
+ An object of class \code{BayesFactor}. A \code{BayesFactor} object has
+ four attributes. They are: \code{BF.mat} an \eqn{M \times M}{M by M}
+ matrix in which element \eqn{i,j}{i,j} contains the Bayes factor for
+ model \eqn{i}{i} relative to model \eqn{j}{j}; \code{BF.log.mat} an
+ \eqn{M \times M}{M by M} matrix in which element \eqn{i,j}{i,j}
+ contains the natural log of the Bayes factor for model \eqn{i}{i}
+ relative to model \eqn{j}{j}; \code{BF.logmarglike} an \eqn{M}{M}
+ vector containing the log marginal likelihoods for models 1 through
+ \eqn{M}{M}; and \code{BF.call} an \eqn{M}{M} element list containing
+ the calls used to fit models 1 through \eqn{M}{M}.
+}
+
+\examples{
+\dontrun{
+data(birthwt)
+
+model1 <- MCMCregress(bwt~age+lwt+as.factor(race) + smoke + ht,
+ data=birthwt, b0=c(2700, 0, 0, -500, -500,
+ -500, -500),
+ B0=c(1e-6, .01, .01, 1.6e-5, 1.6e-5, 1.6e-5,
+ 1.6e-5), c0=10, d0=4500000,
+ marginal.likelihood="Chib95", mcmc=10000)
+
+model2 <- MCMCregress(bwt~age+lwt+as.factor(race) + smoke,
+ data=birthwt, b0=c(2700, 0, 0, -500, -500,
+ -500),
+ B0=c(1e-6, .01, .01, 1.6e-5, 1.6e-5, 1.6e-5),
+ c0=10, d0=4500000,
+ marginal.likelihood="Chib95", mcmc=10000)
+
+model3 <- MCMCregress(bwt~as.factor(race) + smoke + ht,
+ data=birthwt, b0=c(2700, -500, -500,
+ -500, -500),
+ B0=c(1e-6, 1.6e-5, 1.6e-5, 1.6e-5,
+ 1.6e-5), c0=10, d0=4500000,
+ marginal.likelihood="Chib95", mcmc=10000)
+
+BF <- BayesFactor(model1, model2, model3)
+print(BF)
+
+}
+}
+
+\concept{Bayes factor}
+\concept{model comparison}
+
+\seealso{\code{\link{MCMCregress}}}
+\keyword{models}
diff --git a/man/MCMCSVDreg.Rd b/man/MCMCSVDreg.Rd
new file mode 100644
index 0000000..5f9a045
--- /dev/null
+++ b/man/MCMCSVDreg.Rd
@@ -0,0 +1,161 @@
+\name{MCMCSVDreg}
+\alias{MCMCSVDreg}
+\title{Markov Chain Monte Carlo for SVD Regression}
+\description{
+ This function generates a sample from the posterior distribution
+ of a linear regression model with Gaussian errors in which the
+ design matrix has been decomposed with singular value
+ decomposition.The sampling is done via the Gibbs sampling algorithm.
+ The user supplies data 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.
+ }
+
+\usage{
+MCMCSVDreg(formula, data=parent.frame(), burnin = 1000, mcmc = 10000,
+ thin=1, verbose = 0, seed = NA, tau2.start = 1,
+ g0 = 0, a0 = 0.001, b0 = 0.001, c0=2, d0=2, w0=1,
+ beta.samp=FALSE, intercept=TRUE, ...)}
+
+\arguments{
+ \item{formula}{Model formula. Predictions are returned for elements
+ of y that are coded as NA.}
+
+ \item{data}{Data frame.}
+
+ \item{burnin}{The number of burn-in iterations for the sampler.}
+
+ \item{mcmc}{The number of MCMC iterations after burnin.}
+
+ \item{thin}{The thinning interval used in the simulation. The number of
+ MCMC iterations must be divisible by this value.}
+
+ \item{verbose}{A switch which determines whether or not the progress of
+ the sampler is printed to the screen. If \code{verbose} is greater
+ than 0 the iteration number, the
+ \eqn{\beta}{beta} vector, and the error variance are printed to
+ the screen every \code{verbose}th iteration.}
+
+ \item{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 \code{rep(12345,6)} is used). The second element of
+ list is a positive substream number. See the MCMCpack
+ specification for more details.}
+
+ \item{tau2.start}{The starting values for the \eqn{\tau^2}{tau^2} vector. Can
+ be either a scalar or a vector. If a scalar is passed then that value
+ will be the starting value for all elements of \eqn{\tau^2}{tau^2}.}
+
+
+
+ \item{g0}{The prior mean of \eqn{\gamma}{gamma}. This can either be a
+ scalar or a
+ column vector with dimension equal to the number of gammas. If this
+ takes a scalar value, then that value will serve as the prior
+ mean for all of the betas.}
+
+
+ \item{a0}{\eqn{a_0/2}{a0/2} is the shape parameter for the inverse
+ Gamma prior on \eqn{\sigma^2}{sigma^2} (the variance of the
+ disturbances). The amount of information in the inverse Gamma prior
+ is something like that from \eqn{a_0}{a0} pseudo-observations.}
+
+ \item{b0}{\eqn{b_0/2}{b0/2} is the scale parameter for the
+ inverse Gamma prior on \eqn{\sigma^2}{sigma^2} (the variance of the
+ disturbances). In constructing the inverse Gamma prior,
+ \eqn{b_0}{b0} acts like the sum of squared errors from the
+ \eqn{a_0}{a0} pseudo-observations.}
+
+ \item{c0}{\eqn{c_0/2}{c0/2} is the shape parameter for the inverse
+ Gamma prior on \eqn{\tau_i^2}{tau[i]^2}.}
+
+ \item{d0}{\eqn{d_0/2}{d0/2} is the scale parameter for the
+ inverse Gamma prior on \eqn{\tau_i^2}{tau[i]^2}.}
+
+ \item{w0}{The prior probability that \eqn{\gamma_i = 0}{gamma[i] = 0}.
+ Can be either a scalar or an \eqn{N}{N} vector where \eqn{N}{N} is the
+ number of observations.}
+
+ \item{beta.samp}{Logical indicating whether the sampled elements of
+ beta should be stored and returned.}
+
+ \item{intercept}{Logical indicating whether the original design matrix
+ should include a constant term.}
+
+ \item{...}{further arguments to be passed}
+}
+
+\value{
+ An mcmc object that contains the posterior sample. This
+ object can be summarized by functions provided by the coda package.
+}
+
+\details{
+ The model takes the following form:
+ \deqn{y = X \beta + \varepsilon}{y = X beta + epsilon}
+ Where the errors are assumed to be iid Gaussian:
+ \deqn{\varepsilon_{i} \sim \mathcal{N}(0, \sigma^2)}{epsilon_i ~ N(0,
+ sigma^2)}.
+
+ Let \eqn{N}{N} denote the number of rows of \eqn{X}{X} and \eqn{P}{P}
+ the number of columns of \eqn{X}{X}. Unlike the standard regression
+ setup where \eqn{N >> P}{N >> P} here it is the case that \eqn{P >>
+ N}{P >> N}.
+
+ To deal with this problem a singular value decomposition of
+ \eqn{X'}{X'} is performed: \eqn{X' = ADF}{X' = ADF} and the regression
+ model becomes
+
+ \deqn{y = F'D \gamma + \varepsilon}{y = F'D gamma + epsilon}
+
+ where \eqn{\gamma = A' \beta}{gamma = A' beta}.
+
+ We assume the following priors:
+
+ \deqn{\sigma^{-2} \sim \mathcal{G}amma(a_0/2, b_0/2)}{sigma^(-2) ~
+ Gamma(a0/2, b0/2)}
+
+ \deqn{\tau^{-2} \sim \mathcal{G}amma(c_0/2, d_0/2)}{tau^(-2) ~
+ Gamma(c0/2, d0/2)}
+
+ \deqn{\gamma_i \sim w0_i \delta_0 + (1-w0_i) \mathcal{N}(g0_i,
+ \sigma^2 \tau^2_i / d_i^2)}{
+ gamma[i] ~ w0[i] delta0 + (1-w0[i] N(g0[i], sigma^2 tau[i]^2/ d[i]^2)}
+
+ where \eqn{\delta_0}{delta0} is a unit point mass at 0 and
+ \eqn{d_i}{d[i]} is the \eqn{i}{i}th diagonal element of \eqn{D}{D}.
+
+ }
+
+ \references{
+ West, Mike; Josheph Nevins; Jeffrey Marks; Rainer Spang, and Harry
+ Zuzan. 2000. "DNA MICROARRAY DATA ANALYSIS AND REGRESSION MODELING
+ FOR GENETIC EXPRESSION PROFILING" Duke ISDS working paper.
+
+ Gottardo, Raphael, and Adrian Raftery. 2004. "Markov chain Monte
+ Carlo with mixtures of singular distributions". Statistics
+ department, University of Washington, Technical report 470.
+
+ Andrew D. Martin, Kevin M. Quinn, and Daniel Pemstein. 2004.
+ \emph{Scythe Statistical Library 1.0.} \url{http://scythe.wustl.edu}.
+
+ Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2002.
+ \emph{Output Analysis and Diagnostics for MCMC (CODA)}.
+ \url{http://www-fis.iarc.fr/coda/}.
+}
+
+
+\examples{
+\dontrun{
+}
+}
+
+\keyword{models}
+
+\seealso{\code{\link[coda]{plot.mcmc}},
+ \code{\link[coda]{summary.mcmc}}, \code{\link[stats]{lm}}}
diff --git a/man/MCMCirt1d.Rd b/man/MCMCirt1d.Rd
index 9077a77..8fbc827 100644
--- a/man/MCMCirt1d.Rd
+++ b/man/MCMCirt1d.Rd
@@ -21,7 +21,8 @@
MCMCirt1d(datamatrix, theta.constraints=list(), burnin = 1000,
mcmc = 20000, thin=1, verbose = 0, seed = NA, theta.start = NA,
alpha.start = NA, beta.start = NA, t0 = 0, T0 = 1, ab0=0, AB0=.25,
- store.item = FALSE, drop.constant.items=TRUE, ... ) }
+ store.item = FALSE, store.ability = TRUE,
+ drop.constant.items=TRUE, ... ) }
\arguments{
\item{datamatrix}{The matrix of data. Must be 0, 1, or missing values.
@@ -108,10 +109,19 @@ MCMCirt1d(datamatrix, theta.constraints=list(), burnin = 1000,
\item{store.item}{A switch that determines whether or not to
store the item parameters for posterior analysis.
- \emph{NOTE: This takes an enormous amount of memory, so
- should only be used if the chain is thinned heavily, or for
- applications with a small number of items}. By default, the
- item parameters are not stored.}
+ \emph{NOTE: In situations with many items storing the item
+ parameters takes an enormous amount of memory, so
+ \code{store.item} should only be \code{FALSE} if the chain is thinned
+ heavily, or for applications with a small number of items}.
+ By default, the item parameters are not stored.}
+
+ \item{store.ability}{A switch that determines whether or not to
+ store the ability parameters for posterior analysis.
+ \emph{NOTE: In situations with many individuals storing the ability
+ parameters takes an enormous amount of memory, so
+ \code{store.ability} should only be \code{TRUE} if the chain is thinned
+ heavily, or for applications with a small number of individuals}.
+ By default, the item parameters are stored.}
\item{drop.constant.items}{A switch that determines whether or not
items that have no variation
diff --git a/man/MCMCirtKdRob.Rd b/man/MCMCirtKdRob.Rd
new file mode 100644
index 0000000..af56250
--- /dev/null
+++ b/man/MCMCirtKdRob.Rd
@@ -0,0 +1,358 @@
+\name{MCMCirtKdRob}
+\alias{MCMCirtKdRob}
+\title{Markov Chain Monte Carlo for Robust K-Dimensional Item Response Theory
+ Model}
+\description{
+ This function generates a posterior sample from a
+ Robust K-dimensional item response theory (IRT) model with logistic
+ link, independent standard normal priors on the subject
+ abilities (ideal points), and independent normal priors on the
+ item parameters. The user supplies data 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.
+}
+
+\usage{
+MCMCirtKdRob(datamatrix, dimensions, item.constraints=list(),
+ ability.constraints=list(), burnin = 500, mcmc = 5000, thin=1,
+ interval.method="step", theta.w=0.5, theta.mp=4,
+ alphabeta.w=1.0, alphabeta.mp=4, delta0.w=NA, delta0.mp=3,
+ delta1.w=NA, delta1.mp=3, verbose = FALSE, seed = NA,
+ theta.start = NA, alphabeta.start = NA,
+ delta0.start = NA, delta1.start = NA,
+ b0 = 0, B0=0, k0=.1, k1=.1, c0=1, d0=1,
+ c1=1, d1=1, store.item=TRUE, store.ability=FALSE,
+ drop.constant.items=TRUE, ... ) }
+
+\arguments{
+ \item{datamatrix}{The matrix of data. Must be 0, 1, or missing values.
+ It is of dimensionality subjects by items.}
+
+ \item{dimensions}{The number of dimensions in the latent space.}
+
+ \item{item.constraints}{List of lists specifying possible equality
+ or simple inequality constraints on the item parameters. A typical
+ entry in the list has one of three forms: \code{rowname=list(d,c)}
+ which will constrain the dth item parameter for the item named
+ rowname to be equal to c, \code{rowname=list(d,"+")} which will
+ constrain the dth item parameter for the item named rowname to be
+ positive, and \code{rowname=list(d, "-")} which will constrain the dth
+ item parameter for the item named rowname to be negative. If
+ datamatrix is a
+ matrix without row names defaults names of ``V1", ``V2", ... , etc
+ will be used. In a \eqn{K}{K}-dimensional model,
+ the first item parameter for
+ item \eqn{i}{i} is the difficulty parameter (\eqn{\alpha_i}{alpha_i}),
+ the second item parameter is the discrimation parameter on dimension
+ 1 (\eqn{\beta_{i,1}}{beta_{i,1}}), the third item parameter is the
+ discrimation parameter on dimension 2
+ (\eqn{\beta_{i,2}}{beta_{i,2}}), ..., and the \eqn{(K+1)}{(K+1)}th
+ item parameter
+ is the discrimation parameter on dimension \eqn{K}{K}
+ (\eqn{\beta_{i,K}}{beta_{i,K}}).
+ The item difficulty parameters (\eqn{\alpha}{alpha}) should
+ generally not be constrained.
+ }
+
+ \item{ability.constraints}{List of lists specifying possible equality
+ or simple inequality constraints on the ability parameters. A typical
+ entry in the list has one of three forms: \code{colname=list(d,c)}
+ which will constrain the dth ability parameter for the subject named
+ colname to be equal to c, \code{colname=list(d,"+")} which will
+ constrain the dth ability parameter for the subject named colname to be
+ positive, and \code{colname=list(d, "-")} which will constrain the dth
+ ability parameter for the subject named colname to be negative. If
+ datamatrix is a
+ matrix without column names defaults names of ``V1", ``V2", ... , etc
+ will be used.}
+
+ \item{burnin}{The number of burn-in iterations for the sampler.}
+
+ \item{mcmc}{The number of iterations for the sampler after burn-in.}
+
+ \item{thin}{The thinning interval used in the simulation. The number of
+ iterations must be divisible by this value.}
+
+ \item{interval.method}{Method for finding the slicing interval. Can
+ be equal to either \code{step} in which case the stepping out
+ algorithm of Neal (2003) is used or \code{doubling} in which case the
+ doubling procedure of Neal (2003) is used. The stepping out method
+ tends to be faster on a per-iteration basis as it typically requires few
+ function calls. The doubling method expands the initial interval
+ more quickly which makes the Markov chain mix somewhat more
+ quickly-- at least in some situations. }
+
+ \item{theta.w}{The initial width of the slice sampling interval for
+ each ability parameter (the elements of \eqn{\theta}{theta})}
+
+ \item{theta.mp}{The parameter governing the maximum possible width of
+ the slice interval for each ability parameter (the elements of
+ \eqn{\theta}{theta}). If \code{interval.method="step"} then the maximum
+ width is \code{theta.w * theta.mp}.
+
+ If \code{interval.method="doubling"}
+ then the maximum width is \code{theta.w * 2^theta.mp}. }
+
+ \item{alphabeta.w}{The initial width of the slice sampling interval for
+ each item parameter (the elements of \eqn{\alpha}{alpha} and
+ \eqn{\beta}{beta})}
+
+ \item{alphabeta.mp}{ The parameter governing the maximum possible width of
+ the slice interval for each item parameters (the elements of
+ \eqn{\alpha}{alpha} and \eqn{\beta}{beta}). If
+ \code{interval.method="step"} then the maximum width is
+ \code{alphabeta.w * alphabeta.mp}.
+
+ If \code{interval.method="doubling"}
+ then the maximum width is \code{alphabeta.w * 2^alphabeta.mp}. }
+
+ \item{delta0.w}{The initial width of the slice sampling interval for
+ \eqn{\delta_0}{delta0}}
+
+ \item{delta0.mp}{The parameter governing the maximum possible width of
+ the slice interval for \eqn{\delta_0}{delta0}. If
+ \code{interval.method="step"} then the maximum width is
+ \code{delta0.w * delta0.mp}. If \code{interval.method="doubling"}
+ then the maximum width is \code{delta0.w * 2^delta0.mp}. }
+
+ \item{delta1.w}{The initial width of the slice sampling interval for
+ \eqn{\delta_1}{delta1}}
+
+ \item{delta1.mp}{The parameter governing the maximum possible width of
+ the slice interval for \eqn{\delta_1}{delta1}. If
+ \code{interval.method="step"} then the maximum width is
+ \code{delta1.w * delta1.mp}. If \code{interval.method="doubling"}
+ then the maximum width is \code{delta1.w * 2^delta1.mp}. }
+
+ \item{verbose}{A switch which determines whether or not the progress of
+ the sampler is printed to the screen. If verbose > 0, the
+ iteration number with be printed to the screen every verbose'th
+ iteration.}
+
+ \item{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 \code{rep(12345,6)} is used). The second element of
+ list is a positive substream number. See the MCMCpack
+ specification for more details.}
+
+ \item{theta.start}{The starting values for the ability parameters
+ \eqn{\theta}{theta}. Can be either a scalar or a matrix with
+ number of rows equal to the number of subjects and number of
+ columns equal to the dimension \eqn{K}{K} of the latent space. If
+ \code{theta.start=NA} then starting values will be chosen that are
+ 0 for unconstrained subjects, -0.5 for subjects with negative
+ inequality constraints and 0.5 for subjects with positive inequality
+ constraints. }
+
+ \item{alphabeta.start}{The starting values for the
+ \eqn{\alpha}{alpha} and \eqn{\beta}{beta} difficulty and
+ discrimination parameters. If \code{alphabeta.start} is set to a
+ scalar the starting value for all unconstrained item parameters will
+ be set to that scalar. If \code{alphabeta.start} is a matrix of
+ dimension \eqn{(K+1) \times items}{(K+1) x items} then the
+ \code{alphabeta.start} matrix is used as the starting values (except
+ for equality-constrained elements). If \code{alphabeta.start} is set
+ to \code{NA} (the default) then starting values for unconstrained
+ elements are set to values generated from a series of proportional
+ odds logistic regression fits, and starting values for inequality
+ constrained elements are set to either 1.0 or -1.0 depending on the
+ nature of the constraints. }
+
+ \item{delta0.start}{The starting value for the
+ \eqn{\delta_0}{delta0} parameter.}
+
+ \item{delta1.start}{The starting value for the
+ \eqn{\delta_1}{delta1} parameter.}
+
+ \item{b0}{The prior means of the
+ \eqn{\alpha}{alpha} and \eqn{\beta}{beta} difficulty and
+ discrimination parameters, stacked for all items.
+ If a scalar is passed, it
+ is used as the prior mean for all items.}
+
+ \item{B0}{The prior precisions (inverse variances) of the
+ independent Normal prior on the item parameters.
+ Can be either a scalar or a matrix of dimension
+ \eqn{(K+1) \times items}{(K+1) x items}.}
+
+ \item{k0}{\eqn{\delta_0}{delta0} is constrained to lie in the interval
+ between 0 and \code{k0}.}
+
+ \item{k1}{\eqn{\delta_1}{delta1} is constrained to lie in the interval
+ between 0 and \code{k1}.}
+
+ \item{c0}{Parameter governing the prior for
+ \eqn{\delta_0}{delta0}. \eqn{\delta_0}{delta0} divided by \code{k0} is
+ assumed to be follow a beta distribution with first parameter
+ \code{c0}.}
+
+ \item{d0}{Parameter governing the prior for
+ \eqn{\delta_0}{delta0}. \eqn{\delta_0}{delta0} divided by \code{k0} is
+ assumed to be follow a beta distribution with second parameter
+ \code{d0}.}
+
+ \item{c1}{Parameter governing the prior for
+ \eqn{\delta_1}{delta1}. \eqn{\delta_1}{delta1} divided by \code{k1} is
+ assumed to be follow a beta distribution with first parameter
+ \code{c1}.}
+
+ \item{d1}{Parameter governing the prior for
+ \eqn{\delta_1}{delta1}. \eqn{\delta_1}{delta1} divided by \code{k1} is
+ assumed to be follow a beta distribution with second parameter
+ \code{d1}.}
+
+ \item{store.item}{A switch that determines whether or not to
+ store the item parameters for posterior analysis.
+ \emph{NOTE: This typically takes an enormous amount of memory, so
+ should only be used if the chain is thinned heavily, or for
+ applications with a small number of items}. By default, the
+ item parameters are not stored.}
+
+ \item{store.ability}{A switch that determines whether or not to store
+ the subject abilities for posterior analysis. By default, the
+ item parameters are all stored.}
+
+ \item{drop.constant.items}{A switch that determines whether or not
+ items that have no variation
+ should be deleted before fitting the model. Default = TRUE.}
+
+ \item{...}{further arguments to be passed}
+}
+
+\value{
+ An mcmc object that contains the posterior sample. This
+ object can be summarized by functions provided by the coda package.
+}
+
+\details{
+ \code{MCMCirtKd} simulates from the posterior using
+ the slice sampling algorithm of Neal (2003).
+ 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. We assume that each
+ subject has an subject ability (ideal point) denoted
+ \eqn{\theta_j}{theta_j} \eqn{(K \times 1)}{(K x 1)},
+ and that each item has a scalar difficulty
+ parameter \eqn{\alpha_i}{alpha_i} and discrimination parameter
+ \eqn{\beta_i}{beta_i} \eqn{(K \times 1)}{(K x 1)}.
+ The observed choice by subject
+ \eqn{j}{j} on item \eqn{i}{i} is the observed data matrix which
+ is \eqn{(I \times J)}{(I * J)}.
+
+ The probability that subject \eqn{j}{j} answers item \eqn{i}{i}
+ correctly is assumed to be:
+
+ \deqn{\pi_{ij} = \delta_0 + (1 - \delta_0 - \delta_1)
+ /(1+\exp(\alpha_i - \beta_i \theta_j))}{pi_{ij} =
+ delta0 + (1 - delta0 - delta1) / (1 + exp(alpha_i - beta_i * theta_j))}
+
+ This model was discussed in Bafumi et al. (2005).
+
+ We assume the following priors. For the subject abilities (ideal points)
+ we assume independent standard Normal priors:
+ \deqn{\theta_{j,k} \sim \mathcal{N}(0,1)}{theta_j,k ~ N(0, 1)}
+ These cannot be changed by the user.
+ For each item parameter, we assume independent Normal priors:
+ \deqn{\left[\alpha_i, \beta_i \right]' \sim \mathcal{N}_{(K+1)}
+ (b_{0,i},B_{0,i})}{[alpha_i beta_i]' ~ N_(K+1) (b_0,i, B_0,i)}
+ Where \eqn{B_{0,i}}{B_0,i} is a diagonal matrix.
+ One can specify a separate prior mean and precision
+ for each item parameter. We also assume \eqn{\delta_0 / k_0 \sim
+ \mathcal{B}eta(c_0, d_0)}{delta0/k0 ~ Beta(c0, d0)} and
+ \eqn{\delta_1 / k_1 \sim
+ \mathcal{B}eta(c_1, d_1)}{delta1/k1 ~ Beta(c1, d1)}.
+
+ The model is identified by constraints on the item parameters and / or
+ ability parameters. See Rivers (2004) for a discussion of
+ identification of IRT models.
+
+ As is the case with all measurement models, make sure that you have plenty
+ of free memory, especially when storing the item parameters.
+}
+
+\references{
+ James H. Albert. 1992. ``Bayesian Estimation of Normal Ogive Item Response
+ Curves Using Gibbs Sampling." \emph{Journal of Educational Statistics}.
+ 17: 251-269.
+
+ Joseph Bafumi, Andrew Gelman, David K. Park, and Noah
+ Kaplan. 2005. ``Practical Issues in Implementing and Understanding
+ Bayesian Ideal Point Estimation.'' \emph{Political Analysis}.
+
+ Joshua Clinton, Simon Jackman, and Douglas Rivers. 2000. ``The Statistical
+ Analysis of Legislative Behavior: A Unified Approach." Paper presented at
+ the Annual Meeting of the Political Methodology Society.
+
+ Simon Jackman. 2001. ``Multidimensional Analysis of Roll Call Data
+ via Bayesian Simulation.'' \emph{Political Analysis.} 9: 227-241.
+
+ Valen E. Johnson and James H. Albert. 1999. ``Ordinal Data Modeling."
+ Springer: New York.
+
+ Andrew D. Martin, Kevin M. Quinn, and Daniel Pemstein. 2004.
+ \emph{Scythe Statistical Library 1.0.} \url{http://scythe.wustl.edu}.
+
+ Radford Neal. 2003. ``Slice Sampling'' (with discussion). \emph{Annals of
+ Statistics}, 31: 705-767.
+
+ Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2002.
+ \emph{Output Analysis and Diagnostics for MCMC (CODA)}.
+ \url{http://www-fis.iarc.fr/coda/}.
+
+ Douglas Rivers. 2004. ``Identification of Multidimensional Item-Response
+ Models." Stanford University, typescript.
+}
+
+
+\examples{
+ \dontrun{
+ ## Court example with ability (ideal point) and
+ ## item (case) constraints
+ data(SupremeCourt)
+ post1 <- MCMCirtKdRob(t(SupremeCourt), dimensions=1,
+ burnin=500, mcmc=5000, thin=1,
+ B0=.25, store.item=TRUE, store.ability=TRUE,
+ ability.constraints=list("Thomas"=list(1,"+"),
+ "Stevens"=list(1,-4)),
+ item.constraints=list("1"=list(2,"-")),
+ verbose=50)
+ plot(post1)
+ summary(post1)
+
+ ## Senate example with ability (ideal point) constraints
+ data(Senate)
+ namestring <- as.character(Senate$member)
+ namestring[78] <- "CHAFEE1"
+ namestring[79] <- "CHAFEE2"
+ namestring[59] <- "SMITH.NH"
+ namestring[74] <- "SMITH.OR"
+ rownames(Senate) <- namestring
+ post2 <- MCMCirtKdRob(Senate[,6:677], dimensions=1,
+ burnin=1000, mcmc=5000, thin=1,
+ ability.constraints=list("KENNEDY"=list(1,-4),
+ "HELMS"=list(1, 4), "ASHCROFT"=list(1,"+"),
+ "BOXER"=list(1,"-"), "KERRY"=list(1,"-"),
+ "HATCH"=list(1,"+")),
+ B0=0.1, store.ability=TRUE, store.item=FALSE,
+ verbose=5, k0=0.15, k1=0.15,
+ delta0.start=0.13, delta1.start=0.13)
+
+ plot(post2)
+ summary(post2)
+ }
+}
+
+\keyword{models}
+
+\seealso{\code{\link[coda]{plot.mcmc}},\code{\link[coda]{summary.mcmc}},
+\code{\link[MCMCpack]{MCMCirt1d}}, \code{\link[MCMCpack]{MCMCirtKd}}
+}
+
diff --git a/man/MCMClogit.Rd b/man/MCMClogit.Rd
index c0ab626..1636a37 100644
--- a/man/MCMClogit.Rd
+++ b/man/MCMClogit.Rd
@@ -13,7 +13,8 @@
\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, user.prior.density=NULL, logfun=TRUE, ...) }
+ b0 = 0, B0 = 0, user.prior.density=NULL, logfun=TRUE,
+ marginal.likelihood=c("none", "Laplace"), ...) }
\arguments{
\item{formula}{Model formula.}
@@ -84,6 +85,11 @@ MCMClogit(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
returns the natural log of a density function (TRUE) or a density
(FALSE).}
+ \item{marginal.likelihood}{How should the marginal likelihood be
+ calculated? Options are: \code{none} in which case the marginal
+ likelihood will not be calculated or \code{Laplace} in which case the
+ Laplace approximation (see Kass and Raftery, 1995) is used.}
+
\item{\ldots}{further arguments to be passed}
}
diff --git a/man/MCMCpoisson.Rd b/man/MCMCpoisson.Rd
index d2ef0aa..40154c9 100644
--- a/man/MCMCpoisson.Rd
+++ b/man/MCMCpoisson.Rd
@@ -13,7 +13,7 @@
\usage{
MCMCpoisson(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, marginal.likelihood = c("none", "Laplace"), ...) }
\arguments{
\item{formula}{Model formula.}
@@ -70,7 +70,12 @@ MCMCpoisson(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
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{marginal.likelihood}{How should the marginal likelihood be
+ calculated? Options are: \code{none} in which case the marginal
+ likelihood will not be calculated or \code{Laplace} in which case the
+ Laplace approximation (see Kass and Raftery, 1995) is used.}
+
\item{\ldots}{further arguments to be passed}
}
diff --git a/man/MCMCprobit.Rd b/man/MCMCprobit.Rd
index 3aa6b8a..6986d6a 100644
--- a/man/MCMCprobit.Rd
+++ b/man/MCMCprobit.Rd
@@ -13,7 +13,8 @@
\usage{
MCMCprobit(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
thin = 1, verbose = 0, seed = NA, beta.start = NA,
- b0 = 0, B0 = 0, bayes.resid = FALSE, ...) }
+ b0 = 0, B0 = 0, bayes.resid = FALSE,
+ marginal.likelihood=c("none", "Laplace"), ...) }
\arguments{
\item{formula}{Model formula.}
@@ -71,7 +72,12 @@ MCMCprobit(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
giving the observation numbers for which latent residuals should be
calculated and returned. TRUE will return draws of
latent residuals for all observations.}
-
+
+ \item{marginal.likelihood}{How should the marginal likelihood be
+ calculated? Options are: \code{none} in which case the marginal
+ likelihood will not be calculated or \code{Laplace} in which case the
+ Laplace approximation (see Kass and Raftery, 1995) is used.}
+
\item{...}{further arguments to be passed}
}
diff --git a/man/MCMCregress.Rd b/man/MCMCregress.Rd
index 1ccb195..a563f26 100644
--- a/man/MCMCregress.Rd
+++ b/man/MCMCregress.Rd
@@ -15,7 +15,8 @@
\usage{
MCMCregress(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
thin = 1, verbose = 0, seed = NA, beta.start = NA,
- b0 = 0, B0 = 0, c0 = 0.001, d0 = 0.001, ...) }
+ b0 = 0, B0 = 0, c0 = 0.001, d0 = 0.001,
+ marginal.likelihood = c("none", "Laplace", "Chib95"), ...) }
\arguments{
\item{formula}{Model formula.}
@@ -76,7 +77,13 @@ MCMCregress(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
disturbances). In constructing the inverse Gamma prior,
\eqn{d_0}{d0} acts like the sum of squared errors from the
\eqn{c_0}{c0} pseudo-observations.}
-
+
+ \item{marginal.likelihood}{How should the marginal likelihood be
+ calculated? Options are: \code{none} in which case the marginal
+ likelihood will not be calculated, \code{Laplace} in which case the
+ Laplace approximation (see Kass and Raftery, 1995) is used, and
+ \code{Chib95} in which case the method of Chib (1995) is used.}
+
\item{...}{further arguments to be passed}
}
@@ -110,7 +117,14 @@ MCMCregress(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
as the first block in the sampler.
}
-\references{
+ \references{
+ Siddhartha Chib. 1995. "Marginal Likelihood from the Gibbs Output."
+ \emph{Journal of the American Statistical Association}. 90:
+ 1313-1321.
+
+ Robert E. Kass and Adrian E. Raftery. 1995. "Bayes Factors."
+ \emph{Journal of the American Statistical Association}. 90: 773-795.
+
Andrew D. Martin, Kevin M. Quinn, and Daniel Pemstein. 2004.
\emph{Scythe Statistical Library 1.0.} \url{http://scythe.wustl.edu}.
diff --git a/man/MCMCtobit.Rd b/man/MCMCtobit.Rd
index 907f301..3d21a5e 100644
--- a/man/MCMCtobit.Rd
+++ b/man/MCMCtobit.Rd
@@ -143,7 +143,8 @@ MCMCtobit(formula, data = parent.frame(), below = 0, above = Inf,
\code{\link[survival]{survreg}},
\code{\link[MCMCpack]{MCMCregress}}}
-\examples{
+\examples{
+\dontrun{
library(survival)
example(tobin)
summary(tfit)
@@ -152,6 +153,7 @@ tfit.mcmc <- MCMCtobit(durable ~ age + quant, data=tobin, mcmc=30000,
plot(tfit.mcmc)
raftery.diag(tfit.mcmc)
summary(tfit.mcmc)
+}
}
\keyword{models}
diff --git a/man/MCbinomialbeta.Rd b/man/MCbinomialbeta.Rd
new file mode 100644
index 0000000..ad554fd
--- /dev/null
+++ b/man/MCbinomialbeta.Rd
@@ -0,0 +1,60 @@
+\name{MCbinomialbeta}
+\alias{MCbinomialbeta}
+\title{Monte Carlo Simulation from a Binomial Likelihood with a Beta Prior}
+
+\description{
+ This function generates a sample from the posterior distribution
+ of a binomial likelihood with a Beta prior.
+ }
+
+\usage{
+MCbinomialbeta(y, n, alpha=1, beta=1, mc=1000, ...)
+}
+
+\arguments{
+ \item{y}{The number of successes in the independent Bernoulli trials.}
+
+ \item{n}{The number of independent Bernoulli trials.}
+
+ \item{alpha}{Beta prior distribution alpha parameter.}
+
+ \item{beta}{Beta prior distribution beta parameter.}
+
+ \item{mc}{The number of Monte Carlo draws to make.}
+
+ \item{...}{further arguments to be passed}
+}
+
+\value{
+ An mcmc object that contains the posterior sample. This
+ object can be summarized by functions provided by the coda package.
+}
+
+\details{
+ \code{MCbinomialbeta} directly simulates from the posterior distribution.
+ This model is designed primarily for instructional use. \eqn{\pi}{pi} is
+ the probability of success for each independent Bernoulli trial. We assume
+ a conjugate Beta prior:
+ \deqn{\pi \sim \mathcal{B}eta(\alpha, \beta)}{pi ~ Beta(alpha, beta)}
+ \eqn{y} is the number of successes in \eqn{n} trials.
+ By default, a uniform prior is used.
+ }
+
+\examples{
+\dontrun{
+posterior <- MCbinomialbeta(3,12,mc=5000)
+summary(posterior)
+plot(posterior)
+grid <- seq(0,1,0.01)
+plot(grid, dbeta(grid, 1, 1), type="l", col="red", lwd=3, ylim=c(0,3.6),
+ xlab="pi", ylab="density")
+lines(density(posterior), col="blue", lwd=3)
+legend(.75, 3.6, c("prior", "posterior"), lwd=3, col=c("red", "blue"))
+}
+}
+
+\keyword{models}
+
+\seealso{\code{\link[coda]{plot.mcmc}},
+ \code{\link[coda]{summary.mcmc}}}
+
diff --git a/man/MCmultinomdirichlet.Rd b/man/MCmultinomdirichlet.Rd
new file mode 100644
index 0000000..295cda7
--- /dev/null
+++ b/man/MCmultinomdirichlet.Rd
@@ -0,0 +1,56 @@
+\name{MCmultinomdirichlet}
+\alias{MCmultinomdirichlet}
+\title{Monte Carlo Simulation from a Multinomial Likelihood with a
+ Dirichlet Prior}
+
+\description{
+ This function generates a sample from the posterior distribution
+ of a multinomial likelihood with a Dirichlet prior.
+ }
+
+\usage{
+MCmultinomdirichlet(y, alpha0, mc=1000, ...)
+}
+
+\arguments{
+ \item{y}{A vector of data (number of successes for each category).}
+
+ \item{alpha0}{The vector of parameters of the Dirichlet prior.}
+
+ \item{mc}{The number of Monte Carlo draws to make.}
+
+ \item{...}{further arguments to be passed}
+
+}
+
+\value{
+ An mcmc object that contains the posterior sample. This
+ object can be summarized by functions provided by the coda package.
+}
+
+\details{
+ \code{MCmultinomdirichlet} directly simulates from the posterior distribution.
+ This model is designed primarily for instructional use. \eqn{\pi}{pi}
+ is the parameter of interest of the multinomial distribution. It is of
+ dimension \eqn{(d \times 1)}{(d x 1)}. We assume
+ a conjugate Dirichlet prior:
+ \deqn{\pi \sim \mathcal{D}irichlet(\alpha_0)}{pi ~ Dirichlet(alpha0)}
+ \eqn{y} is a \eqn{(d \times 1)}{(d x 1)} vector of observed data.
+ }
+
+\examples{
+\dontrun{
+## Example from Gelman, et. al. (1995, p. 78)
+posterior <- MCmultinomdirichlet(c(727,583,137), c(1,1,1), mc=10000)
+bush.dukakis.diff <- posterior[,1] - posterior[,2]
+cat("Pr(Bush > Dukakis): ",
+ sum(bush.dukakis.diff > 0) / length(bush.dukakis.diff), "\n")
+hist(bush.dukakis.diff)
+}
+}
+
+\keyword{models}
+
+\seealso{\code{\link[coda]{plot.mcmc}},
+ \code{\link[coda]{summary.mcmc}}}
+
diff --git a/man/MCnormalnormal.Rd b/man/MCnormalnormal.Rd
new file mode 100644
index 0000000..86ca069
--- /dev/null
+++ b/man/MCnormalnormal.Rd
@@ -0,0 +1,62 @@
+\name{MCnormalnormal}
+\alias{MCnormalnormal}
+\title{Monte Carlo Simulation from a Normal Likelihood (with known variance) with a Normal Prior}
+
+\description{
+ This function generates a sample from the posterior distribution
+ of a Normal likelihood (with known variance) with a Normal prior.
+ }
+
+\usage{
+MCnormalnormal(y, sigma2, mu0, tau20, mc=1000, ...)
+}
+
+\arguments{
+ \item{y}{The data.}
+
+ \item{sigma2}{The known variance of y.}
+
+ \item{mu0}{The prior mean of mu.}
+
+ \item{tau20}{The prior variance of mu.}
+
+ \item{mc}{The number of Monte Carlo draws to make.}
+
+ \item{...}{further arguments to be passed}
+
+}
+
+\value{
+ An mcmc object that contains the posterior sample. This
+ object can be summarized by functions provided by the coda package.
+}
+
+\details{
+ \code{MCnormalnormal} directly simulates from the posterior distribution.
+ This model is designed primarily for instructional use. \eqn{\mu}{mu}
+ is the parameter of interest of the Normal distribution.
+ We assume
+ a conjugate normal prior:
+ \deqn{\mu \sim \mathcal{N}(\mu_0, \tau^2_0)}{mu ~ N(mu0, tau20)}
+ \eqn{y} is a vector of observed data.
+ }
+
+\examples{
+\dontrun{
+y <- c(2.65, 1.80, 2.29, 2.11, 2.27, 2.61, 2.49, 0.96, 1.72, 2.40)
+posterior <- MCMCpack:::MCnormalnormal(y, 1, 0, 1, 5000)
+summary(posterior)
+plot(posterior)
+grid <- seq(-3,3,0.01)
+plot(grid, dnorm(grid, 0, 1), type="l", col="red", lwd=3, ylim=c(0,1.4),
+ xlab="mu", ylab="density")
+lines(density(posterior), col="blue", lwd=3)
+legend(-3, 1.4, c("prior", "posterior"), lwd=3, col=c("red", "blue"))
+}
+}
+
+\keyword{models}
+
+\seealso{\code{\link[coda]{plot.mcmc}},
+ \code{\link[coda]{summary.mcmc}}}
+
diff --git a/man/MCpoissongamma.Rd b/man/MCpoissongamma.Rd
new file mode 100644
index 0000000..baf139e
--- /dev/null
+++ b/man/MCpoissongamma.Rd
@@ -0,0 +1,60 @@
+\name{MCpoissongamma}
+\alias{MCpoissongamma}
+\title{Monte Carlo Simulation from a Poisson Likelihood with a Gamma Prior}
+
+\description{
+ This function generates a sample from the posterior distribution
+ of a Poisson likelihood with a Gamma prior.
+ }
+
+\usage{
+MCpoissongamma(y, alpha, beta, mc=1000, ...)
+}
+
+\arguments{
+ \item{y}{A vector of counts (must be non-negative).}
+
+ \item{alpha}{Gamma prior distribution shape parameter.}
+
+ \item{beta}{Gamma prior distribution scale parameter.}
+
+ \item{mc}{The number of Monte Carlo draws to make.}
+
+ \item{...}{further arguments to be passed}
+
+}
+
+\value{
+ An mcmc object that contains the posterior sample. This
+ object can be summarized by functions provided by the coda package.
+}
+
+\details{
+ \code{MCpoissongamma} directly simulates from the posterior distribution.
+ This model is designed primarily for instructional use. \eqn{\lambda}{lambda}
+ is the parameter of interest of the Poisson distribution.
+ We assume
+ a conjugate Gamma prior:
+ \deqn{\lambda \sim \mathcal{G}amma(\alpha, \beta)}{lambda ~ Gamma(alpha, beta)}
+ \eqn{y} is a vector of counts.
+ }
+
+\examples{
+\dontrun{
+data(quine)
+posterior <- MCpoissongamma(quine$Days, 15, 1, 5000)
+summary(posterior)
+plot(posterior)
+grid <- seq(14,18,0.01)
+plot(grid, dgamma(grid, 15, 1), type="l", col="red", lwd=3, ylim=c(0,1.3),
+ xlab="lambda", ylab="density")
+lines(density(posterior), col="blue", lwd=3)
+legend(17, 1.3, c("prior", "posterior"), lwd=3, col=c("red", "blue"))
+}
+}
+
+\keyword{models}
+
+\seealso{\code{\link[coda]{plot.mcmc}},
+ \code{\link[coda]{summary.mcmc}}}
+
diff --git a/man/PostProbMod.Rd b/man/PostProbMod.Rd
new file mode 100644
index 0000000..7838d10
--- /dev/null
+++ b/man/PostProbMod.Rd
@@ -0,0 +1,66 @@
+\name{PostProbMod}
+\alias{PostProbMod}
+
+
+\title{Calculate Posterior Probability of Model}
+\description{
+This function takes an object of class \code{BayesFactor} and calculates
+the posterior probability that each model under study is correct given
+that one of the models under study is correct.
+
+}
+\usage{
+PostProbMod(BF, prior.probs=1)
+}
+
+\arguments{
+ \item{BF}{An object of class \code{BayesFactor}.}
+
+ \item{prior.probs}{The prior probabilities that each model is
+ correct. Can be either a scalar or array. Must be positive. If the
+ sum of the prior probabilities is not equal to 1 prior.probs will be
+ normalized so that it does sum to unity.}
+}
+
+\value{
+ An array holding the posterior probabilities that each model under
+ study is correct given that one of the models under study is correct.
+}
+
+\examples{
+\dontrun{
+data(birthwt)
+
+post1 <- MCMCregress(bwt~age+lwt+as.factor(race) + smoke + ht,
+ data=birthwt, b0=c(2700, 0, 0, -500, -500,
+ -500, -500),
+ B0=c(1e-6, .01, .01, 1.6e-5, 1.6e-5, 1.6e-5,
+ 1.6e-5), c0=10, d0=4500000,
+ marginal.likelihood="Chib95", mcmc=10000)
+
+post2 <- MCMCregress(bwt~age+lwt+as.factor(race) + smoke,
+ data=birthwt, b0=c(2700, 0, 0, -500, -500,
+ -500),
+ B0=c(1e-6, .01, .01, 1.6e-5, 1.6e-5, 1.6e-5),
+ c0=10, d0=4500000,
+ marginal.likelihood="Chib95", mcmc=10000)
+
+post3 <- MCMCregress(bwt~as.factor(race) + smoke + ht,
+ data=birthwt, b0=c(2700, -500, -500,
+ -500, -500),
+ B0=c(1e-6, 1.6e-5, 1.6e-5, 1.6e-5,
+ 1.6e-5), c0=10, d0=4500000,
+ marginal.likelihood="Chib95", mcmc=10000)
+
+BF <- BayesFactor(post1, post2, post3)
+mod.probs <- PostProbMod(BF)
+print(mod.probs)
+}
+}
+
+\concept{Bayes factor}
+\concept{model comparison}
+
+\seealso{\code{\link{MCMCregress}}}
+\keyword{models}
+
diff --git a/man/procrust.Rd b/man/procrust.Rd
new file mode 100644
index 0000000..f222a82
--- /dev/null
+++ b/man/procrust.Rd
@@ -0,0 +1,51 @@
+\name{procrustes}
+\alias{procrustes}
+
+
+\title{Procrustes Transformation}
+\description{
+This function performs a Procrustes transformation on a matrix \code{X} to
+minimize the squared distance between \code{X} and another matrix \code{Xstar}.
+}
+\usage{
+procrustes(X, Xstar, translation=FALSE, dilation=FALSE)
+}
+
+\arguments{
+ \item{X}{The matrix to be transformed.}
+
+ \item{Xstar}{The target matrix.}
+
+ \item{translation}{logical value indicating whether \code{X} should be
+ translated.}
+
+ \item{dilation}{logical value indicating whether \code{X} should be
+ dilated.}
+}
+
+\value{
+ A list containing: \code{X.new} the matrix that is the Procrustes
+ transformed version of \code{X}, \code{R} the rotation matrix,
+ \code{tt} the translation vector, and \code{s} the scale factor.
+}
+
+\details{
+ \code{R}, \code{tt}, and \code{s} are chosen so that
+
+ \deqn{s X R + 1 tt' \approx X^*}{s X R + 1 tt' approximately Xstar}.
+
+ \code{X.new} is given by:
+
+ \deqn{X_{new} = s X R + 1 tt'}{X.new = s X R + 1 tt'}.
+ }
+
+
+\references{
+ Borg and Groenen. 1997. \emph{Modern Multidimensional Scaling}. New
+ York: Springer. pp. 340-342.
+ }
+
+
+\seealso{\code{\link{MCMCirtKd}}}
+\keyword{manip}
+
diff --git a/man/xpnd.Rd b/man/xpnd.Rd
index dd7d997..f26a1fb 100644
--- a/man/xpnd.Rd
+++ b/man/xpnd.Rd
@@ -13,7 +13,8 @@
\arguments{
\item{x}{A list of elements to expand into symmetric matrix.}
- \item{nrow}{The number of rows (and columns) in the returned matrix.}
+ \item{nrow}{The number of rows (and columns) in the returned matrix.
+ Look into the details.}
}
\value{
@@ -25,13 +26,21 @@
covariance matrices. Note that R stores matrices in column major
order, and that the items in \code{x} will be recycled to fill the
matrix if need be.
+
+ The number of rows can be specified or automatically computed from
+ the number of elements in a given object via
+ \eqn{(-1 + \sqrt{(1 + 8 * length(x))}) / 2}.
}
\examples{
xpnd(c(1,2,3,4,4,5,6,7,8,9),4)
+ xpnd(c(1,2,3,4,4,5,6,7,8,9))
}
\keyword{manip}
+\concept{triangular}
+
\seealso{\code{\link{vech}}}
+
diff --git a/src/MCMCSVDreg.cc b/src/MCMCSVDreg.cc
new file mode 100644
index 0000000..bd77052
--- /dev/null
+++ b/src/MCMCSVDreg.cc
@@ -0,0 +1,235 @@
+// MCMCSVDreg.R samples from the posterior distribution of a Gaussian
+// linear regression model in which the X matrix has been decomposed
+// with an SVD. Useful for prediction when number of columns of X
+// is (possibly much) greater than the number of rows of X.
+//
+// See West, Mike. 2000. "Bayesian Regression in the 'Large p, Small n'
+// Paradigm." Duke ISDS Discussion Paper 2000-22.
+//
+// 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 9/9/2005
+//
+
+#include "matrix.h"
+#include "distributions.h"
+#include "stat.h"
+#include "la.h"
+#include "ide.h"
+#include "smath.h"
+#include "MCMCrng.h"
+
+#include <iostream>
+
+#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" {
+
+ // simulate from posterior distribution and return an mcmc by parameters
+ // matrix of the posterior sample
+ void MCMCSVDreg(double *sampledata, const int *samplerow,
+ const int *samplecol,
+ const double *Ydata, const int *Yrow, const int *Ycol,
+ const int *Ymiss,
+ const double *Adata, const int *Arow, const int *Acol,
+ const double *Ddata, const int *Drow, const int *Dcol,
+ const double *Fdata, const int *Frow, const int *Fcol,
+ const int *burnin, const int *mcmc,
+ const int *thin, const int *lecuyer,
+ const int *seedarray,
+ const int *lecuyerstream, const int *verbose,
+ const double *taustartdata, const int *taustartrow,
+ const int *taustartcol,
+ const double *g0data,
+ const int *g0row, const int *g0col,
+ const double *a0, const double *b0,
+ const double* c0, const double* d0,
+ const double* w0,
+ const int* betasamp) {
+
+
+ // pull together Matrix objects
+ Matrix <double> y = r2scythe(*Yrow, *Ycol, Ydata);
+ Matrix <double> A = r2scythe(*Arow, *Acol, Adata);
+ Matrix <double> D = r2scythe(*Drow, *Dcol, Ddata);
+ Matrix <double> F = r2scythe(*Frow, *Fcol, Fdata);
+ Matrix <double> g0 = r2scythe(*g0row, *g0col, g0data);
+
+ // define constants
+ const int tot_iter = *burnin + *mcmc; // total number of mcmc iterations
+ const int nstore = *mcmc / *thin; // number of draws to store
+ const int k = *Arow;
+ const int n = *Yrow;
+ Matrix <double> FtD = t(F) * D;
+ Matrix <double> DinvF = invpd(D) * F;
+ Matrix <double> Dg0 = D * g0;
+ double dsquared[n];
+ for (int i=0; i<n; ++i){
+ dsquared[i] = std::pow(D(i,i), 2);
+ }
+ int holder = 0;
+ for (int i=0; i<n; ++i){
+ if (Ymiss[i] == 1){
+ ++holder;
+ }
+ }
+ const int nYmiss = holder;
+
+ // storage matrices
+ Matrix<double> beta_store;
+ if (*betasamp == 1){
+ beta_store = Matrix<double>(k, nstore);
+ }
+ Matrix<double> gamma_store(n, nstore);
+ Matrix<double> Y_store(nYmiss, nstore);
+ Matrix<double> sigma2_store(1, nstore);
+ Matrix<double> tau2_store(n, nstore);
+
+ // initialize rng stream
+ rng *stream = MCMCpack_get_rng(*lecuyer, seedarray, *lecuyerstream);
+
+
+ // set starting values
+ double tau2[n];
+ for (int i=0; i<n; ++i){
+ tau2[i] = taustartdata[i];
+ }
+ Matrix <double> gamma(n, 1);
+ double sigma2 = 1;
+ Matrix <double> beta;
+ if (*betasamp == 1){
+ beta = Matrix<double>(k, 1);
+ }
+
+
+
+ /////////////////// Gibbs sampler ///////////////////
+ int count = 0;
+ for (int iter = 0; iter < tot_iter; ++iter) {
+
+ // sample [sigma2 | Y, A, D, F, tau2]
+ Matrix <double> Fy = F * y;
+ Fy = Fy - Dg0;
+ double q = 0.0;
+ for (int i=0; i<n; ++i){
+ q += (std::pow(Fy[i], 2) / (1 + tau2[i]));
+ }
+ sigma2 = stream->rigamma((*a0+n)*0.5, (*b0 + q)*0.5);
+
+
+ // sample [gamma | Y, A, D, F, sigma2, tau2]
+ // w0[i] is assumed to be the prior prob that gamma[i] == 0
+ Matrix <double> gammahat = DinvF * y;
+ for (int i=0; i<n; ++i){
+ double mstar = (g0[i] + tau2[i]*gammahat[i]) / (1 + tau2[i]);
+ double vstar = (sigma2 * tau2[i]) / (dsquared[i] * (tau2[i] + 1));
+ double wstar = 1.0 - (1.0 - w0[i]) /
+ (1.0 - w0[i] + w0[i] * dnorm(0.0, mstar, std::sqrt(vstar)));
+ if (stream->runif() < wstar){
+ gamma[i] = 0.0;
+ }
+ else {
+ gamma[i] = stream->rnorm(mstar, std::sqrt(vstar));
+ }
+ }
+
+
+ // sample [tau2 | Y, A, D, F, gamma, sigma2]
+ for (int i=0; i<n; ++i){
+ double gamg2 = std::pow((gamma[i] - g0[i]), 2);
+ tau2[i] = stream->rigamma((1+c0[i])*0.5,
+ (gamg2 * (dsquared[i] / sigma2) + d0[i])
+ *0.5);
+ }
+
+
+ // sample [y[miss] | Y[obs], A, D, F, gamma, sigma2, tau2]
+ Matrix <double> eta = FtD * gamma;
+ for (int i=0; i<n; ++i){
+ if (Ymiss[i] == 1){
+ y[i] = stream->rnorm(eta[i], std::sqrt(sigma2));
+ }
+ }
+
+
+ // optional sample [beta | Y, A, D, F, gamma, sigma2, tau2]
+ if (*betasamp == 1){
+ beta = A * gamma;
+ }
+
+
+
+
+
+ // store draws in storage matrix (or matrices)
+ if (iter >= *burnin && (iter % *thin == 0)) {
+ sigma2_store[count] = sigma2;
+ int Ymiss_count = 0;
+ for (int i = 0; i<n; ++i){
+ gamma_store(i, count) = gamma[i];
+ tau2_store(i, count) = tau2[i];
+ if (Ymiss[i] == 1){
+ Y_store(Ymiss_count, count) = y[i];
+ ++Ymiss_count;
+ }
+ }
+ if (*betasamp == 1){
+ for (int i=0; i<k; ++i){
+ beta_store(i, count) = beta[i];
+ }
+ }
+ ++count;
+ }
+
+
+ // print output to stdout
+ if(*verbose > 0 && iter % *verbose == 0) {
+ Rprintf("\n\nMCMCSVDreg iteration %i of %i \n",
+ (iter+1), tot_iter);
+ Rprintf("gamma = \n");
+ for (int j=0; j<n; ++j)
+ Rprintf("%10.5f\n", gamma[j]);
+ Rprintf("tau2 = \n");
+ for (int j=0; j<n; ++j)
+ Rprintf("%10.5f\n", tau2[j]);
+ Rprintf("sigma2 = %10.5f\n", sigma2);
+ }
+
+ R_CheckUserInterrupt(); // allow user interrupts
+
+ } // end MCMC loop
+
+ delete stream; // clean up random number stream
+
+ // load draws into sample array
+ Matrix <double> storagematrix = cbind(t(Y_store), t(gamma_store));
+ storagematrix = cbind(storagematrix, t(tau2_store));
+ storagematrix = cbind(storagematrix, t(sigma2_store));
+ if (*betasamp == 1){
+ storagematrix = cbind(storagematrix, t(beta_store));
+ }
+
+ const int size = *samplerow * *samplecol;
+ for(int i = 0; i < size; ++i)
+ sampledata[i] = storagematrix[i];
+
+ } // end MCMCregress
+} // end extern "C"
diff --git a/src/MCMCfcds.cc b/src/MCMCfcds.cc
index 55ad9d8..ccb2349 100644
--- a/src/MCMCfcds.cc
+++ b/src/MCMCfcds.cc
@@ -311,7 +311,7 @@ namespace SCYTHE {
const Matrix<double> AB0ab0 = AB0 * ab0;
// perform update
- const Matrix<double> Ttheta_star = t(cbind(-1.0*ones<double>(J,1),theta)); // only needed for option 2
+ //const Matrix<double> Ttheta_star = t(cbind(-1.0*ones<double>(J,1),theta)); // only needed for option 2
const Matrix<double> tpt(2,2);
for (int i=0; i<J; ++i){
const double theta_i = theta[i];
@@ -357,8 +357,8 @@ namespace SCYTHE {
const double T0t0 = T0*t0;
const Matrix<double> alpha = eta(_, 0);
const Matrix<double> beta = eta(_, 1);
- const Matrix<double> tbeta = t(beta); // only neede for option 2
- const Matrix<double> talpha = t(alpha); // only needed for option 2
+ //const Matrix<double> tbeta = t(beta); // only neede for option 2
+ //const Matrix<double> talpha = t(alpha); // only needed for option 2
// calculate the posterior variance outside the justice specific loop
double theta_post_var = T0;
diff --git a/src/MCMCirt1d.cc b/src/MCMCirt1d.cc
index c1848bc..bd84d14 100644
--- a/src/MCMCirt1d.cc
+++ b/src/MCMCirt1d.cc
@@ -4,6 +4,7 @@
// ADM and KQ 1/15/2003
// ADM 7/28/2004 [updated to new Scythe version]
// completely rewritten and optimized for the 1-d case 8/2/2004 KQ
+// storage changed to save memory KQ 1/27/2006
#include "matrix.h"
#include "distributions.h"
@@ -60,7 +61,8 @@ extern "C"{
const double* thetaineqdata,
const int* thetaineqrow,
const int* thetaineqcol,
- const int* store
+ const int* storei,
+ const int* storea
) {
using namespace SCYTHE; //Added by Matthew S. Fasman on 11/04/2004
@@ -93,9 +95,15 @@ extern "C"{
// storage matrices (row major order)
- Matrix<double> theta_store = Matrix<double>(nsamp, J);
- Matrix<double> eta_store = Matrix<double>(nsamp, K*2);
-
+ Matrix<double> theta_store;
+ Matrix<double> eta_store;
+ if (*storea == 1){
+ theta_store = Matrix<double>(nsamp, J);
+ }
+ if (*storei == 1){
+ eta_store = Matrix<double>(nsamp, K*2);
+ }
+
// starting values
Matrix<double> eta = cbind(alpha, beta);
Matrix<double> Z = Matrix<double>(J,K);
@@ -127,15 +135,19 @@ extern "C"{
if ((iter >= burnin[0]) && ((iter % thin[0]==0))) {
// store ideal points
- for (int l=0; l<J; ++l)
- theta_store(count, l) = theta[l];
-
+ if (*storea == 1){
+ for (int l=0; l<J; ++l)
+ theta_store(count, l) = theta[l];
+ }
+
// store bill parameters
- for (int l=0; l<K*2; ++l)
- eta_store(count, l) = eta[l];
+ if (*storei == 1){
+ for (int l=0; l<K*2; ++l)
+ eta_store(count, l) = eta[l];
+ }
count++;
}
-
+
R_CheckUserInterrupt(); // allow user interrupts
} // end Gibbs loop
@@ -144,9 +156,12 @@ extern "C"{
// return output
Matrix<double> output;
- if(*store == 0) {
+ if(*storei == 0 && *storea == 1) {
output = theta_store;
}
+ else if (*storei == 1 && *storea == 0){
+ output = eta_store;
+ }
else {
output = cbind(theta_store, eta_store);
}
diff --git a/src/MCMCirtKdRob.cc b/src/MCMCirtKdRob.cc
new file mode 100644
index 0000000..89fe595
--- /dev/null
+++ b/src/MCMCirtKdRob.cc
@@ -0,0 +1,1062 @@
+// MCMCirtKdRob.cc is C++ code to estimate a robust K-dimensional
+// item response model
+//
+// 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
+//
+// 2/18/2005 KQ
+//
+
+#include <iostream>
+
+#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;
+
+
+/* Equal probability sampling; without-replacement case */
+// pulled from R-2.0.1/src/main/random.c lines 352-364
+// slightly modified by KQ 2/21/2005
+// x: n array of original indices running from 0 to (n-1)
+// y: k array of samled indices
+// k: length of y (must be <= n)
+// n: length of x (must be >= k)
+static void SampleNoReplace(const int k, int n, int *y, int *x,
+ rng* stream){
+
+ for (int i = 0; i < n; i++)
+ x[i] = i;
+ for (int i = 0; i < k; i++) {
+ int j = static_cast<int>(n * stream->runif());
+ y[i] = x[j];
+ x[j] = x[--n];
+ }
+}
+
+
+
+// full conditional distribution for a single element of Lambda
+// this single element is Lambda(rowindex, colindex)
+static double Lambda_logfcd(const double& lam_ij,
+ const Matrix<int>& X,
+ const Matrix<double>& Lambda,
+ const Matrix<double>& theta,
+ const double& delta0,
+ const double& delta1,
+ const Matrix<double>& Lambda_prior_mean,
+ const Matrix<double>& Lambda_prior_prec,
+ const Matrix<double>& Lambda_ineq,
+ const Matrix<double>& theta_ineq,
+ const double& k0,
+ const double& k1,
+ const double& c0,
+ const double& d0,
+ const double& c1,
+ const double& d1,
+ const int& rowindex,
+ const int& colindex){
+
+
+ const int D = Lambda.cols();
+
+ // check to see if inequality constraint is satisfied and
+ // evaluate prior
+ /*
+ double logprior = 0.0;
+ if (Lambda_ineq(rowindex,colindex) * lam_ij < 0){
+ return log(0.0);
+ }
+ if (Lambda_prior_prec(rowindex,colindex) != 0){
+ logprior += lndnorm(lam_ij,
+ Lambda_prior_mean(rowindex,colindex),
+ sqrt(1.0 / Lambda_prior_prec(rowindex,colindex)));
+ }
+ */
+ // prior is uniform on hypersphere with radius 10
+ if (Lambda_ineq(rowindex,colindex) * lam_ij < 0){
+ return log(0.0);
+ }
+ double absqdist = 0.0;
+ for (int i=0; i<D; ++i){
+ if (i==colindex){
+ absqdist += ::pow(lam_ij, 2);
+ }
+ else{
+ absqdist += ::pow(Lambda(rowindex,i), 2);
+ }
+ }
+ if (absqdist > 100.0){
+ return log(0.0);
+ }
+ const double logprior = 0.0;
+
+
+ // likelihood
+ double loglike = 0.0;
+ const int N = X.rows();
+ for (int i=0; i<N; ++i){
+ if (X(i,rowindex) != -999){
+ double eta = 0.0;
+ for (int j=0; j<D; ++j){
+ if (j==colindex){
+ eta += theta(i,j) * lam_ij;
+ }
+ else{
+ eta += theta(i,j) * Lambda(rowindex,j);
+ }
+ }
+ const double p = delta0 + (1 - delta0 - delta1) *
+ (1.0 / (1.0 + exp(-1*eta)));
+ loglike += X(i,rowindex) * log(p) + (1.0 - X(i,rowindex)) * log(1.0 - p);
+ }
+ }
+
+ return (loglike + logprior);
+}
+
+
+
+// full conditional for a single element of theta
+// this single element is theta(rowindex, colindex)
+static double theta_logfcd(const double& t_ij,
+ const Matrix<int>& X,
+ const Matrix<double>& Lambda,
+ const Matrix<double>& theta,
+ const double& delta0,
+ const double& delta1,
+ const Matrix<double>& Lambda_prior_mean,
+ const Matrix<double>& Lambda_prior_prec,
+ const Matrix<double>& Lambda_ineq,
+ const Matrix<double>& theta_ineq,
+ const double& k0,
+ const double& k1,
+ const double& c0,
+ const double& d0,
+ const double& c1,
+ const double& d1,
+ const int& rowindex,
+ const int& colindex){
+
+ const int D = Lambda.cols();
+ // evaluate prior
+ /*
+ if (theta_ineq(rowindex,colindex-1) * t_ij < 0){
+ return log(0.0);
+ }
+ const double logprior = lndnorm(t_ij, 0.0, 1.0);
+ */
+ // prior is uniform on unit circle
+ if (theta_ineq(rowindex,colindex-1) * t_ij < 0){
+ return log(0.0);
+ }
+ double tsqdist = 0.0;
+ for (int i=0; i<(D-1); ++i){
+ if (i==(colindex-1)){
+ tsqdist += ::pow(t_ij, 2);
+ }
+ else{
+ tsqdist += ::pow(theta(rowindex,(i-1)), 2);
+ }
+ }
+ if (tsqdist > 1.0){
+ return log(0.0);
+ }
+ const double logprior = 1.0;
+
+
+ // likelihood
+ double loglike = 0.0;
+ const int K = X.cols();
+ for (int i=0; i<K; ++i){
+ if (X(rowindex,i) != -999){
+ double eta = 0.0;
+ for (int j=0; j<D; ++j){
+ if (j==colindex){
+ eta += t_ij * Lambda(i,j);
+ }
+ else{
+ eta += theta(rowindex,j) * Lambda(i,j);
+ }
+ }
+ const double p = delta0 + (1 - delta0 - delta1) *
+ (1.0 / (1.0 + exp(-1*eta)));
+ loglike += X(rowindex,i) * log(p) + (1.0 - X(rowindex,i)) * log(1.0 - p);
+ }
+ }
+
+ return (loglike + logprior);
+}
+
+
+
+// full conditional for delta0
+static double delta0_logfcd(const double& delta0,
+ const Matrix<int>& X,
+ const Matrix<double>& Lambda,
+ const Matrix<double>& theta,
+ const double& junk,
+ const double& delta1,
+ const Matrix<double>& Lambda_prior_mean,
+ const Matrix<double>& Lambda_prior_prec,
+ const Matrix<double>& Lambda_ineq,
+ const Matrix<double>& theta_ineq,
+ const double& k0,
+ const double& k1,
+ const double& c0,
+ const double& d0,
+ const double& c1,
+ const double& d1,
+ const int& rowindex,
+ const int& colindex){
+
+ // evaluate prior
+ if (delta0 >=k0 || delta0 <=0){
+ return log(0.0);
+ }
+ const double logprior = lndbeta1(delta0 * (1.0/k0), c0, d0);
+
+ // likelihood
+ double loglike = 0.0;
+ const int D = Lambda.cols();
+ const int N = X.rows();
+ const int K = X.cols();
+ for (int i=0; i<N; ++i){
+ for (int k=0; k<K; ++k){
+ if (X(i,k) != -999){
+ double eta = 0.0;
+ for (int j=0; j<D; ++j){
+ eta += theta(i,j) * Lambda(k,j);
+ }
+ const double p = delta0 + (1 - delta0 - delta1) *
+ (1.0 / (1.0 + exp(-1*eta)));
+ loglike += X(i,k) * log(p) + (1.0 - X(i,k)) * log(1.0 - p);
+ }
+ }
+ }
+
+ return (loglike + logprior);
+}
+
+
+
+// full conditional for delta1
+static double delta1_logfcd(const double& delta1,
+ const Matrix<int>& X,
+ const Matrix<double>& Lambda,
+ const Matrix<double>& theta,
+ const double& delta0,
+ const double& junk,
+ const Matrix<double>& Lambda_prior_mean,
+ const Matrix<double>& Lambda_prior_prec,
+ const Matrix<double>& Lambda_ineq,
+ const Matrix<double>& theta_ineq,
+ const double& k0,
+ const double& k1,
+ const double& c0,
+ const double& d0,
+ const double& c1,
+ const double& d1,
+ const int& rowindex,
+ const int& colindex){
+
+ // evaluate prior
+ if (delta1 >=k1 || delta1 <=0){
+ return log(0.0);
+ }
+ const double logprior = lndbeta1(delta1 * (1.0/k1), c1, d1);
+
+ // likelihood
+ double loglike = 0.0;
+ const int D = Lambda.cols();
+ const int N = X.rows();
+ const int K = X.cols();
+ for (int i=0; i<N; ++i){
+ for (int k=0; k<K; ++k){
+ if (X(i,k) != -999){
+ double eta = 0.0;
+ for (int j=0; j<D; ++j){
+ eta += theta(i,j) * Lambda(k,j);
+ }
+ const double p = delta0 + (1 - delta0 - delta1) *
+ (1.0 / (1.0 + exp(-1*eta)));
+ loglike += X(i,k) * log(p) + (1.0 - X(i,k)) * log(1.0 - p);
+ }
+ }
+ }
+
+ return (loglike + logprior);
+}
+
+
+// Radford Neal's (2000) doubling procedure coded for a logdensity
+static void doubling(double (*logfun)(const double&,
+ const Matrix<int>&,
+ const Matrix<double>&,
+ const Matrix<double>&,
+ const double&,
+ const double&,
+ const Matrix<double>&,
+ const Matrix<double>&,
+ const Matrix<double>&,
+ const Matrix<double>&,
+ const double&,
+ const double&,
+ const double&,
+ const double&,
+ const double&,
+ const double&,
+ const int&,
+ const int&),
+ const Matrix<int>& X,
+ const Matrix<double>& Lambda,
+ const Matrix<double>& theta,
+ const double& delta0,
+ const double& delta1,
+ const Matrix<double>& Lambda_prior_mean,
+ const Matrix<double>& Lambda_prior_prec,
+ const Matrix<double>& Lambda_ineq,
+ const Matrix<double>& theta_ineq,
+ const double& k0,
+ const double& k1,
+ const double& c0,
+ const double& d0,
+ const double& c1,
+ const double& d1,
+ const int& rowindex,
+ const int& colindex,
+ const double& z,
+ const double& w, const int& p,
+ rng* stream, double& L, double& R,
+ const int& param){
+
+ const double U = stream->runif();
+ double x0;
+ if (param==0){ // Lambda
+ x0 = Lambda(rowindex, colindex);
+ }
+ else if (param==1){ // theta
+ x0 = theta(rowindex, colindex);
+ }
+ else if (param==2){ // delta0
+ x0 = delta0;
+ }
+ else if (param==3){ // delta1
+ x0 = delta1;
+ }
+ else {
+ Rprintf("\nERROR: param not in {0,1,2,3} in doubling().\n");
+ exit(1);
+ }
+
+ L = x0 - w*U;
+ R = L + w;
+ int K = p;
+ while (K > 0 &&
+ (z < logfun(L, X, Lambda, theta, delta0, delta1, Lambda_prior_mean,
+ Lambda_prior_prec, Lambda_ineq, theta_ineq,
+ k0, k1, c0, d0, c1, d1, rowindex, colindex) |
+ z < logfun(R, X, Lambda, theta, delta0, delta1, Lambda_prior_mean,
+ Lambda_prior_prec, Lambda_ineq, theta_ineq,
+ k0, k1, c0, d0, c1, d1, rowindex, colindex))){
+ double V = stream->runif();
+ if (V < 0.5){
+ L = L - (R - L);
+ }
+ else {
+ R = R + (R - L);
+ }
+ --K;
+ }
+
+}
+
+
+// Radford Neal's (2000) stepping out procedure coded for a logdensity
+static void StepOut(double (*logfun)(const double&,
+ const Matrix<int>&,
+ const Matrix<double>&,
+ const Matrix<double>&,
+ const double&,
+ const double&,
+ const Matrix<double>&,
+ const Matrix<double>&,
+ const Matrix<double>&,
+ const Matrix<double>&,
+ const double&,
+ const double&,
+ const double&,
+ const double&,
+ const double&,
+ const double&,
+ const int&,
+ const int&),
+ const Matrix<int>& X,
+ const Matrix<double>& Lambda,
+ const Matrix<double>& theta,
+ const double& delta0,
+ const double& delta1,
+ const Matrix<double>& Lambda_prior_mean,
+ const Matrix<double>& Lambda_prior_prec,
+ const Matrix<double>& Lambda_ineq,
+ const Matrix<double>& theta_ineq,
+ const double& k0,
+ const double& k1,
+ const double& c0,
+ const double& d0,
+ const double& c1,
+ const double& d1,
+ const int& rowindex,
+ const int& colindex,
+ const double& z,
+ const double& w, const int& m,
+ rng* stream, double& L, double& R,
+ const int& param){
+
+ const double U = stream->runif();
+ double x0;
+ if (param==0){ // Lambda
+ x0 = Lambda(rowindex, colindex);
+ }
+ else if (param==1){ // theta
+ x0 = theta(rowindex, colindex);
+ }
+ else if (param==2){ // delta0
+ x0 = delta0;
+ }
+ else if (param==3){ // delta1
+ x0 = delta1;
+ }
+ else {
+ Rprintf("\nERROR: param not in {0,1,2,3} in StepOut().\n");
+ exit(1);
+ }
+
+
+ L = x0 - w*U;
+ R = L + w;
+ const double V = stream->runif();
+ int J = static_cast<int>(m*V);
+ int K = (m-1) - J;
+ while (J > 0 &&
+ (z < logfun(L, X, Lambda, theta, delta0, delta1, Lambda_prior_mean,
+ Lambda_prior_prec, Lambda_ineq, theta_ineq,
+ k0, k1, c0, d0, c1, d1, rowindex, colindex))){
+ L = L - w;
+ J = J - 1;
+ }
+ while (K > 0 &&
+ (z < logfun(R, X, Lambda, theta, delta0, delta1, Lambda_prior_mean,
+ Lambda_prior_prec, Lambda_ineq, theta_ineq,
+ k0, k1, c0, d0, c1, d1, rowindex, colindex))){
+ R = R + w;
+ K = K - 1;
+ }
+
+}
+
+
+
+
+// Radford Neal's (2000) Accept procedure coded for a logdensity
+static const bool Accept(double (*logfun)(const double&,
+ const Matrix<int>&,
+ const Matrix<double>&,
+ const Matrix<double>&,
+ const double&,
+ const double&,
+ const Matrix<double>&,
+ const Matrix<double>&,
+ const Matrix<double>&,
+ const Matrix<double>&,
+ const double&,
+ const double&,
+ const double&,
+ const double&,
+ const double&,
+ const double&,
+ const int&,
+ const int&),
+ const Matrix<int>& X,
+ const Matrix<double>& Lambda,
+ const Matrix<double>& theta,
+ const double& delta0,
+ const double& delta1,
+ const Matrix<double>& Lambda_prior_mean,
+ const Matrix<double>& Lambda_prior_prec,
+ const Matrix<double>& Lambda_ineq,
+ const Matrix<double>& theta_ineq,
+ const double& k0,
+ const double& k1,
+ const double& c0,
+ const double& d0,
+ const double& c1,
+ const double& d1,
+ const int& rowindex,
+ const int& colindex,
+ const double& z,
+ const double& w, const double& x0,
+ const double& x1, const double& L, const double& R){
+
+ double Lhat = L;
+ double Rhat = R;
+ bool D = false;
+
+ while ((Rhat - Lhat ) > 1.1 * w){
+ double M = (Lhat + Rhat) / 2.0;
+ if ( (x0 < M && x1 >= M) || (x0 >= M && x1 < M)){
+ D = true;
+ }
+ if (x1 < M){
+ Rhat = M;
+ }
+ else {
+ Lhat = M;
+ }
+
+ if (D && z >= logfun(Lhat, X, Lambda, theta, delta0, delta1,
+ Lambda_prior_mean, Lambda_prior_prec,
+ Lambda_ineq, theta_ineq, k0, k1, c0, d0,
+ c1, d1, rowindex, colindex) &&
+ z >= logfun(Rhat, X, Lambda, theta, delta0, delta1, Lambda_prior_mean,
+ Lambda_prior_prec, Lambda_ineq, theta_ineq,
+ k0, k1, c0, d0, c1, d1, rowindex, colindex)){
+ return(false);
+ }
+ }
+ return(true);
+}
+
+
+// Radford Neal's (2000) shrinkage procedure coded for a log density
+// assumes the doubling procedure has been used to find L and R
+static double shrinkageDoubling(double (*logfun)(const double&,
+ const Matrix<int>&,
+ const Matrix<double>&,
+ const Matrix<double>&,
+ const double&,
+ const double&,
+ const Matrix<double>&,
+ const Matrix<double>&,
+ const Matrix<double>&,
+ const Matrix<double>&,
+ const double&,
+ const double&,
+ const double&,
+ const double&,
+ const double&,
+ const double&,
+ const int&,
+ const int&),
+ const Matrix<int>& X,
+ const Matrix<double>& Lambda,
+ const Matrix<double>& theta,
+ const double& delta0,
+ const double& delta1,
+ const Matrix<double>& Lambda_prior_mean,
+ const Matrix<double>& Lambda_prior_prec,
+ const Matrix<double>& Lambda_ineq,
+ const Matrix<double>& theta_ineq,
+ const double& k0,
+ const double& k1,
+ const double& c0,
+ const double& d0,
+ const double& c1,
+ const double& d1,
+ const int& rowindex,
+ const int& colindex,
+ const double& z, const double& w,
+ rng* stream, const double& L, const double& R,
+ const int& param){
+
+ double Lbar = L;
+ double Rbar = R;
+ double x0;
+ if (param==0){ // Lambda
+ x0 = Lambda(rowindex, colindex);
+ }
+ else if (param==1){ // theta
+ x0 = theta(rowindex, colindex);
+ }
+ else if (param==2){ // delta0
+ x0 = delta0;
+ }
+ else if (param==3){ // delta1
+ x0 = delta1;
+ }
+ else {
+ Rprintf("\nERROR: param not in {0,1,2,3} in shrinkageDoubling().\n");
+ exit(1);
+ }
+
+ for (;;){
+ const double U = stream->runif();
+ const double x1 = Lbar + U*(Rbar - Lbar);
+ if (z < logfun(x1, X, Lambda, theta, delta0, delta1, Lambda_prior_mean,
+ Lambda_prior_prec, Lambda_ineq, theta_ineq,
+ k0, k1, c0, d0, c1, d1, rowindex, colindex) &&
+ Accept(logfun, X, Lambda, theta, delta0, delta1, Lambda_prior_mean,
+ Lambda_prior_prec, Lambda_ineq, theta_ineq,
+ k0, k1, c0, d0, c1, d1, rowindex, colindex, z, w,
+ x0, x1, L, R)){
+ return(x1);
+ }
+ if (x1 < x0){
+ Lbar = x1;
+ }
+ else {
+ Rbar = x1;
+ }
+ } // end infinite loop
+}
+
+
+
+// Radford Neal's (2000) shrinkage procedure coded for a log density
+// assumes the stepping out procedure has been used to find L and R
+static double shrinkageStep(double (*logfun)(const double&,
+ const Matrix<int>&,
+ const Matrix<double>&,
+ const Matrix<double>&,
+ const double&,
+ const double&,
+ const Matrix<double>&,
+ const Matrix<double>&,
+ const Matrix<double>&,
+ const Matrix<double>&,
+ const double&,
+ const double&,
+ const double&,
+ const double&,
+ const double&,
+ const double&,
+ const int&,
+ const int&),
+ const Matrix<int>& X,
+ const Matrix<double>& Lambda,
+ const Matrix<double>& theta,
+ const double& delta0,
+ const double& delta1,
+ const Matrix<double>& Lambda_prior_mean,
+ const Matrix<double>& Lambda_prior_prec,
+ const Matrix<double>& Lambda_ineq,
+ const Matrix<double>& theta_ineq,
+ const double& k0,
+ const double& k1,
+ const double& c0,
+ const double& d0,
+ const double& c1,
+ const double& d1,
+ const int& rowindex,
+ const int& colindex,
+ const double& z, const double& w,
+ rng* stream, const double& L, const double& R,
+ const int& param){
+
+ double Lbar = L;
+ double Rbar = R;
+ double x0;
+ if (param==0){ // Lambda
+ x0 = Lambda(rowindex, colindex);
+ }
+ else if (param==1){ // theta
+ x0 = theta(rowindex, colindex);
+ }
+ else if (param==2){ // delta0
+ x0 = delta0;
+ }
+ else if (param==3){ // delta1
+ x0 = delta1;
+ }
+ else {
+ Rprintf("\nERROR: param not in {0,1,2,3} in shrinkageDoubling().\n");
+ exit(1);
+ }
+
+ for (;;){
+ const double U = stream->runif();
+ const double x1 = Lbar + U*(Rbar - Lbar);
+ if (z < logfun(x1, X, Lambda, theta, delta0, delta1, Lambda_prior_mean,
+ Lambda_prior_prec, Lambda_ineq, theta_ineq,
+ k0, k1, c0, d0, c1, d1, rowindex, colindex) ){
+ return(x1);
+ }
+ if (x1 < x0){
+ Lbar = x1;
+ }
+ else {
+ Rbar = x1;
+ }
+ } // end infinite loop
+}
+
+
+
+
+extern "C"{
+
+// function called by R to fit model
+ void
+ irtKdRobpost (double* sampledata, const int* samplerow,
+ const int* samplecol,
+ const int* Xdata, const int* Xrow, const int* Xcol,
+ const int* burnin, const int* mcmc, const int* thin,
+ const int *lecuyer, const int *seedarray,
+ const int *lecuyerstream, const int* verbose,
+ const int* method_step,
+ const double* theta_w, const int* theta_p,
+ const double* lambda_w, const int* lambda_p,
+ const double* delta0_w, const int* delta0_p,
+ const double* delta1_w, const int* delta1_p,
+ const double * delta0start, const double* delta1start,
+ const double* Lamstartdata, const int* Lamstartrow,
+ const int* Lamstartcol,
+ const double* thetstartdata, const int* thetstartrow,
+ const int* thetstartcol,
+ const double* Lameqdata, const int* Lameqrow,
+ const int* Lameqcol,
+ const double* Lamineqdata, const int* Lamineqrow,
+ const int* Lamineqcol,
+ const double* theteqdata, const int* theteqrow,
+ const int* theteqcol,
+ const double* thetineqdata, const int* thetineqrow,
+ const int* thetineqcol,
+ const double* Lampmeandata, const int* Lampmeanrow,
+ const int* Lampmeancol,
+ const double* Lampprecdata, const int* Lampprecrow,
+ const int* Lamppreccol,
+ const double* k0, const double* k1,
+ const double* c0, const double* c1,
+ const double* d0, const double* d1,
+ const int* storeitem,
+ const int* storeability
+ ) {
+
+ // put together matrices
+ const Matrix<int> X = r2scythe(*Xrow, *Xcol, Xdata);
+ Matrix<double> Lambda = r2scythe(*Lamstartrow, *Lamstartcol, Lamstartdata);
+ Matrix<double> theta = r2scythe(*thetstartrow, *thetstartcol, thetstartdata);
+ const Matrix<double> Lambda_eq = r2scythe(*Lameqrow, *Lameqcol, Lameqdata);
+ const Matrix<double> Lambda_ineq = r2scythe(*Lamineqrow, *Lamineqcol,
+ Lamineqdata);
+ const Matrix<double> theta_eq = r2scythe(*theteqrow, *theteqcol,
+ theteqdata);
+ const Matrix<double> theta_ineq = r2scythe(*thetineqrow, *thetineqcol,
+ thetineqdata);
+ const Matrix<double> Lambda_prior_mean = r2scythe(*Lampmeanrow,
+ *Lampmeancol,
+ Lampmeandata);
+ const Matrix<double> Lambda_prior_prec = r2scythe(*Lampprecrow,
+ *Lamppreccol,
+ Lampprecdata);
+
+
+ // initialize rng stream
+ rng *stream = MCMCpack_get_rng(*lecuyer, seedarray, *lecuyerstream);
+
+ // constants
+ const int K = X.cols(); // number of items
+ const int N = X.rows(); // number of subjects
+ const int D = Lambda.cols(); // number of dimensions + 1
+ const int tot_iter = *burnin + *mcmc;
+ const int nsamp = *mcmc / *thin;
+ // const Matrix<double> Lambda_free_indic = Matrix<double>(K, D);
+ //for (int i=0; i<(K*D); ++i){
+ // if (Lambda_eq[i] == -999) Lambda_free_indic[i] = 1.0;
+ //}
+
+ // starting values
+ // Matrix<double> theta = Matrix<double>(N, D-1);
+ theta = cbind(-1*ones<double>(N,1), theta);
+ double delta0 = *delta0start;
+ double delta1 = *delta1start;
+
+ // index arrays (used for the random order of the sampling)
+ int K_array[K];
+ int N_array[N];
+ int D_array[D];
+ int Dm1_array[D-1];
+ int K_inds_perm[K];
+ int N_inds_perm[N];
+ int D_inds_perm[D];
+ int Dm1_inds_perm[D-1];
+
+
+ // storage matrices (row major order)
+ Matrix<double> Lambda_store;
+ if (storeitem[0]==1){
+ Lambda_store = Matrix<double>(nsamp, K*D);
+ }
+ Matrix<double> theta_store;
+ if (*storeability==1){
+ theta_store = Matrix<double>(nsamp, N*D);
+ }
+ Matrix<double> delta0_store = Matrix<double>(nsamp, 1);
+ Matrix<double> delta1_store = Matrix<double>(nsamp, 1);
+
+ ///////////////////
+ // Slice Sampler //
+ ///////////////////
+ int count = 0;
+ for (int iter=0; iter < tot_iter; ++iter){
+
+
+ double L, R, w, funval, z;
+ int p;
+
+ // sample theta
+ int param = 1;
+ SampleNoReplace(N, N, N_inds_perm, N_array, stream);
+ for (int ii=0; ii<N; ++ii){
+ int i = N_inds_perm[ii];
+ SampleNoReplace(D-1, D-1, Dm1_inds_perm, Dm1_array, stream);
+ for (int jj=0; jj<(D-1); ++jj){
+ int j = Dm1_inds_perm[jj]+1;
+ if (theta_eq(i,j-1) == -999){
+ w = *theta_w;
+ p = *theta_p;
+ L = -1.0;
+ R = 1.0;
+ funval = theta_logfcd(theta(i,j), X, Lambda,
+ theta, delta0,
+ delta1, Lambda_prior_mean,
+ Lambda_prior_prec,
+ Lambda_ineq, theta_ineq,
+ *k0, *k1, *c0, *d0,
+ *c1, *d1, i, j);
+ z = funval - stream->rexp(1.0);
+ if (*method_step == 1){
+ StepOut(&theta_logfcd, X, Lambda, theta, delta0, delta1,
+ Lambda_prior_mean, Lambda_prior_prec, Lambda_ineq,
+ theta_ineq, *k0, *k1, *c0, *d0, *c1, *d1, i, j, z,
+ w, p, stream, L, R, param);
+ theta(i,j) = shrinkageStep(&theta_logfcd, X, Lambda, theta,
+ delta0, delta1, Lambda_prior_mean,
+ Lambda_prior_prec, Lambda_ineq,
+ theta_ineq, *k0, *k1,
+ *c0, *d0, *c1,
+ *d1, i, j, z, w, stream,
+ L, R, param);
+ }
+ else{
+ doubling(&theta_logfcd, X, Lambda, theta, delta0, delta1,
+ Lambda_prior_mean, Lambda_prior_prec, Lambda_ineq,
+ theta_ineq, *k0, *k1, *c0, *d0, *c1, *d1, i, j, z,
+ w, p, stream, L, R, param);
+ theta(i,j) = shrinkageDoubling(&theta_logfcd, X, Lambda, theta,
+ delta0, delta1, Lambda_prior_mean,
+ Lambda_prior_prec, Lambda_ineq,
+ theta_ineq, *k0, *k1,
+ *c0, *d0, *c1,
+ *d1, i, j, z, w, stream,
+ L, R, param);
+ }
+ }
+ }
+ }
+
+
+ // sample Lambda
+ param = 0;
+ SampleNoReplace(K, K, K_inds_perm, K_array, stream);
+ for (int ii=0; ii<K; ++ii){
+ int i = K_inds_perm[ii];
+ SampleNoReplace(D, D, D_inds_perm, D_array, stream);
+ for (int jj=0; jj<D; ++jj){
+ int j = D_inds_perm[jj];
+ if (Lambda_eq(i,j) == -999){
+ w = *lambda_w;
+ p = *lambda_p;
+ L = -1.0;
+ R = 1.0;
+ funval = Lambda_logfcd(Lambda(i,j), X, Lambda,
+ theta, delta0,
+ delta1, Lambda_prior_mean,
+ Lambda_prior_prec,
+ Lambda_ineq, theta_ineq, *k0, *k1, *c0, *d0,
+ *c1, *d1, i, j);
+ z = funval - stream->rexp(1.0);
+ if (*method_step == 1){
+ StepOut(&Lambda_logfcd, X, Lambda, theta, delta0, delta1,
+ Lambda_prior_mean, Lambda_prior_prec, Lambda_ineq,
+ theta_ineq, *k0, *k1, *c0, *d0, *c1, *d1, i, j, z,
+ w, p, stream, L, R, param);
+ Lambda(i,j) = shrinkageStep(&Lambda_logfcd, X, Lambda, theta,
+ delta0, delta1, Lambda_prior_mean,
+ Lambda_prior_prec, Lambda_ineq,
+ theta_ineq, *k0, *k1, *c0, *d0,
+ *c1, *d1, i, j, z, w, stream,
+ L, R, param);
+ }
+ else{
+ doubling(&Lambda_logfcd, X, Lambda, theta, delta0, delta1,
+ Lambda_prior_mean, Lambda_prior_prec, Lambda_ineq,
+ theta_ineq, *k0, *k1, *c0, *d0, *c1, *d1, i, j, z,
+ w, p, stream, L, R, param);
+ Lambda(i,j) = shrinkageDoubling(&Lambda_logfcd, X, Lambda, theta,
+ delta0, delta1, Lambda_prior_mean,
+ Lambda_prior_prec, Lambda_ineq,
+ theta_ineq, *k0, *k1, *c0, *d0,
+ *c1, *d1, i, j, z, w, stream,
+ L, R, param);
+ }
+ }
+ }
+ }
+
+ // sample delta0
+ param = 2;
+ w = *delta0_w;
+ p = *delta0_p;
+ L = -1.0;
+ R = 1.0;
+ funval = delta0_logfcd(delta0, X, Lambda,
+ theta, delta0,
+ delta1, Lambda_prior_mean,
+ Lambda_prior_prec,
+ Lambda_ineq, theta_ineq,
+ *k0, *k1, *c0, *d0,
+ *c1, *d1, 0, 0);
+ z = funval - stream->rexp(1.0);
+ if (*method_step == 1){
+ StepOut(&delta0_logfcd, X, Lambda, theta, delta0, delta1,
+ Lambda_prior_mean, Lambda_prior_prec, Lambda_ineq,
+ theta_ineq, *k0, *k1, *c0, *d0, *c1, *d1, 0, 0, z,
+ w, p, stream, L, R, param);
+ delta0 = shrinkageStep(&delta0_logfcd, X, Lambda, theta,
+ delta0, delta1, Lambda_prior_mean,
+ Lambda_prior_prec, Lambda_ineq, theta_ineq,
+ *k0, *k1, *c0, *d0, *c1, *d1, 0, 0,
+ z, w, stream, L, R, param);
+ }
+ else{
+ doubling(&delta0_logfcd, X, Lambda, theta, delta0, delta1,
+ Lambda_prior_mean, Lambda_prior_prec, Lambda_ineq,
+ theta_ineq, *k0, *k1, *c0, *d0, *c1, *d1, 0, 0, z,
+ w, p, stream, L, R, param);
+ delta0 = shrinkageDoubling(&delta0_logfcd, X, Lambda, theta,
+ delta0, delta1, Lambda_prior_mean,
+ Lambda_prior_prec, Lambda_ineq, theta_ineq,
+ *k0, *k1, *c0, *d0, *c1, *d1, 0, 0,
+ z, w, stream, L, R, param);
+ }
+ // sample delta1
+ param = 3;
+ w = *delta1_w;
+ p = *delta1_p;
+ L = -1.0;
+ R = 1.0;
+ funval = delta1_logfcd(delta1, X, Lambda,
+ theta, delta0,
+ delta1, Lambda_prior_mean,
+ Lambda_prior_prec,
+ Lambda_ineq, theta_ineq, *k0, *k1, *c0, *d0,
+ *c1, *d1, 0, 0);
+ z = funval - stream->rexp(1.0);
+ if (*method_step == 1){
+ StepOut(&delta1_logfcd, X, Lambda, theta, delta0, delta1,
+ Lambda_prior_mean, Lambda_prior_prec, Lambda_ineq,
+ theta_ineq, *k0, *k1, *c0, *d0, *c1, *d1, 0, 0, z,
+ w, p, stream, L, R, param);
+ delta1 = shrinkageStep(&delta1_logfcd, X, Lambda, theta,
+ delta0, delta1, Lambda_prior_mean,
+ Lambda_prior_prec, Lambda_ineq, theta_ineq,
+ *k0, *k1, *c0, *d0, *c1, *d1, 0, 0,
+ z, w, stream, L, R, param);
+ }
+ else{
+ doubling(&delta1_logfcd, X, Lambda, theta, delta0, delta1,
+ Lambda_prior_mean, Lambda_prior_prec, Lambda_ineq,
+ theta_ineq, *k0, *k1, *c0, *d0, *c1, *d1, 0, 0, z,
+ w, p, stream, L, R, param);
+ delta1 = shrinkageDoubling(&delta1_logfcd, X, Lambda, theta,
+ delta0, delta1, Lambda_prior_mean,
+ Lambda_prior_prec, Lambda_ineq, theta_ineq,
+ *k0, *k1, *c0, *d0, *c1, *d1, 0, 0,
+ z, w, stream, L, R, param);
+ }
+
+ // print results to screen
+ if (verbose[0] > 0 && iter % verbose[0] == 0){
+ Rprintf("\n\nMCMCirtKdRob iteration %i of %i \n", (iter+1), tot_iter);
+ }
+
+ // store results
+ if ((iter >= burnin[0]) && ((iter % thin[0]==0))) {
+
+ // store Lambda
+ if (storeitem[0]==1){
+ Matrix<double> Lambda_store_vec = reshape(Lambda,1,K*D);
+ for (int l=0; l<K*D; ++l)
+ Lambda_store(count, l) = Lambda_store_vec[l];
+ }
+
+ // store theta
+ if (storeability[0]==1){
+ Matrix<double> theta_store_vec = reshape(theta, 1, N*D);
+ for (int l=0; l<N*D; ++l)
+ theta_store(count, l) = theta_store_vec[l];
+ }
+
+ // store delta0 and delta1
+ delta0_store[count] = delta0;
+ delta1_store[count] = delta1;
+
+ count++;
+ }
+
+ // allow user interrupts
+ void R_CheckUserInterrupt(void);
+ } // end MCMC loop
+
+ delete stream; // clean up random number stream
+
+ // return output
+ Matrix<double> output = delta0_store;
+ output = cbind(output, delta1_store);
+ if (*storeitem == 1){
+ output = cbind(output, Lambda_store);
+ }
+ if(*storeability == 1) {
+ output = cbind(output, theta_store);
+ }
+
+ int size = *samplerow * *samplecol;
+ for (int i=0; i<size; ++i)
+ sampledata[i] = output[i];
+
+
+}
+
+}
+
+
+
+
+
+
+
+
diff --git a/src/MCMCregress.cc b/src/MCMCregress.cc
index 1c277fa..26efd10 100644
--- a/src/MCMCregress.cc
+++ b/src/MCMCregress.cc
@@ -43,6 +43,15 @@
using namespace SCYTHE;
using namespace std;
+
+static double digamma(const double& theta, const double& a, const double& b){
+ double logf = a * log(b) - lngammafn(a) + -(a+1) * log(theta) + -b/theta;
+
+ return exp(logf);
+ //pow(b, a) / gammafn(a) * pow(theta, -(a+1)) * exp(-b/theta);
+}
+
+
extern "C" {
// simulate from posterior density and return an mcmc by parameters
@@ -57,7 +66,8 @@ extern "C" {
const int *betastartcol, const double *b0data,
const int *b0row, const int *b0col,
const double *B0data, const int *B0row,
- const int *B0col, const double *c0, const double *d0) {
+ const int *B0col, const double *c0, const double *d0,
+ double* logmarglikeholder, const int* chib) {
// pull together Matrix objects
Matrix <double> Y = r2scythe(*Yrow, *Ycol, Ydata);
@@ -114,13 +124,73 @@ extern "C" {
} // end MCMC loop
- delete stream; // clean up random number stream
+ if (*chib == 1){
+ // marginal likelihood calculation stuff starts here
+ const double sigma2star = meanc(t(sigmamatrix))[0];
+ double sigma2fcdsum = 0.0;
+ // second set of Gibbs scans
+ for (int iter = 0; iter < tot_iter; ++iter) {
+ double sigma2 = NormIGregress_sigma2_draw (X, Y, beta, *c0,
+ *d0, stream);
+ beta = NormNormregress_beta_draw (XpX, XpY, b0, B0, sigma2,
+ stream);
+ const Matrix <double> e = gaxpy(X, (-1*beta), Y);
+ const Matrix <double> SSE = crossprod (e);
+ const double c_post = (*c0 + X.rows ()) * 0.5;
+ const double d_post = (*d0 + SSE[0]) * 0.5;
+ sigma2fcdsum += digamma(sigma2star, c_post, d_post);
+
+ // print output to stdout
+ if(*verbose > 0 && iter % *verbose == 0) {
+ Rprintf("\n\nMCMCregress (reduced) iteration %i of %i \n",
+ (iter+1), tot_iter);
+ }
+
+ R_CheckUserInterrupt(); // allow user interrupts
+
+ } // end MCMC loop
+ double sigma2fcdmean = sigma2fcdsum / static_cast<double>(tot_iter);
+
+
+ const Matrix<double> betastar = t(meanc(t(betamatrix)));
+ const double sig2_inv = 1.0 / sigma2star;
+ const Matrix <double> sig_beta = invpd (B0 + XpX * sig2_inv);
+ const Matrix <double> betahat = sig_beta *
+ gaxpy(B0, b0, XpY*sig2_inv);
+ const double logbetafcd = lndmvn(betastar, betahat, sig_beta);
+
+
+
+ // calculate loglikelihood at (betastar, sigma2star)
+ double sigmastar = sqrt(sigma2star);
+ Matrix<double> eta = X * betastar;
+ double loglike = 0.0;
+ for (int i=0; i<X.rows(); ++i){
+ loglike += lndnorm(Y[i], eta[i], sigmastar);
+ }
+
+ // calculate log prior ordinate
+ double logprior = log(digamma(sigma2star, *c0/2.0, *d0/2.0)) +
+ lndmvn(betastar, b0, invpd(B0));
+
+
+ // put pieces together and print the marginal likelihood
+ double logmarglike = loglike + logprior - logbetafcd -
+ log(sigma2fcdmean);
+
+
+ logmarglikeholder[0] = logmarglike;
+
+ }
+ delete stream; // clean up random number stream
+
// load draws into sample array
Matrix <double> storeagematrix = cbind (t (betamatrix), t (sigmamatrix));
const int size = *samplerow * *samplecol;
for(int i = 0; i < size; ++i)
- sampledata[i] = storeagematrix[i];
-
+ sampledata[i] = storeagematrix[i];
+
+
} // end MCMCregress
} // end extern "C"
diff --git a/src/Makevars b/src/Makevars
index 3bff6ac..174a0f5 100644
--- a/src/Makevars
+++ b/src/Makevars
@@ -1,2 +1 @@
PKG_CXXFLAGS = -DSCYTHE_COMPILE_DIRECT -DSCYTHE_NO_RANGE -DHAVE_TRUNC
-
diff --git a/src/Makevars.in b/src/Makevars.in
index f71ce14..5c12b4a 100644
--- a/src/Makevars.in
+++ b/src/Makevars.in
@@ -1,2 +1 @@
PKG_CXXFLAGS = -DSCYTHE_COMPILE_DIRECT -DSCYTHE_NO_RANGE @MV_HAVE_IEEEFP_H@ @MV_HAVE_TRUNC@
-
--
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