[r-cran-mcmcpack] 03/90: Imported Upstream version 0.6-1
Andreas Tille
tille at debian.org
Fri Dec 16 09:07:31 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 fa4a6c0b50b1b1712fb361b2bc79f6ff7f34001e
Author: Andreas Tille <tille at debian.org>
Date: Fri Dec 16 08:07:01 2016 +0100
Imported Upstream version 0.6-1
---
DESCRIPTION | 8 +-
HISTORY | 30 +-
INDEX | 31 +-
LICENSE | 2 +-
NAMESPACE | 4 +-
R/MCMCfactanal.R | 2 +-
R/MCMCirtKd.R | 2 +-
R/MCMCmixfactanal.R | 30 +-
R/MCMCmnl.R | 292 ++++++++
R/MCMCordfactanal.R | 18 +-
R/MCMCtobit.R | 65 ++
R/automate.R | 3 +-
R/distn.R | 7 +
R/hidden.R | 11 +-
R/zzz.R | 2 +-
README | 6 +-
config.status | 712 -------------------
data/Nethvote.rda | Bin 0 -> 167791 bytes
man/MCMCdynamicEI.Rd | 2 +-
man/MCMCfactanal.Rd | 8 +-
man/MCMChierEI.Rd | 2 +-
man/MCMCirt1d.Rd | 5 +-
man/MCMCirtKd.Rd | 19 +-
man/MCMClogit.Rd | 2 +-
man/MCMCmixfactanal.Rd | 33 +-
man/MCMCmnl.Rd | 266 +++++++
man/MCMCoprobit.Rd | 2 +-
man/MCMCordfactanal.Rd | 8 +-
man/MCMCpanel.Rd | 2 +-
man/MCMCpoisson.Rd | 2 +-
man/MCMCprobit.Rd | 2 +-
man/MCMCtobit.Rd | 155 ++++
man/Nethvote.Rd | 56 ++
man/PErisk.Rd | 6 +-
man/choicevar.Rd | 25 +
man/tomog.Rd | 2 +-
man/wishart.Rd | 4 +-
src/MCMCdistn.cc | 466 ++++++++++++
src/MCMChierEI.cc | 16 +-
src/MCMCirt1d.cc | 2 +
src/MCMCmixfactanal.cc | 76 +-
src/MCMCmnlMH.cc | 171 +++++
src/MCMCmnlslice.cc | 335 +++++++++
src/MCMCoprobit.cc | 380 +++++-----
src/MCMCordfactanal.cc | 34 +-
src/MCMCregress.cc | 22 +-
src/MCMCrng.cc | 4 +-
src/MCMCrng.h | 4 +-
src/{MCMCregress.cc => MCMCtobit.cc} | 112 +--
src/Makevars | 2 +-
src/Makevars.in | 2 +-
src/optimize.cc | 1284 +++++++++++++++++-----------------
src/rng.cc | 194 ++---
src/smath.cc | 80 +--
src/smath.h | 24 +-
55 files changed, 3138 insertions(+), 1896 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index dfc8e96..46fc5b5 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,11 +1,11 @@
Package: MCMCpack
-Version: 0.5-2
-Date: 2004-8-17
+Version: 0.6-1
+Date: 2005-2-18
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 (>= 1.9.0), coda (>= 0.8-1), MASS
+Depends: R (>= 2.0.0), coda (>= 0.8-1), MASS
Description: This package contains functions for posterior
simulation for a number of statistical models. All simulation
is done in compiled C++ written in the Scythe Statistical
@@ -17,4 +17,4 @@ Description: This package contains functions for posterior
sampling algorithm, and tools for visualization.
License: GPL version 2 or newer
URL: http://mcmcpack.wustl.edu
-Packaged: Tue Aug 17 13:21:51 2004; adm
+Packaged: Sun Feb 20 22:22:48 2005; hornik
diff --git a/HISTORY b/HISTORY
index 7abe47c..fa62d93 100644
--- a/HISTORY
+++ b/HISTORY
@@ -2,6 +2,32 @@
// Changes and Bug Fixes
//
+0.5-2 to 0.6-1
+ * added a tobit model for a linear model with censoring in MCMCtobit()
+ * added multinomial logit model in MCMCmnl()
+ * added vote data from the Netherlands to illustrate MCMCmnl()
+ * added choicevar() function to specify choice-specific variables in
+ multinomial choice models
+ * fixed some data-handling issues in MCMCmixfactanal() and
+ MCMCordfactanal() [thanks to Ben Goodrich for isolating these and
+ providing patches]
+ * fixed the Metropolis-Hastings step for the Cowles algorithm for
+ cutpoints in MCMCoprobit(), MCMCordfactanal(), and MCMCmixfactanal()
+ [thanks to Alexander Raach for isolating the problem and providing
+ a patch]
+ * removed gcc specific compilation flags from Makevars.in (per the
+ request of Brian Ripley and Kurt Hornik)
+ * added informative message for templates created in auto.Scythe.call()
+ * modified vector.tune() hidden function
+ * fixed Neal's shrinkage procedure in MCMChierEI.cc
+ * a number of editorial fixes, updating for new calendar year, etc.
+ * removed functions acosh, asinh, atanh, and expm1 from smath.h and
+ smath.cc so cross-compilation will work
+
+//
+// Old Changes and Bug Fixes
+//
+
MCMCpack 0.5-1 was a major revision of MCMCpack. The entire package was
been essentially re-written using the new development environment (documented
in the MCMCpack specification) and the new Scythe Statistical Library 1.0.
@@ -43,10 +69,6 @@ exhaustive.
* MCMCirt1d() has a new interface with new types of constraints--
sampling for this model is also now much faster.
-//
-// Old Changes and Bug Fixes
-//
-
0.4-8 to 0.4-9
* Fixed a minor Scythe issue to fix error found by gcc 3.4.
diff --git a/INDEX b/INDEX
index 7aa12a1..a6eb58c 100644
--- a/INDEX
+++ b/INDEX
@@ -1,39 +1,46 @@
Dirichlet The Dirichlet Distribution
InvGamma The Inverse Gamma Distribution
InvWishart The Inverse Wishart Distribution
-MCMCdynamicEI Markov chain Monte Carlo for Quinn's Dynamic
+MCMCdynamicEI Markov Chain Monte Carlo for Quinn's Dynamic
Ecological Inference Model
-MCMCfactanal Markov chain Monte Carlo for Normal Theory
+MCMCfactanal Markov Chain Monte Carlo for Normal Theory
Factor Analysis Model
-MCMChierEI Markov chain Monte Carlo for Wakefield's
+MCMChierEI Markov Chain Monte Carlo for Wakefield's
Hierarchial Ecological Inference Model
-MCMCirt1d Markov chain Monte Carlo for One Dimensional
+MCMCirt1d Markov Chain Monte Carlo for One Dimensional
Item Response Theory Model
-MCMCirtKd Markov chain Monte Carlo for K-Dimensional
+MCMCirtKd Markov Chain Monte Carlo for K-Dimensional
Item Response Theory Model
-MCMClogit Markov chain Monte Carlo for Logistic
+MCMClogit Markov Chain Monte Carlo for Logistic
Regression
MCMCmetrop1R Metropolis Sampling from User-Written R
function
-MCMCmixfactanal Markov chain Monte Carlo for Mixed Data Factor
+MCMCmixfactanal Markov Chain Monte Carlo for Mixed Data Factor
Analysis Model
-MCMCoprobit Markov chain Monte Carlo for Ordered Probit
+MCMCmnl Markov Chain Monte Carlo for Multinomial
+ Logistic Regression
+MCMCoprobit Markov Chain Monte Carlo for Ordered Probit
Regression
-MCMCordfactanal Markov chain Monte Carlo for Ordinal Data
+MCMCordfactanal Markov Chain Monte Carlo for Ordinal Data
Factor Analysis Model
-MCMCpanel Markov chain Monte Carlo for the General
+MCMCpanel Markov Chain Monte Carlo for the General
Linear Panel Model
-MCMCpoisson Markov chain Monte Carlo for Poisson
+MCMCpoisson Markov Chain Monte Carlo for Poisson
Regression
-MCMCprobit Markov chain Monte Carlo for Probit Regression
+MCMCprobit Markov Chain Monte Carlo for Probit Regression
MCMCregress Markov Chain Monte Carlo for Gaussian Linear
Regression
+MCMCtobit Markov Chain Monte Carlo for Gaussian Linear
+ Regression with a Censored Dependent Variable
+Nethvote Dutch Voting Behavior in 1989
NoncenHypergeom The Noncentral Hypergeometric Distribution
PErisk Political Economic Risk Data from 62 Countries
in 1987
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
read.Scythe Read a Matrix from a File written by Scythe
tomogplot Tomography Plot
diff --git a/LICENSE b/LICENSE
index 9ff4eee..ac61a0c 100644
--- a/LICENSE
+++ b/LICENSE
@@ -1,5 +1,5 @@
Markov chain Monte Carlo Package (MCMCpack)
-Copyright (C) 2003, 2004 Andrew D. Martin and Kevin M. Quinn
+Copyright (C) 2003-2005 Andrew D. Martin and Kevin M. Quinn
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
diff --git a/NAMESPACE b/NAMESPACE
index 32d40fe..ac1ab79 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -3,6 +3,7 @@ import(coda)
import(MASS)
export(
+ choicevar,
ddirichlet,
dinvgamma,
diwish,
@@ -17,12 +18,14 @@ export(
MCMClogit,
MCMCmetrop1R,
MCMCmixfactanal,
+ MCMCmnl,
MCMCoprobit,
MCMCordfactanal,
MCMCpanel,
MCMCpoisson,
MCMCprobit,
MCMCregress,
+ MCMCtobit,
rdirichlet,
read.Scythe,
rinvgamma,
@@ -34,4 +37,3 @@ export(
write.Scythe,
xpnd
)
-
diff --git a/R/MCMCfactanal.R b/R/MCMCfactanal.R
index 0c1d9e2..77c568f 100644
--- a/R/MCMCfactanal.R
+++ b/R/MCMCfactanal.R
@@ -26,7 +26,7 @@
"MCMCfactanal" <-
function(x, factors, lambda.constraints=list(),
- data=parent.environment(), burnin = 1000, mcmc = 20000,
+ data=parent.frame(), burnin = 1000, mcmc = 20000,
thin=1, verbose = FALSE, seed = NA,
lambda.start = NA, psi.start = NA,
l0=0, L0=0, a0=0.001, b0=0.001,
diff --git a/R/MCMCirtKd.R b/R/MCMCirtKd.R
index 5d90394..4c9115d 100644
--- a/R/MCMCirtKd.R
+++ b/R/MCMCirtKd.R
@@ -21,7 +21,7 @@
store.item=FALSE, store.ability=TRUE,
drop.constantvars=TRUE, ... ) {
- datamatrix <- t(as.matrix(datamatrix))
+ datamatrix <- as.matrix(datamatrix)
post <- MCMCordfactanal(x=datamatrix, factors=dimensions,
lambda.constraints=item.constraints,
diff --git a/R/MCMCmixfactanal.R b/R/MCMCmixfactanal.R
index f745651..efa88b7 100644
--- a/R/MCMCmixfactanal.R
+++ b/R/MCMCmixfactanal.R
@@ -27,7 +27,7 @@
"MCMCmixfactanal" <-
function(x, factors, lambda.constraints=list(),
- data=parent.environment(), burnin = 1000, mcmc = 20000,
+ data=parent.frame(), burnin = 1000, mcmc = 20000,
thin=1, tune=NA, verbose = FALSE, seed = NA,
lambda.start = NA, psi.start=NA,
l0=0, L0=0, a0=0.001, b0=0.001,
@@ -46,6 +46,7 @@
mf$store.lambda <- mf$store.scores <- mf$std.var <- mf$... <- NULL
mf$drop.unused.levels <- TRUE
mf[[1]] <- as.name("model.frame")
+ mf$na.action <- 'na.pass'
mf <- eval(mf, sys.frame(sys.parent()))
attributes(mt)$intercept <- 0
Xterm.length <- length(attr(mt, "variables"))
@@ -60,11 +61,7 @@
ncat[i] <- -999
X[is.na(X[,i]), i] <- -999
}
- else if (is.ordered(X[,i])){
- ncat[i] <- nlevels(X[,i])
- X[,i] <- as.integer(X[,i])
- X[is.na(X[,i]), i] <- -999
- }
+ else if (is.ordered(X[, i])) {
ncat[i] <- nlevels(X[, i])
temp <- as.integer(X[,i])
temp <- ifelse(is.na(X[,i]) | (X[,i] == "<NA>"), -999, temp)
X[, i] <- temp
}
else {
stop("Manifest variable ", dimnames(X)[[2]][i],
" neither ordered factor nor numeric variable.\n")
@@ -226,6 +223,7 @@
gamma[2,i] <- 0
gamma[3:ncat[i],i] <- (polr.out$zeta[2:(ncat[i]-1)] -
polr.out$zeta[1])*.588
+
gamma[ncat[i]+1,i] <- 300
}
}
@@ -252,7 +250,9 @@
sample <- matrix(data=0, mcmc/thin, K*(factors+1)+(factors+1)*N +
length(gamma)+K)
}
-
+
+ accepts <- matrix(0, K, 1)
+
# Call the C++ code to do the real work
posterior <- .C("mixfactanalpost",
samdata = as.double(sample),
@@ -301,14 +301,20 @@
b0col = as.integer(ncol(b0)),
storelambda = as.integer(store.lambda),
storescores = as.integer(store.scores),
- accepts = as.integer(0),
+ accepts = as.integer(accepts),
+ acceptsrow = as.integer(nrow(accepts)),
+ acceptscol = as.integer(ncol(accepts)),
PACKAGE="MCMCpack"
)
- cat(" overall acceptance rate = ",
- posterior$accepts / ((posterior$burnin+posterior$mcmc)*n.ord.ge3), "\n")
-
-
+ accepts <- matrix(posterior$accepts, posterior$acceptsrow,
+ posterior$acceptscol, byrow=TRUE)
+ rownames(accepts) <- X.names
+ colnames(accepts) <- ""
+ cat("\n\nAcceptance rates:\n")
+ print(t(accepts) / (posterior$burnin+posterior$mcmc), digits=2,
+ width=6)
+
# put together matrix and build MCMC object to return
sample <- matrix(posterior$samdata, posterior$samrow, posterior$samcol,
byrow=TRUE)
diff --git a/R/MCMCmnl.R b/R/MCMCmnl.R
new file mode 100644
index 0000000..879da08
--- /dev/null
+++ b/R/MCMCmnl.R
@@ -0,0 +1,292 @@
+## MCMCmnl.R samples from the posterior density of a multinomial
+## logit model using Metropolis-Hastings.
+##
+## KQ 12/22/2004
+
+## parse formula and return a list that contains the model response
+## matrix as element one, the model matrix as element two,
+## the column names of X as element three, the rownames of
+## X and y as element four, and the number of choices in the largest
+## choice set in element five
+"parse.formula.mnl" <- function(formula, data, baseline=NULL,
+ intercept=TRUE){
+
+ ## 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 <- mf$baseline <- NULL
+ mf$drop.unused.levels <- TRUE
+ mf[[1]] <- as.name("model.frame")
+
+ mf <- eval(mf, sys.frame(sys.parent()))
+ if (!intercept){
+ attributes(mt)$intercept <- 0
+ }
+
+
+ ## deal with Y
+ Y <- as.matrix(model.response(mf, "numeric")) # Y matrix
+ if (ncol(Y)==1){
+ Y <- factor(Y)
+ number.choices <- length(unique(Y))
+ choice.names <- sort(unique(Y))
+ Ymat <- matrix(NA, length(Y), number.choices)
+ colnames(Ymat) <- choice.names
+ for (i in 1:(number.choices)){
+ Ymat[,i] <- as.numeric(Y==choice.names[i])
+ }
+ }
+ else{
+ ## this block will allow for nonconstant choice sets
+ number.choices <- ncol(Y)
+ Ytemp <- Y
+ Ytemp[Y== -999] <- NA
+ if ( min(unique(array(Y)) %in% c(-999,0,1))==0 ||
+ min(apply(Ytemp, 1, sum, na.rm=TRUE) == rep(1, nrow(Y)))==0){
+ stop("Y is a matrix but it is not composed of 0/1/-999 values\n and/or rows do not sum to 1\n")
+ }
+ Ymat <- Y
+ choice.names <- colnames(Y)
+ }
+ colnames(Ymat) <- choice.names
+ rownames(Ymat) <- 1:nrow(Ymat)
+
+ #Y.long <- matrix(t(Ymat), length(Ymat), 1)
+ #colnames(Y.long) <- "Y"
+ #rownames(Y.long) <- rep(1:nrow(Ymat), rep(length(choice.names), nrow(Ymat)))
+ #rownames(Y.long) <- paste(rownames(Y.long), choice.names, sep=".")
+ #group.id <- rep(1:nrow(Ymat), rep(ncol(Ymat), nrow(Ymat)))
+
+ ## deal with X
+ ## null model support
+ X <- if (!is.empty.model(mt)) model.matrix(mt, mf, contrasts)
+ X <- as.matrix(X) # X matrix
+ xvars <- dimnames(X)[[2]] # X variable names
+ xobs <- dimnames(X)[[1]] # X observation names
+
+ if (is.null(baseline)){
+ baseline <- choice.names[1]
+ }
+ if (! baseline %in% choice.names){
+ stop("'baseline' not consistent with observed choice levels in y\n")
+ }
+
+ ## deal with choice specific covariates
+ choicevar.indic <- rep(FALSE, length(xvars)) # indicators for choice
+ # specific variables
+ choicevar.indic[grep("^choicevar\\(", xvars)] <- TRUE
+ if (sum(choicevar.indic) > 0){
+ cvarname1.vec <- rep(NA, sum(choicevar.indic))
+ cvarname2.vec <- rep(NA, sum(choicevar.indic))
+ counter <- 0
+ for (i in 1:length(xvars)){
+ if (choicevar.indic[i]){
+ counter <- counter + 1
+ vn2 <- strsplit(xvars[i], ",")
+ vn3 <- strsplit(vn2[[1]], "\\(")
+ vn4 <- strsplit(vn3[[3]], "=")
+ cvarname1 <- vn3[[2]][1]
+ cvarname1 <- strsplit(cvarname1, "\"")[[1]]
+ cvarname1 <- cvarname1[length(cvarname1)]
+ cvarname2 <- vn4[[1]][length(vn4[[1]])]
+ cvarname2 <- strsplit(cvarname2, "\"")[[1]][2]
+ if (! cvarname2 %in% choice.names){
+ stop("choicelevel that was set in choicevar() not consistent with\n observed choice levels in y")
+ }
+ cvarname1.vec[counter] <- cvarname1
+ cvarname2.vec[counter] <- cvarname2
+ xvars[i] <- paste(cvarname1, cvarname2, sep=".")
+ }
+ }
+
+ X.cho <- X[, choicevar.indic]
+ X.cho.mat <- matrix(NA, length(choice.names)*nrow(X),
+ length(unique(cvarname1.vec)))
+ rownames(X.cho.mat) <- rep(rownames(X), rep(length(choice.names), nrow(X)))
+ rownames(X.cho.mat) <- paste(rownames(X.cho.mat), choice.names, sep=".")
+ colnames(X.cho.mat) <- unique(cvarname1.vec)
+ choice.names.n <- rep(choice.names, nrow(X))
+ for (j in 1:length(unique(cvarname1.vec))){
+ for (i in 1:length(cvarname2.vec)){
+ if (colnames(X.cho.mat)[j] == cvarname1.vec[i]){
+ X.cho.mat[choice.names.n==cvarname2.vec[i], j] <-
+ X.cho[,i]
+ }
+ }
+ }
+ }
+
+
+ ## deal with individual specific covariates
+ xvars.ind.mat <- rep(xvars[!choicevar.indic],
+ rep(length(choice.names),
+ sum(!choicevar.indic)))
+ xvars.ind.mat <- paste(xvars.ind.mat, choice.names, sep=".")
+
+ if (sum(!choicevar.indic) > 0){
+ X.ind.mat <- X[,!choicevar.indic] %x% diag(length(choice.names))
+ colnames(X.ind.mat) <- xvars.ind.mat
+ rownames(X.ind.mat) <- rep(rownames(X), rep(length(choice.names), nrow(X)))
+ rownames(X.ind.mat) <- paste(rownames(X.ind.mat), choice.names, sep=".")
+
+ ## delete columns correpsonding to the baseline choice
+ ivarname1 <- strsplit(xvars.ind.mat, "\\.")
+ ivar.keep.indic <- rep(NA, ncol(X.ind.mat))
+ for (i in 1:ncol(X.ind.mat)){
+ ivar.keep.indic[i] <- ivarname1[[i]][length(ivarname1[[i]])] != baseline
+ }
+ X.ind.mat <- X.ind.mat[,ivar.keep.indic]
+ }
+
+ if (sum(choicevar.indic) > 0 & sum(!choicevar.indic) > 0){
+ X <- cbind(X.cho.mat, X.ind.mat)
+ }
+ else if (sum(!choicevar.indic) > 0){
+ X <- X.ind.mat
+ }
+ else if (sum(choicevar.indic) > 0){
+ X <- X.cho.mat
+ }
+ else {
+ stop("X matrix appears to have neither choice-specific nor individual-specific variables.\n")
+ }
+ #Y <- Y.long
+ xvars <- colnames(X)
+ xobs <- rownames(X)
+
+ return(list(Ymat, X, xvars, xobs, number.choices))
+
+}
+
+## dummy function used to handle choice-specific covariates
+"choicevar" <- function(var, varname, choicelevel){
+ junk1 <- varname
+ junk2 <- choicelevel
+ return(var)
+}
+
+
+
+## MNL log-likelihood function (used to get starting values)
+"mnl.loglike" <- function(beta, Y, X){
+ k <- ncol(X)
+ numer <- exp(X %*% beta)
+ numer[Y== -999] <- NA
+ numer.mat <- matrix(numer, nrow(Y), ncol(Y), byrow=TRUE)
+ denom <- apply(numer.mat, 1, sum, na.rm=TRUE)
+ choice.probs <- numer.mat / denom
+ Yna <- Y
+ Yna[Y == -999] <- NA
+ log.like.mat <- log(choice.probs) * Yna
+ log.like <- sum(apply(log.like.mat, 1, sum, na.rm=TRUE))
+ return(log.like)
+}
+
+
+
+
+"MCMCmnl" <-
+ function(formula, baseline=NULL, data = parent.frame(),
+ burnin = 1000, mcmc = 10000, thin=1,
+ mcmc.method = "MH", tune = 1.1, verbose = FALSE, seed = NA,
+ beta.start = NA, b0 = 0, B0 = 0, ...) {
+
+ ## 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 matrix
+ holder <- parse.formula.mnl(formula=formula, baseline=baseline,
+ data=data)
+ Y <- holder[[1]]
+ ## check to make sure baseline category is always available in choiceset
+ if (is.null(baseline)){
+ if (max(Y[,1] == -999) == 1){
+ stop("Baseline choice not available in all choicesets.\n Respecify baseline category and try again.\n")
+ }
+ }
+ else{
+ if (max(Y[,baseline] == -999) == 1){
+ stop("Baseline choice not available in all choicesets.\n Respecify baseline category and try again.\n")
+ }
+ }
+ X <- holder[[2]]
+ xnames <- holder[[3]]
+ xobs <- holder[[4]]
+ number.choices <- holder[[5]]
+ K <- ncol(X) # number of covariates
+
+ ## form the tuning parameter
+ tune <- vector.tune(tune, K)
+ beta.init <- rep(0, K)
+ cat("Calculating MLEs and large sample var-cov matrix.\n")
+ cat("This may take a moment...\n")
+ optim.out <- optim(beta.init, mnl.loglike, method="BFGS",
+ control=list(fnscale=-1),
+ hessian=TRUE, Y=Y, X=X)
+ V <- solve(-1*optim.out$hessian)
+
+ ## starting values and priors
+ if (is.na(beta.start) || is.null(beta.start)){
+ beta.start <- matrix(optim.out$par, K, 1)
+ }
+ else if(is.null(dim(beta.start))) {
+ beta.start <- matrix(beta.start, K, 1)
+ }
+ else if (length(beta.start != K)){
+ stop("beta.start not of appropriate dimension\n")
+ }
+ mvn.prior <- form.mvn.prior(b0, B0, K)
+ b0 <- mvn.prior[[1]]
+ B0 <- mvn.prior[[2]]
+
+ ## define holder for posterior density 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),
+ mcmc=as.integer(mcmc), thin=as.integer(thin),
+ tune=tune, lecuyer=as.integer(lecuyer),
+ seedarray=as.integer(seed.array),
+ lecuyerstream=as.integer(lecuyer.stream),
+ verbose=as.integer(verbose),
+ betastart=beta.start,
+ b0=b0, B0=B0,
+ V=V)
+
+ ## put together matrix and build MCMC object to return
+ output <- form.mcmc.object(posterior, names=xnames,
+ title="MCMCmnl Posterior Density Sample")
+ }
+ else if (mcmc.method=="slice"){
+ ## call C++ code to draw sample
+ auto.Scythe.call(output.object="posterior", cc.fun.name="MCMCmnlslice",
+ sample=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),
+ lecuyerstream=as.integer(lecuyer.stream),
+ verbose=as.integer(verbose), betastart=beta.start,
+ b0=b0, B0=B0, V=V)
+
+ ## put together matrix and build MCMC object to return
+ output <- form.mcmc.object(posterior, names=xnames,
+ title="MCMCmnl Posterior Density Sample")
+
+ }
+
+ return(output)
+
+ }
+
+
+
diff --git a/R/MCMCordfactanal.R b/R/MCMCordfactanal.R
index ea772b5..ed75b0c 100644
--- a/R/MCMCordfactanal.R
+++ b/R/MCMCordfactanal.R
@@ -25,7 +25,7 @@
"MCMCordfactanal" <-
function(x, factors, lambda.constraints=list(),
- data=parent.environment(), burnin = 1000, mcmc = 20000,
+ data=parent.frame(), burnin = 1000, mcmc = 20000,
thin=1, tune=NA, verbose = FALSE, seed = NA,
lambda.start = NA, l0=0, L0=0,
store.lambda=TRUE, store.scores=FALSE,
@@ -87,6 +87,7 @@
mf$... <- NULL
mf$drop.unused.levels <- TRUE
mf[[1]] <- as.name("model.frame")
+ mf$na.action <- 'na.pass'
mf <- eval(mf, sys.frame(sys.parent()))
attributes(mt)$intercept <- 0
Xterm.length <- length(attr(mt, "variables"))
@@ -236,6 +237,8 @@
length(gamma))
}
+
+ accepts <- matrix(0, K, 1)
## Call the C++ code to do the real work
posterior <- .C("ordfactanalpost",
@@ -276,13 +279,20 @@
Lamppreccol = as.integer(ncol(Lambda.prior.prec)),
storelambda = as.integer(store.lambda),
storescores = as.integer(store.scores),
- accepts = as.integer(0),
+ accepts = as.integer(accepts),
+ acceptsrow = as.integer(nrow(accepts)),
+ acceptscol = as.integer(ncol(accepts)),
outswitch = as.integer(case.switch),
PACKAGE="MCMCpack"
)
if(case.switch==1) {
- cat(" overall acceptance rate = ",
- posterior$accepts / ((posterior$burnin+posterior$mcmc)*K), "\n")
+ accepts <- matrix(posterior$accepts, posterior$acceptsrow,
+ posterior$acceptscol, byrow=TRUE)
+ rownames(accepts) <- X.names
+ colnames(accepts) <- ""
+ cat("\n\nAcceptance rates:\n")
+ print(t(accepts) / (posterior$burnin+posterior$mcmc), digits=2,
+ width=6)
}
# put together matrix and build MCMC object to return
diff --git a/R/MCMCtobit.R b/R/MCMCtobit.R
new file mode 100644
index 0000000..21476bb
--- /dev/null
+++ b/R/MCMCtobit.R
@@ -0,0 +1,65 @@
+"MCMCtobit" <-
+function(formula, data=parent.frame(), below = 0.0, above = Inf,
+ burnin = 1000, mcmc = 10000,
+ thin=1, verbose = FALSE, seed = NA, beta.start = NA,
+ b0 = 0, B0 = 0, c0 = 0.001, d0 = 0.001, ...) {
+
+ # checks
+ check.offset(list(...))
+ check.mcmc.parameters(burnin, mcmc, thin)
+ if (!is.numeric(below) | !is.numeric(above)) {
+ cat("Error: Censoring points must be numeric, which includes +-Inf.\n")
+ stop("Please respecify and call ", calling.function(), " again.",
+ call.=FALSE)
+ }
+ if (below >= above) {
+ cat("Error: Truncation points are logically inconsistent.\n")
+ stop("Please respecify and call ", calling.function(), " again.",
+ call.=FALSE)
+ }
+
+ # convert infinite values to finite approximations
+ if(is.infinite(below)) below <- .Machine$double.xmax*-1
+ if(is.infinite(above)) above <- .Machine$double.xmax
+
+ # seeds
+ seeds <- form.seeds(seed)
+ lecuyer <- seeds[[1]]
+ seed.array <- seeds[[2]]
+ lecuyer.stream <- seeds[[3]]
+
+ # 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
+ 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
+ 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),
+ above=as.double(above),
+ 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), betastart=beta.start,
+ b0=b0, B0=B0, c0=as.double(c0), d0=as.double(d0))
+
+ # pull together matrix and build MCMC object to return
+ output <- form.mcmc.object(posterior,
+ names=c(xnames, "sigma2"),
+ title="MCMCtobit Posterior Density Sample")
+ return(output)
+ }
diff --git a/R/automate.R b/R/automate.R
index 7123e65..46949b9 100644
--- a/R/automate.R
+++ b/R/automate.R
@@ -250,7 +250,8 @@
cat("\nThe call to .C is:\n")
draw.sample.call <- parse(text=paste(output.object, " <- ", call))
print(draw.sample.call)
- cat("\nAUTOMATIC TEMPLATE FILE CREATION SUCCEEDED.\n\n")
+ cat("\nAUTOMATIC TEMPLATE FILE CREATION SUCCEEDED.\n")
+ cat("All created templated files placed in ", getwd(), ".\n", sep="")
invokeRestart("abort")
}
diff --git a/R/distn.R b/R/distn.R
index d086212..05ac9db 100644
--- a/R/distn.R
+++ b/R/distn.R
@@ -1,4 +1,11 @@
########## Density Functions and Random Number Generators ##########
+#
+# note: Matthew Fasman has been working on replacing these with
+# functions coded in C++. These have been moved to the tmp
+# directory. We need to get the R seeding issues worked out
+# before these functions can be replaced. In addition, they
+# need to be documented one distribution at a time (see the
+# existing documentation for an example).
##
## Wishart
diff --git a/R/hidden.R b/R/hidden.R
index 9db2bac..6a963e1 100644
--- a/R/hidden.R
+++ b/R/hidden.R
@@ -479,7 +479,7 @@
output <- mcmc(data=holder, start=1,
end=posterior.object$mcmc,
thin=posterior.object$thin)
- varnames(output) <- names
+ varnames(output) <- as.list(names)
attr(output,"title") <- title
return(output)
}
@@ -703,7 +703,12 @@
cat("Error: Vector tuning parameter cannot contain negative values.\n")
stop("Please respecify and call ", calling.function(), " again.",
call.=FALSE)
- }
- return(diag(as.double(mcmc.tune)))
+ }
+ if (length(mcmc.tune)==1){
+ return(matrix(mcmc.tune, 1, 1))
+ }
+ else{
+ return(diag(as.double(mcmc.tune)))
+ }
}
diff --git a/R/zzz.R b/R/zzz.R
index 1c8f07f..417fc49 100644
--- a/R/zzz.R
+++ b/R/zzz.R
@@ -1,6 +1,6 @@
.onAttach <- function(...) {
cat("##\n## Markov Chain Monte Carlo Package (MCMCpack)\n")
- cat("## Copyright (C) 2003, 2004, Andrew D. Martin and Kevin M. Quinn\n")
+ cat("## Copyright (C) 2003-2005, Andrew D. Martin and Kevin M. Quinn\n")
cat("##\n## Support provided by the U.S. National Science Foundation\n")
cat("## (Grants SES-0350646 and SES-0350613)\n##\n")
require(coda, quietly=TRUE)
diff --git a/README b/README
index deb012e..eca6e6c 100644
--- a/README
+++ b/README
@@ -1,8 +1,8 @@
//
-// MCMCpack, Version 0.5-2
+// MCMCpack, Version 0.6-1
//
-Release Date: August 17, 2004
+Release Date: February 18, 2005
This package contains functions for posterior simulation for a number of
statistical models. All simulation is done in compiled C++ written in
@@ -18,7 +18,7 @@ sampling algorithm, and tools for visualization.
//
Andrew D. Martin <admartin at wustl.edu>
-Kevin M. Quinn <kquinn at latte.harvard.edu>
+Kevin M. Quinn <kevin_quinn at harvard.edu>
//
// Compilation
diff --git a/config.status b/config.status
deleted file mode 100755
index 7826a67..0000000
--- a/config.status
+++ /dev/null
@@ -1,712 +0,0 @@
-#! /bin/sh
-# Generated by configure.
-# Run this file to recreate the current configuration.
-# Compiler output produced by configure, useful for debugging
-# configure, is in config.log if it exists.
-
-debug=false
-ac_cs_recheck=false
-ac_cs_silent=false
-SHELL=${CONFIG_SHELL-/bin/sh}
-## --------------------- ##
-## M4sh Initialization. ##
-## --------------------- ##
-
-# Be Bourne compatible
-if test -n "${ZSH_VERSION+set}" && (emulate sh) >/dev/null 2>&1; then
- emulate sh
- NULLCMD=:
- # Zsh 3.x and 4.x performs word splitting on ${1+"$@"}, which
- # is contrary to our usage. Disable this feature.
- alias -g '${1+"$@"}'='"$@"'
-elif test -n "${BASH_VERSION+set}" && (set -o posix) >/dev/null 2>&1; then
- set -o posix
-fi
-DUALCASE=1; export DUALCASE # for MKS sh
-
-# Support unset when possible.
-if ( (MAIL=60; unset MAIL) || exit) >/dev/null 2>&1; then
- as_unset=unset
-else
- as_unset=false
-fi
-
-
-# Work around bugs in pre-3.0 UWIN ksh.
-$as_unset ENV MAIL MAILPATH
-PS1='$ '
-PS2='> '
-PS4='+ '
-
-# NLS nuisances.
-for as_var in \
- LANG LANGUAGE LC_ADDRESS LC_ALL LC_COLLATE LC_CTYPE LC_IDENTIFICATION \
- LC_MEASUREMENT LC_MESSAGES LC_MONETARY LC_NAME LC_NUMERIC LC_PAPER \
- LC_TELEPHONE LC_TIME
-do
- if (set +x; test -z "`(eval $as_var=C; export $as_var) 2>&1`"); then
- eval $as_var=C; export $as_var
- else
- $as_unset $as_var
- fi
-done
-
-# Required to use basename.
-if expr a : '\(a\)' >/dev/null 2>&1; then
- as_expr=expr
-else
- as_expr=false
-fi
-
-if (basename /) >/dev/null 2>&1 && test "X`basename / 2>&1`" = "X/"; then
- as_basename=basename
-else
- as_basename=false
-fi
-
-
-# Name of the executable.
-as_me=`$as_basename "$0" ||
-$as_expr X/"$0" : '.*/\([^/][^/]*\)/*$' \| \
- X"$0" : 'X\(//\)$' \| \
- X"$0" : 'X\(/\)$' \| \
- . : '\(.\)' 2>/dev/null ||
-echo X/"$0" |
- sed '/^.*\/\([^/][^/]*\)\/*$/{ s//\1/; q; }
- /^X\/\(\/\/\)$/{ s//\1/; q; }
- /^X\/\(\/\).*/{ s//\1/; q; }
- s/.*/./; q'`
-
-
-# PATH needs CR, and LINENO needs CR and PATH.
-# Avoid depending upon Character Ranges.
-as_cr_letters='abcdefghijklmnopqrstuvwxyz'
-as_cr_LETTERS='ABCDEFGHIJKLMNOPQRSTUVWXYZ'
-as_cr_Letters=$as_cr_letters$as_cr_LETTERS
-as_cr_digits='0123456789'
-as_cr_alnum=$as_cr_Letters$as_cr_digits
-
-# The user is always right.
-if test "${PATH_SEPARATOR+set}" != set; then
- echo "#! /bin/sh" >conf$$.sh
- echo "exit 0" >>conf$$.sh
- chmod +x conf$$.sh
- if (PATH="/nonexistent;."; conf$$.sh) >/dev/null 2>&1; then
- PATH_SEPARATOR=';'
- else
- PATH_SEPARATOR=:
- fi
- rm -f conf$$.sh
-fi
-
-
- as_lineno_1=$LINENO
- as_lineno_2=$LINENO
- as_lineno_3=`(expr $as_lineno_1 + 1) 2>/dev/null`
- test "x$as_lineno_1" != "x$as_lineno_2" &&
- test "x$as_lineno_3" = "x$as_lineno_2" || {
- # Find who we are. Look in the path if we contain no path at all
- # relative or not.
- case $0 in
- *[\\/]* ) as_myself=$0 ;;
- *) as_save_IFS=$IFS; IFS=$PATH_SEPARATOR
-for as_dir in $PATH
-do
- IFS=$as_save_IFS
- test -z "$as_dir" && as_dir=.
- test -r "$as_dir/$0" && as_myself=$as_dir/$0 && break
-done
-
- ;;
- esac
- # We did not find ourselves, most probably we were run as `sh COMMAND'
- # in which case we are not to be found in the path.
- if test "x$as_myself" = x; then
- as_myself=$0
- fi
- if test ! -f "$as_myself"; then
- { { echo "$as_me:$LINENO: error: cannot find myself; rerun with an absolute path" >&5
-echo "$as_me: error: cannot find myself; rerun with an absolute path" >&2;}
- { (exit 1); exit 1; }; }
- fi
- case $CONFIG_SHELL in
- '')
- as_save_IFS=$IFS; IFS=$PATH_SEPARATOR
-for as_dir in /bin$PATH_SEPARATOR/usr/bin$PATH_SEPARATOR$PATH
-do
- IFS=$as_save_IFS
- test -z "$as_dir" && as_dir=.
- for as_base in sh bash ksh sh5; do
- case $as_dir in
- /*)
- if ("$as_dir/$as_base" -c '
- as_lineno_1=$LINENO
- as_lineno_2=$LINENO
- as_lineno_3=`(expr $as_lineno_1 + 1) 2>/dev/null`
- test "x$as_lineno_1" != "x$as_lineno_2" &&
- test "x$as_lineno_3" = "x$as_lineno_2" ') 2>/dev/null; then
- $as_unset BASH_ENV || test "${BASH_ENV+set}" != set || { BASH_ENV=; export BASH_ENV; }
- $as_unset ENV || test "${ENV+set}" != set || { ENV=; export ENV; }
- CONFIG_SHELL=$as_dir/$as_base
- export CONFIG_SHELL
- exec "$CONFIG_SHELL" "$0" ${1+"$@"}
- fi;;
- esac
- done
-done
-;;
- esac
-
- # Create $as_me.lineno as a copy of $as_myself, but with $LINENO
- # uniformly replaced by the line number. The first 'sed' inserts a
- # line-number line before each line; the second 'sed' does the real
- # work. The second script uses 'N' to pair each line-number line
- # with the numbered line, and appends trailing '-' during
- # substitution so that $LINENO is not a special case at line end.
- # (Raja R Harinath suggested sed '=', and Paul Eggert wrote the
- # second 'sed' script. Blame Lee E. McMahon for sed's syntax. :-)
- sed '=' <$as_myself |
- sed '
- N
- s,$,-,
- : loop
- s,^\(['$as_cr_digits']*\)\(.*\)[$]LINENO\([^'$as_cr_alnum'_]\),\1\2\1\3,
- t loop
- s,-$,,
- s,^['$as_cr_digits']*\n,,
- ' >$as_me.lineno &&
- chmod +x $as_me.lineno ||
- { { echo "$as_me:$LINENO: error: cannot create $as_me.lineno; rerun with a POSIX shell" >&5
-echo "$as_me: error: cannot create $as_me.lineno; rerun with a POSIX shell" >&2;}
- { (exit 1); exit 1; }; }
-
- # Don't try to exec as it changes $[0], causing all sort of problems
- # (the dirname of $[0] is not the place where we might find the
- # original and so on. Autoconf is especially sensible to this).
- . ./$as_me.lineno
- # Exit status is that of the last command.
- exit
-}
-
-
-case `echo "testing\c"; echo 1,2,3`,`echo -n testing; echo 1,2,3` in
- *c*,-n*) ECHO_N= ECHO_C='
-' ECHO_T=' ' ;;
- *c*,* ) ECHO_N=-n ECHO_C= ECHO_T= ;;
- *) ECHO_N= ECHO_C='\c' ECHO_T= ;;
-esac
-
-if expr a : '\(a\)' >/dev/null 2>&1; then
- as_expr=expr
-else
- as_expr=false
-fi
-
-rm -f conf$$ conf$$.exe conf$$.file
-echo >conf$$.file
-if ln -s conf$$.file conf$$ 2>/dev/null; then
- # We could just check for DJGPP; but this test a) works b) is more generic
- # and c) will remain valid once DJGPP supports symlinks (DJGPP 2.04).
- if test -f conf$$.exe; then
- # Don't use ln at all; we don't have any links
- as_ln_s='cp -p'
- else
- as_ln_s='ln -s'
- fi
-elif ln conf$$.file conf$$ 2>/dev/null; then
- as_ln_s=ln
-else
- as_ln_s='cp -p'
-fi
-rm -f conf$$ conf$$.exe conf$$.file
-
-if mkdir -p . 2>/dev/null; then
- as_mkdir_p=:
-else
- test -d ./-p && rmdir ./-p
- as_mkdir_p=false
-fi
-
-as_executable_p="test -f"
-
-# Sed expression to map a string onto a valid CPP name.
-as_tr_cpp="eval sed 'y%*$as_cr_letters%P$as_cr_LETTERS%;s%[^_$as_cr_alnum]%_%g'"
-
-# Sed expression to map a string onto a valid variable name.
-as_tr_sh="eval sed 'y%*+%pp%;s%[^_$as_cr_alnum]%_%g'"
-
-
-# IFS
-# We need space, tab and new line, in precisely that order.
-as_nl='
-'
-IFS=" $as_nl"
-
-# CDPATH.
-$as_unset CDPATH
-
-exec 6>&1
-
-# Open the log real soon, to keep \$[0] and so on meaningful, and to
-# report actual input values of CONFIG_FILES etc. instead of their
-# values after options handling. Logging --version etc. is OK.
-exec 5>>config.log
-{
- echo
- sed 'h;s/./-/g;s/^.../## /;s/...$/ ##/;p;x;p;x' <<_ASBOX
-## Running $as_me. ##
-_ASBOX
-} >&5
-cat >&5 <<_CSEOF
-
-This file was extended by $as_me, which was
-generated by GNU Autoconf 2.59. Invocation command line was
-
- CONFIG_FILES = $CONFIG_FILES
- CONFIG_HEADERS = $CONFIG_HEADERS
- CONFIG_LINKS = $CONFIG_LINKS
- CONFIG_COMMANDS = $CONFIG_COMMANDS
- $ $0 $@
-
-_CSEOF
-echo "on `(hostname || uname -n) 2>/dev/null | sed 1q`" >&5
-echo >&5
-config_files=" src/Makevars"
-
-ac_cs_usage="\
-\`$as_me' instantiates files from templates according to the
-current configuration.
-
-Usage: $0 [OPTIONS] [FILE]...
-
- -h, --help print this help, then exit
- -V, --version print version number, then exit
- -q, --quiet do not print progress messages
- -d, --debug don't remove temporary files
- --recheck update $as_me by reconfiguring in the same conditions
- --file=FILE[:TEMPLATE]
- instantiate the configuration file FILE
-
-Configuration files:
-$config_files
-
-Report bugs to <bug-autoconf at gnu.org>."
-ac_cs_version="\
-config.status
-configured by ./configure, generated by GNU Autoconf 2.59,
- with options \"\"
-
-Copyright (C) 2003 Free Software Foundation, Inc.
-This config.status script is free software; the Free Software Foundation
-gives unlimited permission to copy, distribute and modify it."
-srcdir=.
-# If no file are specified by the user, then we need to provide default
-# value. By we need to know if files were specified by the user.
-ac_need_defaults=:
-while test $# != 0
-do
- case $1 in
- --*=*)
- ac_option=`expr "x$1" : 'x\([^=]*\)='`
- ac_optarg=`expr "x$1" : 'x[^=]*=\(.*\)'`
- ac_shift=:
- ;;
- -*)
- ac_option=$1
- ac_optarg=$2
- ac_shift=shift
- ;;
- *) # This is not an option, so the user has probably given explicit
- # arguments.
- ac_option=$1
- ac_need_defaults=false;;
- esac
-
- case $ac_option in
- # Handling of the options.
- -recheck | --recheck | --rechec | --reche | --rech | --rec | --re | --r)
- ac_cs_recheck=: ;;
- --version | --vers* | -V )
- echo "$ac_cs_version"; exit 0 ;;
- --he | --h)
- # Conflict between --help and --header
- { { echo "$as_me:$LINENO: error: ambiguous option: $1
-Try \`$0 --help' for more information." >&5
-echo "$as_me: error: ambiguous option: $1
-Try \`$0 --help' for more information." >&2;}
- { (exit 1); exit 1; }; };;
- --help | --hel | -h )
- echo "$ac_cs_usage"; exit 0 ;;
- --debug | --d* | -d )
- debug=: ;;
- --file | --fil | --fi | --f )
- $ac_shift
- CONFIG_FILES="$CONFIG_FILES $ac_optarg"
- ac_need_defaults=false;;
- --header | --heade | --head | --hea )
- $ac_shift
- CONFIG_HEADERS="$CONFIG_HEADERS $ac_optarg"
- ac_need_defaults=false;;
- -q | -quiet | --quiet | --quie | --qui | --qu | --q \
- | -silent | --silent | --silen | --sile | --sil | --si | --s)
- ac_cs_silent=: ;;
-
- # This is an error.
- -*) { { echo "$as_me:$LINENO: error: unrecognized option: $1
-Try \`$0 --help' for more information." >&5
-echo "$as_me: error: unrecognized option: $1
-Try \`$0 --help' for more information." >&2;}
- { (exit 1); exit 1; }; } ;;
-
- *) ac_config_targets="$ac_config_targets $1" ;;
-
- esac
- shift
-done
-
-ac_configure_extra_args=
-
-if $ac_cs_silent; then
- exec 6>/dev/null
- ac_configure_extra_args="$ac_configure_extra_args --silent"
-fi
-
-if $ac_cs_recheck; then
- echo "running /bin/sh ./configure " $ac_configure_extra_args " --no-create --no-recursion" >&6
- exec /bin/sh ./configure $ac_configure_extra_args --no-create --no-recursion
-fi
-
-for ac_config_target in $ac_config_targets
-do
- case "$ac_config_target" in
- # Handling of arguments.
- "src/Makevars" ) CONFIG_FILES="$CONFIG_FILES src/Makevars" ;;
- *) { { echo "$as_me:$LINENO: error: invalid argument: $ac_config_target" >&5
-echo "$as_me: error: invalid argument: $ac_config_target" >&2;}
- { (exit 1); exit 1; }; };;
- esac
-done
-
-# If the user did not use the arguments to specify the items to instantiate,
-# then the envvar interface is used. Set only those that are not.
-# We use the long form for the default assignment because of an extremely
-# bizarre bug on SunOS 4.1.3.
-if $ac_need_defaults; then
- test "${CONFIG_FILES+set}" = set || CONFIG_FILES=$config_files
-fi
-
-# Have a temporary directory for convenience. Make it in the build tree
-# simply because there is no reason to put it here, and in addition,
-# creating and moving files from /tmp can sometimes cause problems.
-# Create a temporary directory, and hook for its removal unless debugging.
-$debug ||
-{
- trap 'exit_status=$?; rm -rf $tmp && exit $exit_status' 0
- trap '{ (exit 1); exit 1; }' 1 2 13 15
-}
-
-# Create a (secure) tmp directory for tmp files.
-
-{
- tmp=`(umask 077 && mktemp -d -q "./confstatXXXXXX") 2>/dev/null` &&
- test -n "$tmp" && test -d "$tmp"
-} ||
-{
- tmp=./confstat$$-$RANDOM
- (umask 077 && mkdir $tmp)
-} ||
-{
- echo "$me: cannot create a temporary directory in ." >&2
- { (exit 1); exit 1; }
-}
-
-
-#
-# CONFIG_FILES section.
-#
-
-# No need to generate the scripts if there are no CONFIG_FILES.
-# This happens for instance when ./config.status config.h
-if test -n "$CONFIG_FILES"; then
- # Protect against being on the right side of a sed subst in config.status.
- sed 's/,@/@@/; s/@,/@@/; s/,;t t$/@;t t/; /@;t t$/s/[\\&,]/\\&/g;
- s/@@/,@/; s/@@/@,/; s/@;t t$/,;t t/' >$tmp/subs.sed <<\CEOF
-s, at SHELL@,/bin/sh,;t t
-s, at PATH_SEPARATOR@,:,;t t
-s, at PACKAGE_NAME@,,;t t
-s, at PACKAGE_TARNAME@,,;t t
-s, at PACKAGE_VERSION@,,;t t
-s, at PACKAGE_STRING@,,;t t
-s, at PACKAGE_BUGREPORT@,,;t t
-s, at exec_prefix@,${prefix},;t t
-s, at prefix@,/usr/local,;t t
-s, at program_transform_name@,s,x,x,,;t t
-s, at bindir@,${exec_prefix}/bin,;t t
-s, at sbindir@,${exec_prefix}/sbin,;t t
-s, at libexecdir@,${exec_prefix}/libexec,;t t
-s, at datadir@,${prefix}/share,;t t
-s, at sysconfdir@,${prefix}/etc,;t t
-s, at sharedstatedir@,${prefix}/com,;t t
-s, at localstatedir@,${prefix}/var,;t t
-s, at libdir@,${exec_prefix}/lib,;t t
-s, at includedir@,${prefix}/include,;t t
-s, at oldincludedir@,/usr/include,;t t
-s, at infodir@,${prefix}/info,;t t
-s, at mandir@,${prefix}/man,;t t
-s, at build_alias@,,;t t
-s, at host_alias@,,;t t
-s, at target_alias@,,;t t
-s, at DEFS@,-DPACKAGE_NAME=\"\" -DPACKAGE_TARNAME=\"\" -DPACKAGE_VERSION=\"\" -DPACKAGE_STRING=\"\" -DPACKAGE_BUGREPORT=\"\" -DSTDC_HEADERS=1 -DHAVE_SYS_TYPES_H=1 -DHAVE_SYS_STAT_H=1 -DHAVE_STDLIB_H=1 -DHAVE_STRING_H=1 -DHAVE_MEMORY_H=1 -DHAVE_STRINGS_H=1 -DHAVE_INTTYPES_H=1 -DHAVE_STDINT_H=1 -DHAVE_UNISTD_H=1 -DHAVE_TRUNC=1 ,;t t
-s, at ECHO_C@,,;t t
-s, at ECHO_N@,-n,;t t
-s, at ECHO_T@,,;t t
-s, at LIBS@,,;t t
-s, at CXX@,g++,;t t
-s, at CXXFLAGS@,-g -O2,;t t
-s, at LDFLAGS@,,;t t
-s, at CPPFLAGS@,,;t t
-s, at ac_ct_CXX@,,;t t
-s, at EXEEXT@,,;t t
-s, at OBJEXT@,o,;t t
-s, at CC@,gcc,;t t
-s, at CFLAGS@,-g -O2,;t t
-s, at ac_ct_CC@,gcc,;t t
-s, at CPP@,gcc -E,;t t
-s, at EGREP@,grep -E,;t t
-s, at MV_HAVE_IEEEFP_H@,,;t t
-s, at MV_HAVE_TRUNC@,-DHAVE_TRUNC,;t t
-s, at LIBOBJS@,,;t t
-s, at LTLIBOBJS@,,;t t
-CEOF
-
- # Split the substitutions into bite-sized pieces for seds with
- # small command number limits, like on Digital OSF/1 and HP-UX.
- ac_max_sed_lines=48
- ac_sed_frag=1 # Number of current file.
- ac_beg=1 # First line for current file.
- ac_end=$ac_max_sed_lines # Line after last line for current file.
- ac_more_lines=:
- ac_sed_cmds=
- while $ac_more_lines; do
- if test $ac_beg -gt 1; then
- sed "1,${ac_beg}d; ${ac_end}q" $tmp/subs.sed >$tmp/subs.frag
- else
- sed "${ac_end}q" $tmp/subs.sed >$tmp/subs.frag
- fi
- if test ! -s $tmp/subs.frag; then
- ac_more_lines=false
- else
- # The purpose of the label and of the branching condition is to
- # speed up the sed processing (if there are no `@' at all, there
- # is no need to browse any of the substitutions).
- # These are the two extra sed commands mentioned above.
- (echo ':t
- /@[a-zA-Z_][a-zA-Z_0-9]*@/!b' && cat $tmp/subs.frag) >$tmp/subs-$ac_sed_frag.sed
- if test -z "$ac_sed_cmds"; then
- ac_sed_cmds="sed -f $tmp/subs-$ac_sed_frag.sed"
- else
- ac_sed_cmds="$ac_sed_cmds | sed -f $tmp/subs-$ac_sed_frag.sed"
- fi
- ac_sed_frag=`expr $ac_sed_frag + 1`
- ac_beg=$ac_end
- ac_end=`expr $ac_end + $ac_max_sed_lines`
- fi
- done
- if test -z "$ac_sed_cmds"; then
- ac_sed_cmds=cat
- fi
-fi # test -n "$CONFIG_FILES"
-
-for ac_file in : $CONFIG_FILES; do test "x$ac_file" = x: && continue
- # Support "outfile[:infile[:infile...]]", defaulting infile="outfile.in".
- case $ac_file in
- - | *:- | *:-:* ) # input from stdin
- cat >$tmp/stdin
- ac_file_in=`echo "$ac_file" | sed 's,[^:]*:,,'`
- ac_file=`echo "$ac_file" | sed 's,:.*,,'` ;;
- *:* ) ac_file_in=`echo "$ac_file" | sed 's,[^:]*:,,'`
- ac_file=`echo "$ac_file" | sed 's,:.*,,'` ;;
- * ) ac_file_in=$ac_file.in ;;
- esac
-
- # Compute @srcdir@, @top_srcdir@, and @INSTALL@ for subdirectories.
- ac_dir=`(dirname "$ac_file") 2>/dev/null ||
-$as_expr X"$ac_file" : 'X\(.*[^/]\)//*[^/][^/]*/*$' \| \
- X"$ac_file" : 'X\(//\)[^/]' \| \
- X"$ac_file" : 'X\(//\)$' \| \
- X"$ac_file" : 'X\(/\)' \| \
- . : '\(.\)' 2>/dev/null ||
-echo X"$ac_file" |
- sed '/^X\(.*[^/]\)\/\/*[^/][^/]*\/*$/{ s//\1/; q; }
- /^X\(\/\/\)[^/].*/{ s//\1/; q; }
- /^X\(\/\/\)$/{ s//\1/; q; }
- /^X\(\/\).*/{ s//\1/; q; }
- s/.*/./; q'`
- { if $as_mkdir_p; then
- mkdir -p "$ac_dir"
- else
- as_dir="$ac_dir"
- as_dirs=
- while test ! -d "$as_dir"; do
- as_dirs="$as_dir $as_dirs"
- as_dir=`(dirname "$as_dir") 2>/dev/null ||
-$as_expr X"$as_dir" : 'X\(.*[^/]\)//*[^/][^/]*/*$' \| \
- X"$as_dir" : 'X\(//\)[^/]' \| \
- X"$as_dir" : 'X\(//\)$' \| \
- X"$as_dir" : 'X\(/\)' \| \
- . : '\(.\)' 2>/dev/null ||
-echo X"$as_dir" |
- sed '/^X\(.*[^/]\)\/\/*[^/][^/]*\/*$/{ s//\1/; q; }
- /^X\(\/\/\)[^/].*/{ s//\1/; q; }
- /^X\(\/\/\)$/{ s//\1/; q; }
- /^X\(\/\).*/{ s//\1/; q; }
- s/.*/./; q'`
- done
- test ! -n "$as_dirs" || mkdir $as_dirs
- fi || { { echo "$as_me:$LINENO: error: cannot create directory \"$ac_dir\"" >&5
-echo "$as_me: error: cannot create directory \"$ac_dir\"" >&2;}
- { (exit 1); exit 1; }; }; }
-
- ac_builddir=.
-
-if test "$ac_dir" != .; then
- ac_dir_suffix=/`echo "$ac_dir" | sed 's,^\.[\\/],,'`
- # A "../" for each directory in $ac_dir_suffix.
- ac_top_builddir=`echo "$ac_dir_suffix" | sed 's,/[^\\/]*,../,g'`
-else
- ac_dir_suffix= ac_top_builddir=
-fi
-
-case $srcdir in
- .) # No --srcdir option. We are building in place.
- ac_srcdir=.
- if test -z "$ac_top_builddir"; then
- ac_top_srcdir=.
- else
- ac_top_srcdir=`echo $ac_top_builddir | sed 's,/$,,'`
- fi ;;
- [\\/]* | ?:[\\/]* ) # Absolute path.
- ac_srcdir=$srcdir$ac_dir_suffix;
- ac_top_srcdir=$srcdir ;;
- *) # Relative path.
- ac_srcdir=$ac_top_builddir$srcdir$ac_dir_suffix
- ac_top_srcdir=$ac_top_builddir$srcdir ;;
-esac
-
-# Do not use `cd foo && pwd` to compute absolute paths, because
-# the directories may not exist.
-case `pwd` in
-.) ac_abs_builddir="$ac_dir";;
-*)
- case "$ac_dir" in
- .) ac_abs_builddir=`pwd`;;
- [\\/]* | ?:[\\/]* ) ac_abs_builddir="$ac_dir";;
- *) ac_abs_builddir=`pwd`/"$ac_dir";;
- esac;;
-esac
-case $ac_abs_builddir in
-.) ac_abs_top_builddir=${ac_top_builddir}.;;
-*)
- case ${ac_top_builddir}. in
- .) ac_abs_top_builddir=$ac_abs_builddir;;
- [\\/]* | ?:[\\/]* ) ac_abs_top_builddir=${ac_top_builddir}.;;
- *) ac_abs_top_builddir=$ac_abs_builddir/${ac_top_builddir}.;;
- esac;;
-esac
-case $ac_abs_builddir in
-.) ac_abs_srcdir=$ac_srcdir;;
-*)
- case $ac_srcdir in
- .) ac_abs_srcdir=$ac_abs_builddir;;
- [\\/]* | ?:[\\/]* ) ac_abs_srcdir=$ac_srcdir;;
- *) ac_abs_srcdir=$ac_abs_builddir/$ac_srcdir;;
- esac;;
-esac
-case $ac_abs_builddir in
-.) ac_abs_top_srcdir=$ac_top_srcdir;;
-*)
- case $ac_top_srcdir in
- .) ac_abs_top_srcdir=$ac_abs_builddir;;
- [\\/]* | ?:[\\/]* ) ac_abs_top_srcdir=$ac_top_srcdir;;
- *) ac_abs_top_srcdir=$ac_abs_builddir/$ac_top_srcdir;;
- esac;;
-esac
-
-
-
- # Let's still pretend it is `configure' which instantiates (i.e., don't
- # use $as_me), people would be surprised to read:
- # /* config.h. Generated by config.status. */
- if test x"$ac_file" = x-; then
- configure_input=
- else
- configure_input="$ac_file. "
- fi
- configure_input=$configure_input"Generated from `echo $ac_file_in |
- sed 's,.*/,,'` by configure."
-
- # First look for the input files in the build tree, otherwise in the
- # src tree.
- ac_file_inputs=`IFS=:
- for f in $ac_file_in; do
- case $f in
- -) echo $tmp/stdin ;;
- [\\/$]*)
- # Absolute (can't be DOS-style, as IFS=:)
- test -f "$f" || { { echo "$as_me:$LINENO: error: cannot find input file: $f" >&5
-echo "$as_me: error: cannot find input file: $f" >&2;}
- { (exit 1); exit 1; }; }
- echo "$f";;
- *) # Relative
- if test -f "$f"; then
- # Build tree
- echo "$f"
- elif test -f "$srcdir/$f"; then
- # Source tree
- echo "$srcdir/$f"
- else
- # /dev/null tree
- { { echo "$as_me:$LINENO: error: cannot find input file: $f" >&5
-echo "$as_me: error: cannot find input file: $f" >&2;}
- { (exit 1); exit 1; }; }
- fi;;
- esac
- done` || { (exit 1); exit 1; }
-
- if test x"$ac_file" != x-; then
- { echo "$as_me:$LINENO: creating $ac_file" >&5
-echo "$as_me: creating $ac_file" >&6;}
- rm -f "$ac_file"
- fi
- sed "/^[ ]*VPATH[ ]*=/{
-s/:*\$(srcdir):*/:/;
-s/:*\${srcdir}:*/:/;
-s/:*@srcdir@:*/:/;
-s/^\([^=]*=[ ]*\):*/\1/;
-s/:*$//;
-s/^[^=]*=[ ]*$//;
-}
-
-:t
-/@[a-zA-Z_][a-zA-Z_0-9]*@/!b
-s, at configure_input@,$configure_input,;t t
-s, at srcdir@,$ac_srcdir,;t t
-s, at abs_srcdir@,$ac_abs_srcdir,;t t
-s, at top_srcdir@,$ac_top_srcdir,;t t
-s, at abs_top_srcdir@,$ac_abs_top_srcdir,;t t
-s, at builddir@,$ac_builddir,;t t
-s, at abs_builddir@,$ac_abs_builddir,;t t
-s, at top_builddir@,$ac_top_builddir,;t t
-s, at abs_top_builddir@,$ac_abs_top_builddir,;t t
-" $ac_file_inputs | (eval "$ac_sed_cmds") >$tmp/out
- rm -f $tmp/stdin
- if test x"$ac_file" != x-; then
- mv $tmp/out $ac_file
- else
- cat $tmp/out
- rm -f $tmp/out
- fi
-
-done
-
-{ (exit 0); exit 0; }
diff --git a/data/Nethvote.rda b/data/Nethvote.rda
new file mode 100644
index 0000000..ea8570d
Binary files /dev/null and b/data/Nethvote.rda differ
diff --git a/man/MCMCdynamicEI.Rd b/man/MCMCdynamicEI.Rd
index b93496e..1f000d6 100644
--- a/man/MCMCdynamicEI.Rd
+++ b/man/MCMCdynamicEI.Rd
@@ -1,6 +1,6 @@
\name{MCMCdynamicEI}
\alias{MCMCdynamicEI}
-\title{Markov chain Monte Carlo for Quinn's Dynamic Ecological
+\title{Markov Chain Monte Carlo for Quinn's Dynamic Ecological
Inference Model}
\description{
MCMCdynamicEI is used to fit Quinn's dynamic ecological inference
diff --git a/man/MCMCfactanal.Rd b/man/MCMCfactanal.Rd
index 29457c5..158972c 100644
--- a/man/MCMCfactanal.Rd
+++ b/man/MCMCfactanal.Rd
@@ -1,6 +1,6 @@
\name{MCMCfactanal}
\alias{MCMCfactanal}
-\title{Markov chain Monte Carlo for Normal Theory Factor Analysis Model}
+\title{Markov Chain Monte Carlo for Normal Theory Factor Analysis Model}
\description{
This function generates a posterior density sample from Normal theory
factor analysis model. Normal priors are assumed on the factor
@@ -13,7 +13,7 @@
\usage{
MCMCfactanal(x, factors, lambda.constraints=list(),
- data=parent.environment(), burnin = 1000, mcmc = 20000,
+ data=parent.frame(), burnin = 1000, mcmc = 20000,
thin=1, verbose = FALSE, seed = NA,
lambda.start = NA, psi.start = NA,
l0=0, L0=0, a0=0.001, b0=0.001,
@@ -148,7 +148,9 @@ MCMCfactanal(x, factors, lambda.constraints=list(),
compiled C++ code to maximize efficiency. Please consult the
coda documentation for a comprehensive list of functions that
can be used to analyze the posterior density sample.
- }
+
+ As is the case with all measurement models, make sure that you have plenty
+ of free memory, especially when storing the scores.
}
\references{
diff --git a/man/MCMChierEI.Rd b/man/MCMChierEI.Rd
index 081821d..d1e010f 100644
--- a/man/MCMChierEI.Rd
+++ b/man/MCMChierEI.Rd
@@ -1,6 +1,6 @@
\name{MCMChierEI}
\alias{MCMChierEI}
-\title{Markov chain Monte Carlo for Wakefield's Hierarchial Ecological
+\title{Markov Chain Monte Carlo for Wakefield's Hierarchial Ecological
Inference Model}
\description{
`MCMChierEI' is used to fit Wakefield's hierarchical ecological inference
diff --git a/man/MCMCirt1d.Rd b/man/MCMCirt1d.Rd
index c6866c9..74ceede 100644
--- a/man/MCMCirt1d.Rd
+++ b/man/MCMCirt1d.Rd
@@ -1,6 +1,6 @@
\name{MCMCirt1d}
\alias{MCMCirt1d}
-\title{Markov chain Monte Carlo for One Dimensional Item Response Theory
+\title{Markov Chain Monte Carlo for One Dimensional Item Response Theory
Model}
\description{
This function generates a posterior density sample from a one
@@ -151,6 +151,9 @@ MCMCirt1d(datamatrix, theta.constraints=list(), burnin = 1000,
The model is identified by the proper priors on the item parameters
and constraints placed on the ability parameters.
+
+ As is the case with all measurement models, make sure that you have plenty
+ of free memory, especially when storing the item parameters.
}
\references{
diff --git a/man/MCMCirtKd.Rd b/man/MCMCirtKd.Rd
index d52b888..bc2296c 100644
--- a/man/MCMCirtKd.Rd
+++ b/man/MCMCirtKd.Rd
@@ -1,6 +1,6 @@
\name{MCMCirtKd}
\alias{MCMCirtKd}
-\title{Markov chain Monte Carlo for K-Dimensional Item Response Theory
+\title{Markov Chain Monte Carlo for K-Dimensional Item Response Theory
Model}
\description{
This function generates a posterior density sample from a
@@ -21,7 +21,7 @@ MCMCirtKd(datamatrix, dimensions, item.constraints=list(),
\arguments{
\item{datamatrix}{The matrix of data. Must be 0, 1, or missing values.
- It is of dimensionality items by subjects.}
+ It is of dimensionality subjects by items.}
\item{dimensions}{The number of dimensions in the latent space.}
@@ -165,10 +165,14 @@ MCMCirtKd(datamatrix, dimensions, item.constraints=list(),
The model is identified by the constraints on the item parameters
(see Jackman 2001). The user cannot place constraints on the subect
abilities. This identification scheme differs from that in
- \code{MCMCirt1d}, which uses a single directional constraint on
- one subject ability. However, in our experience, using subject
+ \code{MCMCirt1d}, which uses constraints on
+ the subject abilities to identify the model.
+ In our experience, using subject
ability constraints for models in greater than one dimension does not work
particularly well.
+
+ As is the case with all measurement models, make sure that you have plenty
+ of free memory, especially when storing the item parameters.
}
\references{
@@ -192,6 +196,9 @@ MCMCirtKd(datamatrix, dimensions, item.constraints=list(),
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.
}
@@ -199,7 +206,7 @@ MCMCirtKd(datamatrix, dimensions, item.constraints=list(),
\dontrun{
data(SupremeCourt)
# note that the rownames (the item names) are "1", "2", etc
- posterior1 <- MCMCirtKd(SupremeCourt, dimensions=1,
+ posterior1 <- MCMCirtKd(t(SupremeCourt), dimensions=1,
burnin=5000, mcmc=50000, thin=10,
B0=.25, store.item=TRUE,
item.constraints=list("1"=list(2,"-")))
@@ -211,7 +218,7 @@ MCMCirtKd(datamatrix, dimensions, item.constraints=list(),
rownames(Senate) <- Senate$member
# note that we need to transpose the data to get
# the bills on the rows
- posterior2 <- MCMCirtKd(t(Senate[,6:677]), dimensions=2,
+ posterior2 <- MCMCirtKd(Senate[,6:677], dimensions=2,
burnin=5000, mcmc=50000, thin=10,
item.constraints=list(rc2=list(2,"-"), rc2=c(3,0),
rc3=list(3,"-")),
diff --git a/man/MCMClogit.Rd b/man/MCMClogit.Rd
index dd0070e..d1ae4bd 100644
--- a/man/MCMClogit.Rd
+++ b/man/MCMClogit.Rd
@@ -1,6 +1,6 @@
\name{MCMClogit}
\alias{MCMClogit}
-\title{Markov chain Monte Carlo for Logistic Regression}
+\title{Markov Chain Monte Carlo for Logistic Regression}
\description{
This function generates a posterior density sample
from a logistic regression model using a random walk Metropolis
diff --git a/man/MCMCmixfactanal.Rd b/man/MCMCmixfactanal.Rd
index cfff474..ae7b3eb 100644
--- a/man/MCMCmixfactanal.Rd
+++ b/man/MCMCmixfactanal.Rd
@@ -1,6 +1,6 @@
\name{MCMCmixfactanal}
\alias{MCMCmixfactanal}
-\title{Markov chain Monte Carlo for Mixed Data Factor Analysis Model}
+\title{Markov Chain Monte Carlo for Mixed Data Factor Analysis Model}
\description{
This function generates a posterior density sample from a mixed data
(both continuous and ordinal) factor analysis model. Normal priors are
@@ -14,7 +14,7 @@
\usage{
MCMCmixfactanal(x, factors, lambda.constraints=list(),
- data=parent.environment(), burnin = 1000, mcmc = 20000,
+ data=parent.frame(), burnin = 1000, mcmc = 20000,
thin=1, tune=NA, verbose = FALSE, seed = NA,
lambda.start = NA, psi.start=NA,
l0=0, L0=0, a0=0.001, b0=0.001,
@@ -25,7 +25,8 @@ MCMCmixfactanal(x, factors, lambda.constraints=list(),
\arguments{
\item{x}{A one-sided formula containing the
manifest variables. Ordinal (including dichotomous) variables must
- be coded as ordered factors. NOTE: data input is different in
+ be coded as ordered factors. Each level of these ordered factors must
+ be present in the data passed to the function. NOTE: data input is different in
\code{MCMCmixfactanal} than in either \code{MCMCfactanal} or
\code{MCMCordfactanal}.}
@@ -207,7 +208,10 @@ MCMCmixfactanal(x, factors, lambda.constraints=list(),
efficiency. Please consult the coda documentation for a comprehensive
list of functions that can be used to analyze the posterior density
sample.
- }
+
+ As is the case with all measurement models, make sure that you have plenty
+ of free memory, especially when storing the scores.
+
}
\references{
@@ -232,6 +236,20 @@ MCMCmixfactanal(x, factors, lambda.constraints=list(),
\examples{
\dontrun{
+data(PErisk)
+
+post <- MCMCmixfactanal(~courts+barb2+prsexp2+prscorr2+gdpw2,
+ factors=1, data=PErisk,
+ lambda.constraints = list(courts=list(2,"-")),
+ burnin=5000, mcmc=1000000, thin=50,
+ verbose=TRUE, L0=.25, store.lambda=TRUE,
+ store.scores=TRUE, tune=1.2)
+plot(post)
+summary(post)
+
+
+
+
library(MASS)
data(Cars93)
attach(Cars93)
@@ -259,7 +277,7 @@ new.cars$Origin <- ordered(new.cars$Origin)
-posterior <- MCMCmixfactanal(~log.Price+log.MPG.city+
+post <- MCMCmixfactanal(~log.Price+log.MPG.city+
log.MPG.highway+Cylinders+log.EngineSize+
log.Horsepower+RPM+Length+
Wheelbase+Width+Weight+Origin, data=new.cars,
@@ -267,8 +285,9 @@ posterior <- MCMCmixfactanal(~log.Price+log.MPG.city+
log.Horsepower=c(3,0), weight=list(3,"+")),
factors=2,
burnin=5000, mcmc=500000, thin=100, verbose=TRUE,
- L0=.25, tune=1.5)
-
+ L0=.25, tune=3.0)
+plot(post)
+summary(post)
}
}
diff --git a/man/MCMCmnl.Rd b/man/MCMCmnl.Rd
new file mode 100644
index 0000000..5f03fea
--- /dev/null
+++ b/man/MCMCmnl.Rd
@@ -0,0 +1,266 @@
+\name{MCMCmnl}
+\alias{MCMCmnl}
+
+\title{Markov Chain Monte Carlo for Multinomial Logistic Regression}
+\description{
+ This function generates a posterior density sample
+ from a multinomial logistic regression model using either a random
+ walk Metropolis algorithm or a slice sampler. The user supplies data
+ and priors, and a sample from the posterior density is returned as an
+ mcmc object, which can be subsequently analyzed with functions provided
+ in the coda package.
+}
+
+\usage{
+MCMCmnl(formula, baseline = NULL, data = parent.frame(),
+ burnin = 1000, mcmc = 10000, thin = 1,
+ mcmc.method = "MH", tune = 1.1, verbose = FALSE,
+ seed = NA, beta.start = NA, b0 = 0, B0 = 0, ...) }
+
+\arguments{
+ \item{formula}{Model formula.
+
+ If the choicesets do not vary across individuals,
+ the \code{y} variable should be a
+ factor or numeric variable that gives the observed choice of each
+ individual. If the choicesets do vary across individuals, \code{y} should be
+ a \eqn{n \times p}{n x p} matrix where \eqn{n}{n} is the number of
+ individuals and
+ \eqn{p}{p} is the maximum number of choices in any choiceset.
+ Here each column
+ of \code{y} corresponds to a particular observed choice and the
+ elements of \code{y}
+ should be either \code{0} (not chosen but available), \code{1}
+ (chosen),
+ or \code{-999} (not available).
+
+ Choice-specific covariates have to be entered using the syntax:
+ \code{choicevar(cvar, "var", "choice")} where \code{cvar} is the name
+ of a variable in \code{data}, \code{"var"} is the name of the new
+ variable to be created, and \code{"choice"} is the level of \code{y}
+ that \code{cvar} corresponds to. Specifying each choice-specific
+ covariate will typically require \eqn{p}{p} calls to the
+ \code{choicevar} function in the formula.
+
+ Individual specific covariates can be entered into the formula
+ normally.
+
+ See the examples section below to see the syntax used to fit
+ various models.}
+
+ \item{baseline}{The baseline category of the response
+ variable. \code{baseline} should be set equal to one of the observed
+ levels of the response variable. If left equal to \code{NULL} then the
+ baseline level is set to the alpha-numerically first element of the
+ response variable. If the choicesets vary across individuals, the
+ baseline choice must be in the choiceset of each individual.
+ }
+ \item{data}{The data frame used for the analysis. Each row of the
+ dataframe should correspond to an individual who is making a
+ choice. }
+
+ \item{burnin}{The number of burn-in iterations for the sampler.}
+
+ \item{mcmc}{The number of iterations to run the sampler past burn-in. }
+
+ \item{thin}{The thinning interval used in the simulation. The number of
+ mcmc iterations must be divisible by this value. }
+
+ \item{mcmc.method}{Can be set to either "MH" (default) or "slice" to
+ perform random walk Metropolis sampling or slice sampling
+ respectively.}
+
+ \item{tune}{Metropolis tuning parameter. Can be either a positive
+ scalar or a \eqn{k}{k}-vector, where \eqn{k}{k} is the length of
+ \eqn{\beta}{beta}. Make sure that the
+ acceptance rate is satisfactory (typically between 0.20 and 0.5)
+ before using the posterior density sample for inference. }
+
+
+ \item{verbose}{A switch which determines whether or not the progress of
+ the sampler is printed to the screen. If TRUE, the iteration number,
+ the current beta vector, and the Metropolis acceptance rate are
+ printed to the screen every 500 iterations. }
+
+ \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{beta.start}{The starting value for the \eqn{\beta}{beta} vector.
+ This can either be a scalar or a column vector with dimension equal
+ to the number of betas. If this takes a scalar value, then that
+ value will serve as the starting value for all of the betas. The
+ default value of NA will use the maximum likelihood estimate of
+ \eqn{\beta}{beta} as the starting value. }
+
+ \item{b0}{The prior mean of \eqn{\beta}{beta}. This can either be a
+ scalar or a column vector with dimension equal to the number of
+ betas. If this takes a scalar value, then that value will serve as
+ the prior mean for all of the betas. }
+
+ \item{B0}{The prior precision of \eqn{\beta}{beta}. This can either be
+ a scalar or a square matrix with dimensions equal to the number of
+ betas. If this takes a scalar value, then that value times an
+ identity matrix serves as the prior precision of
+ \eqn{\beta}{beta}. Default value of 0 is equivalent to an improper
+ uniform prior for beta.}
+
+ \item{\dots}{Further arguments to be passed. }
+}
+\value{
+ An mcmc object that contains the posterior density sample. This
+ object can be summarized by functions provided by the coda package.
+}
+\details{
+\code{MCMCmnl} simulates from the posterior density of a multinomial logistic
+ regression model using either a random walk Metropolis algorithm or a
+ univariate slice sampler. The simulation
+ proper is done in compiled C++ code to maximize efficiency. Please consult
+ the coda documentation for a comprehensive list of functions that can be
+ used to analyze the posterior density sample.
+
+ The model takes the following form: \deqn{y_i \sim
+ \mathcal{M}ultinomial(\pi_i)}{y_i ~ Multinomial(pi_i)}
+
+ where:
+ \deqn{\pi_{ij} =
+ \frac{\exp(x_{ij}'\beta)}{\sum_{k=1}^p\exp(x_{ik}'\beta)}}{pi_{ij} =
+ exp(x_{ij}'beta) / [sum_{k=1}^p exp(x_{ik}'beta)]}
+
+ We assume a multivariate Normal prior on \eqn{\beta}{beta}:
+ \deqn{\beta \sim \mathcal{N}(b_0,B_0^{-1})}{beta ~ N(b0,B0^(-1))}
+ The Metropollis proposal distribution is centered at the current value of
+ \eqn{\beta}{beta} and has variance-covariance \eqn{V = T
+ (B_0 + C^{-1})^{-1} T }{V = T (B0 + C^{-1})^{-1} T}, where
+ \eqn{T}{T} is a the diagonal positive definite matrix formed from the
+ \code{tune}, \eqn{B_0}{B0} is the prior precision, and \eqn{C}{C} is
+ the large sample variance-covariance matrix of the MLEs. This last
+ calculation is done via an initial call to \code{optim}.
+}
+
+\references{
+
+ 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/}.
+}
+
+\seealso{\code{\link[coda]{plot.mcmc}},\code{\link[coda]{summary.mcmc}},
+\code{\link[nnet]{multinom}}}
+\examples{
+ \dontrun{
+ data(Nethvote)
+
+ ## just a choice-specific X var
+ post1 <- MCMCmnl(vote ~
+ choicevar(distD66, "sqdist", "D66") +
+ choicevar(distPvdA, "sqdist", "PvdA") +
+ choicevar(distVVD, "sqdist", "VVD") +
+ choicevar(distCDA, "sqdist", "CDA"),
+ baseline="D66", mcmc.method="MH", B0=0,
+ verbose=TRUE, mcmc=100000, thin=10, tune=1.0,
+ data=Nethvote)
+
+ plot(post1)
+ summary(post1)
+
+
+
+ ## just individual-specific X vars
+ post2<- MCMCmnl(vote ~
+ relig + class + income + educ + age + urban,
+ baseline="D66", mcmc.method="MH", B0=0,
+ verbose=TRUE, mcmc=100000, thin=10, tune=0.5,
+ data=Nethvote)
+
+ plot(post2)
+ summary(post2)
+
+
+
+ ## both choice-specific and individual-specific X vars
+ post3 <- MCMCmnl(vote ~
+ choicevar(distD66, "sqdist", "D66") +
+ choicevar(distPvdA, "sqdist", "PvdA") +
+ choicevar(distVVD, "sqdist", "VVD") +
+ choicevar(distCDA, "sqdist", "CDA") +
+ relig + class + income + educ + age + urban,
+ baseline="D66", mcmc.method="MH", B0=0,
+ verbose=TRUE, mcmc=100000, thin=10, tune=0.5,
+ data=Nethvote)
+
+ plot(post3)
+ summary(post3)
+
+
+ ## numeric y variable
+ nethvote$vote <- as.numeric(nethvote$vote)
+ post4 <- MCMCmnl(vote ~
+ choicevar(distD66, "sqdist", "2") +
+ choicevar(distPvdA, "sqdist", "3") +
+ choicevar(distVVD, "sqdist", "4") +
+ choicevar(distCDA, "sqdist", "1") +
+ relig + class + income + educ + age + urban,
+ baseline="2", mcmc.method="MH", B0=0,
+ verbose=TRUE, mcmc=100000, thin=10, tune=0.5,
+ data=Nethvote)
+
+
+ plot(post4)
+ summary(post4)
+
+
+
+ ## Simulated data example with nonconstant choiceset
+ n <- 1000
+ y <- matrix(0, n, 4)
+ colnames(y) <- c("a", "b", "c", "d")
+ xa <- rnorm(n)
+ xb <- rnorm(n)
+ xc <- rnorm(n)
+ xd <- rnorm(n)
+ xchoice <- cbind(xa, xb, xc, xd)
+ z <- rnorm(n)
+ for (i in 1:n){
+ ## randomly determine choiceset (c is always in choiceset)
+ choiceset <- c(3, sample(c(1,2,4), 2, replace=FALSE))
+ numer <- matrix(0, 4, 1)
+ for (j in choiceset){
+ if (j == 3){
+ numer[j] <- exp(xchoice[i, j] )
+ }
+ else {
+ numer[j] <- exp(xchoice[i, j] - z[i] )
+ }
+ }
+ p <- numer / sum(numer)
+ y[i,] <- rmultinom(1, 1, p)
+ y[i,-choiceset] <- -999
+ }
+
+ post5 <- MCMCmnl(y~choicevar(xa, "x", "a") +
+ choicevar(xb, "x", "b") +
+ choicevar(xc, "x", "c") +
+ choicevar(xd, "x", "d") + z,
+ baseline="c", verbose=TRUE,
+ mcmc=100000, thin=10, tune=.85)
+
+ plot(post5)
+ summary(post5)
+
+ }
+}
+\keyword{models}
+
diff --git a/man/MCMCoprobit.Rd b/man/MCMCoprobit.Rd
index 9c375da..4d13702 100644
--- a/man/MCMCoprobit.Rd
+++ b/man/MCMCoprobit.Rd
@@ -1,6 +1,6 @@
\name{MCMCoprobit}
\alias{MCMCoprobit}
-\title{Markov chain Monte Carlo for Ordered Probit Regression}
+\title{Markov Chain Monte Carlo for Ordered Probit Regression}
\description{
This function generates a posterior density sample
from an ordered probit regression model using the data augmentation
diff --git a/man/MCMCordfactanal.Rd b/man/MCMCordfactanal.Rd
index a590921..67e9cb8 100644
--- a/man/MCMCordfactanal.Rd
+++ b/man/MCMCordfactanal.Rd
@@ -1,6 +1,6 @@
\name{MCMCordfactanal}
\alias{MCMCordfactanal}
-\title{Markov chain Monte Carlo for Ordinal Data Factor Analysis Model}
+\title{Markov Chain Monte Carlo for Ordinal Data Factor Analysis Model}
\description{
This function generates a posterior density sample from an ordinal data
factor analysis model. Normal priors are assumed on the factor
@@ -13,7 +13,7 @@
\usage{
MCMCordfactanal(x, factors, lambda.constraints=list(),
- data=parent.environment(), burnin = 1000, mcmc = 20000,
+ data=parent.frame(), burnin = 1000, mcmc = 20000,
thin=1, tune=NA, verbose = FALSE, seed = NA,
lambda.start = NA, l0=0, L0=0,
store.lambda=TRUE, store.scores=FALSE,
@@ -174,8 +174,10 @@ MCMCordfactanal(x, factors, lambda.constraints=list(),
efficiency. Please consult the coda documentation for a comprehensive
list of functions that can be used to analyze the posterior density
sample.
+
+ As is the case with all measurement models, make sure that you have plenty
+ of free memory, especially when storing the scores.
}
-}
\references{
Shawn Treier and Simon Jackman. 2003. ``Democracy as a Latent Variable."
diff --git a/man/MCMCpanel.Rd b/man/MCMCpanel.Rd
index 555a2a6..47781fa 100644
--- a/man/MCMCpanel.Rd
+++ b/man/MCMCpanel.Rd
@@ -1,6 +1,6 @@
\name{MCMCpanel}
\alias{MCMCpanel}
-\title{Markov chain Monte Carlo for the General Linear Panel Model}
+\title{Markov Chain Monte Carlo for the General Linear Panel Model}
\description{
MCMCpanel generates a posterior density sample from a General
Linear Panel Model using Algorithm 2 of Chib and Carlin (1999).
diff --git a/man/MCMCpoisson.Rd b/man/MCMCpoisson.Rd
index bbaca67..dc0a2fe 100644
--- a/man/MCMCpoisson.Rd
+++ b/man/MCMCpoisson.Rd
@@ -1,6 +1,6 @@
\name{MCMCpoisson}
\alias{MCMCpoisson}
-\title{Markov chain Monte Carlo for Poisson Regression}
+\title{Markov Chain Monte Carlo for Poisson Regression}
\description{
This function generates a posterior density sample
from a Poisson regression model using a random walk Metropolis
diff --git a/man/MCMCprobit.Rd b/man/MCMCprobit.Rd
index f470d2e..2ed7219 100644
--- a/man/MCMCprobit.Rd
+++ b/man/MCMCprobit.Rd
@@ -1,6 +1,6 @@
\name{MCMCprobit}
\alias{MCMCprobit}
-\title{Markov chain Monte Carlo for Probit Regression}
+\title{Markov Chain Monte Carlo for Probit Regression}
\description{
This function generates a posterior density sample
from a probit regression model using the data augmentation
diff --git a/man/MCMCtobit.Rd b/man/MCMCtobit.Rd
new file mode 100644
index 0000000..100d9dd
--- /dev/null
+++ b/man/MCMCtobit.Rd
@@ -0,0 +1,155 @@
+\name{MCMCtobit}
+\alias{MCMCtobit}
+\title{Markov Chain Monte Carlo for Gaussian Linear Regression with a Censored Dependent Variable}
+\description{
+ This function generates a posterior density sample
+ from a linear regression model with Gaussian errors using
+ Gibbs sampling (with a multivariate Gaussian prior on the
+ beta vector, and an inverse Gamma prior on the conditional
+ error variance). The dependent variable may be censored
+ from below, from above, or both. The user supplies data
+ and priors, and a sample from the posterior density is
+ returned as an mcmc object, which can be subsequently
+ analyzed with functions provided in the coda package.
+ }
+\usage{
+MCMCtobit(formula, data = parent.frame(), below = 0, above = Inf,
+ burnin = 1000, mcmc = 10000, thin = 1, verbose = FALSE, seed = NA,
+ beta.start = NA, b0 = 0, B0 = 0, c0 = 0.001, d0 = 0.001, ...)
+}
+\arguments{
+ \item{formula}{A model formula.}
+
+ \item{data}{A dataframe.}
+
+ \item{below}{The point at which the dependent variable is
+ censored from below. The default is zero. To censor from
+ above only, specify that below = -Inf.}
+
+ \item{above}{The point at which the dependent variable is
+ censored from above. To censor from below only, use the
+ default value of Inf.}
+
+ \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 TRUE, the iteration number, the
+ \eqn{\beta}{beta} vector, and the conditional error variance is printed to
+ the screen every 1000 iterations.}
+
+ \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{beta.start}{The starting values for the \eqn{\beta}{beta} vector.
+ This can either be a scalar or a column vector with dimension equal to the
+ number of betas. The default value of of NA will use the OLS estimate of
+ \eqn{\beta}{beta} as the starting value. If this is a scalar, that value
+ will serve as the starting value mean for all of the betas.}
+
+ \item{b0}{The prior mean of \eqn{\beta}{beta}. This can either be a
+ scalar or a column vector with dimension equal to the number of betas. If
+ this takes a scalar value, then that value will serve as the prior
+ mean for all of the betas.}
+
+ \item{B0}{The prior precision of \eqn{\beta}{beta}. This can either be a
+ scalar or a square matrix with dimensions equal to the number of
+ betas.
+ If this takes a scalar value, then that value times an
+ identity matrix serves as the prior precision of beta. Default value
+ of 0 is equivalent to an improper uniform prior for beta.}
+
+ \item{c0}{\eqn{c_0/2}{c0/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{c_0}{c0} pseudo-observations.}
+
+ \item{d0}{\eqn{d_0/2}{d0/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{d_0}{d0} acts like the sum of squared errors from the
+ \eqn{c_0}{c0} pseudo-observations.}
+
+ \item{...}{further arguments to be passed}
+
+}
+
+\details{
+ \code{MCMCtobit} simulates from the posterior density using standard
+ Gibbs sampling (a multivariate Normal draw for the betas, and an
+ inverse Gamma draw for the conditional error variance). \code{MCMCtobit}
+ differs from \code{MCMCregress} in that the dependent variable may be
+ censored from below, from above, or both. The simulation
+ proper is done in compiled C++ code to maximize efficiency. Please consult
+ the coda documentation for a comprehensive list of functions that can be
+ used to analyze the posterior density sample.
+
+ The model takes the following form:
+ \deqn{y_i = x_i ' \beta + \varepsilon_{i},}{y_i = x_i'beta + epsilon_i,}
+ where the errors are assumed to be Gaussian:
+ \deqn{\varepsilon_{i} \sim \mathcal{N}(0, \sigma^2).}{epsilon_i ~ N(0,
+ sigma^2).}
+ Let \eqn{c_1} and \eqn{c_2} be the two censoring points, and let
+ \eqn{y_i^\ast}{y_i^star} be the partially observed dependent variable. Then,
+ \deqn{y_i = y_i^{\ast} \texttt{ if } c_1 < y_i^{\ast} < c_2,}{y_i = y_i^star if c_1 < y_i^star < c_2,}
+ \deqn{y_i = c_1 \texttt{ if } c_1 \geq y_i^{\ast},}{y_i = c_1 if c_1 >= y_i^star,}
+ \deqn{y_i = c_2 \texttt{ if } c_2 \leq y_i^{\ast}.}{y_i = c_2 if c_1 <= y_i^star.}
+ We assume standard, semi-conjugate priors:
+ \deqn{\beta \sim \mathcal{N}(b_0,B_0^{-1}),}{beta ~ N(b0,B0^(-1)),}
+ and:
+ \deqn{\sigma^{-2} \sim \mathcal{G}amma(c_0/2, d_0/2),}{sigma^(-2) ~
+ Gamma(c0/2, d0/2),}
+ where \eqn{\beta}{beta} and \eqn{\sigma^{-2}}{sigma^(-2)} are assumed
+ \emph{a priori} independent. Note that only starting values for
+ \eqn{\beta}{beta} are allowed because simulation is done using
+ Gibbs sampling with the conditional error variance
+ as the first block in the sampler.
+ }
+
+\value{
+ An mcmc object that contains the posterior density sample. This
+ object can be summarized by functions provided by the coda package.
+}
+
+\references{
+ 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/}.
+
+ James Tobin. 1958. ``Estimation of relationships for limited dependent
+ variables." \emph{Econometrica.}, 26:24-36.
+}
+
+\author{Ben Goodrich}
+
+\seealso{
+ \code{\link[coda]{plot.mcmc}},
+ \code{\link[coda]{summary.mcmc}},
+ \code{\link[survival]{survreg}},
+ \code{\link[MCMCpack]{MCMCregress}}}
+
+\examples{
+library(survival)
+example(tobin)
+summary(tfit)
+tfit.mcmc <- MCMCtobit(durable ~ age + quant, data=tobin, mcmc=30000, verbose=TRUE)
+plot(tfit.mcmc)
+raftery.diag(tfit.mcmc)
+summary(tfit.mcmc)
+}
+
+\keyword{models}
diff --git a/man/Nethvote.Rd b/man/Nethvote.Rd
new file mode 100644
index 0000000..9df10be
--- /dev/null
+++ b/man/Nethvote.Rd
@@ -0,0 +1,56 @@
+\name{Nethvote}
+\alias{Nethvote}
+\docType{data}
+\title{Dutch Voting Behavior in 1989}
+\description{
+Dutch Voting Behavior in 1989.
+}
+\usage{data(Nethvote)}
+\format{
+ A data frame with 1754 observations and 11 variables from the 1989
+ Dutch Parliamentary Election Study
+ (Anker and Oppenhuis, 1993). Each observation is a survey respondent.
+ These data are a subset of one of five multiply imputed datasets used in
+ Quinn and Martin (2002). For more information see Quinn and Martin
+ (2002).
+ \describe{
+ \item{vote}{A factor giving the self-reported vote choice of each
+ respondent. The levels are CDA (Christen Democratisch Appel), D66
+ (Democraten 66), Pvda (Partij van de Arbeid), and VVD (Volkspartij
+ voor Vrijheid en Democratie).}
+ \item{distD66}{A numeric variable giving the squared ideological distance
+ between the respondent and the D66. Larger values indicate ideological
+ dissimilarity between the respondent and the party.}
+ \item{distPvdA}{A numeric variable giving the squared ideological distance
+ between the respondent and the PvdA. Larger values indicate ideological
+ dissimilarity between the respondent and the party.}
+ \item{distVVD}{A numeric variable giving the squared ideological distance
+ between the respondent and the VVD. Larger values indicate ideological
+ dissimilarity between the respondent and the party.}
+ \item{distCDA}{A numeric variable giving the squared ideological distance
+ between the respondent and the CDA. Larger values indicate ideological
+ dissimilarity between the respondent and the party.}
+ \item{relig}{An indicator variable equal to 0 if the respondent is
+ not religious and 1 if the respondent is religious.}
+ \item{class}{Social class of respondent. 0 is the lowest social class, 4 is
+ the highest social class.}
+ \item{income}{Income of respondent. 0 is lowest and 6 is highest.}
+ \item{educ}{Education of respondent. 0 is lowest and 4 is highest.}
+ \item{age}{Age category of respondent. 0 is lowest and 12 is
+ highest.}
+ \item{urban}{Indicator variable equal to 0 if the respondent is not
+ a resident of an urban area and 1 if the respondent is a resident
+ of an urban area.}
+ }
+}
+\source{
+ H. Anker and E.V. Oppenhuis. 1993. ``Dutch Parliamentary Election
+ Study.'' (computer file). Dutch Electoral Research Foundation and
+ Netherlands Central Bureau of Statistics, Amsterdam.
+}
+\references{
+ K.M. Quinn and A.M. Martin. 2002. ``An Integrated Computational Model
+ of Multiparty Electoral Competition.'' \emph{Statistical Science}. 17:
+ 405-419.
+}
+\keyword{datasets}
diff --git a/man/PErisk.Rd b/man/PErisk.Rd
index df984d6..7027b91 100644
--- a/man/PErisk.Rd
+++ b/man/PErisk.Rd
@@ -5,7 +5,7 @@
\description{
Political Economic Risk Data from 62 Countries in 1987.
}
-\usage{data(CountryRisk)}
+\usage{data(PErisk)}
\format{
A data frame with 62 observations on the following 9 variables. All
data points are from 1987. See Quinn (2004) for more details.
@@ -36,8 +36,8 @@ Political Economic Risk Data from 62 Countries in 1987.
\url{http://www.ssc.upenn.edu/~cheibub/data/}.
Witold J. Henisz. 2002. ``The Political Constraint Index (POLCON)
- Dataset.''
- \url{http://www-management.wharton.upenn.edu/henisz/POLCON/ContactInfo.html}.
+ Dataset.'' \\
+ \url{http://www-management.wharton.upenn.edu/henisz/POLCON/ContactInfo.html}.
Monty G. Marshall, Ted Robert Gurr, and Barbara Harff. 2002. ``State
Failure Task Force Problem Set.''
diff --git a/man/choicevar.Rd b/man/choicevar.Rd
new file mode 100644
index 0000000..5d2df36
--- /dev/null
+++ b/man/choicevar.Rd
@@ -0,0 +1,25 @@
+\name{choicevar}
+\alias{choicevar}
+
+\title{Handle Choice-Specific Covariates in Multinomial Choice Models}
+\description{
+This function handles choice-specific covariates in multinomial choice models.
+See the example for an example of useage.
+}
+\usage{
+choicevar(var, varname, choicelevel)
+}
+
+\arguments{
+ \item{var}{The is the name of the variable in the dataframe.}
+ \item{varname}{The name of the new variable to be created.}
+ \item{choicelevel}{The level of \code{y} that the variable
+ corresponds to.}
+}
+
+\value{
+ The new variable used by the \code{MCMCmnl()} function.
+}
+
+\seealso{\code{\link{MCMCmnl}}}
+\keyword{manip}
diff --git a/man/tomog.Rd b/man/tomog.Rd
index 8a815c8..08b621f 100644
--- a/man/tomog.Rd
+++ b/man/tomog.Rd
@@ -74,7 +74,7 @@ c1 <- (r0 + r1) - c0
tomogplot(r0, r1, c0, c1)
}
-\seealso{\code{\link{MCMCbaselineEI}}, \code{\link{MCMChierEI}},
+\seealso{\code{\link{MCMChierEI}},
\code{\link{MCMCdynamicEI}}, \code{\link{dtomogplot}}
}
diff --git a/man/wishart.Rd b/man/wishart.Rd
index 5f4b274..f931d90 100644
--- a/man/wishart.Rd
+++ b/man/wishart.Rd
@@ -15,13 +15,15 @@
\arguments{
\item{W}{Positive definite matrix W \eqn{(p \times p)}{(p x p)}.}
\item{v}{Wishart degrees of freedom (scalar).}
- \item{S}{Wishart scale matrix \eqn{(p \times p)}{(p x p)}.}}
+ \item{S}{Wishart scale matrix \eqn{(p \times p)}{(p x p)}.}
+}
\value{
\code{dwish} evaluates the density at positive definite matrix W.
\code{rwish} generates one random draw from the distribution.
}
+
\examples{
density <- dwish(matrix(c(2,-.3,-.3,4),2,2), 3, matrix(c(1,.3,.3,1),2,2))
draw <- rwish(3, matrix(c(1,.3,.3,1),2,2))
diff --git a/src/MCMCdistn.cc b/src/MCMCdistn.cc
new file mode 100644
index 0000000..062938b
--- /dev/null
+++ b/src/MCMCdistn.cc
@@ -0,0 +1,466 @@
+// MCMCdistn.cc contains a number of functions that interface the scythe
+// random number generators, density functions, and distribution functions
+// with the MCMCpack R package so that these functions can be called
+// directly from R.
+//
+// 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
+//
+
+#include "matrix.h"
+#include "distributions.h"
+#include "stat.h"
+#include "la.h"
+#include "ide.h"
+#include "smath.h"
+#include "MCMCrng.h"
+
+
+#include <R.h> // needed to use Rprintf()
+#include <R_ext/Utils.h> // needed to allow user interrupts
+
+using namespace SCYTHE;
+using namespace std;
+
+extern "C"{
+
+ void rtnormFromR(double* sampledata, const int* samplesize,
+ const double* mean, const int* meanlength,
+ const double* var, const int* varlength,
+ const double* below, const int* belowlength,
+ const double* above, const int* abovelength,
+ const int *lecuyer, const int *seedarray,
+ const int *lecuyerstream){
+
+ // initialize rng stream
+ rng *stream = MCMCpack_get_rng(*lecuyer, seedarray, *lecuyerstream);
+
+ // individual indexing to do R-style index recycling
+ int m_index = -1;
+ int v_index = -1;
+ int b_index = -1;
+ int a_index = -1;
+
+ // the loop where the random number draws take place
+ for (int i=0; i<*samplesize; ++i){
+
+ // do the R-style index recycling
+ ++m_index;
+ if (m_index==*meanlength){
+ m_index = 0;
+ }
+ ++v_index;
+ if (v_index==*varlength){
+ v_index = 0;
+ }
+ ++b_index;
+ if (b_index==*belowlength){
+ b_index = 0;
+ }
+ ++a_index;
+ if (a_index==*abovelength){
+ a_index = 0;
+ }
+
+ sampledata[i] = stream->rtnorm(mean[m_index], var[v_index],
+ below[b_index], above[a_index]);
+
+ }
+
+ } // end rtnormFromR
+
+ //This next function implements the algorithm for rnchypgeom for a
+ //list of data by calling double rng::rnchypgeom() for each element.
+ void rnchypgeomFromR(const int * number,
+ const double* m1, const double* n1,
+ const double* n2, const double* psi,
+ const double* delta, const int * listLengths,
+ double* result, const int* lecuyer,
+ const int* seedarray, const int* lecuyerstream) {
+
+ //First, we initialize our rng stream.
+ rng * stream = MCMCpack_get_rng(*lecuyer, seedarray, *lecuyerstream);
+
+ //Second, we create some indices for wrap-around. Note that we will
+ //be assuming that the listLenths array has five elements, which it
+ //should if it was called from the R code.
+ int count[5]={-1,-1,-1,-1,-1};
+
+ //Third, we use the function to calculate the return values, while
+ //being sure to use R-like vector wraparound for all of our parameters.
+ for(int i=0; i < *number; ++i) {
+
+ for(int j=0; j < 5; ++j) {
+ count[j] = (count[j] + 1) % listLengths[j];
+ }
+
+ result[i] = stream->rnchypgeom(m1[count[0]], n1[count[1]], n2[count[2]],
+ psi[count[3]], delta[count[4]]);
+ }
+
+ //Finally, we're done, so we just end.
+ }
+
+ //This next function fills the given list from rng::richisq() from the
+ //given value of the nu (degrees of freedom).
+ void richisqFromR(const int * number, const double * nu,
+ const int * listLength, double * result,
+ const int *lecuyer, const int *seedarray,
+ const int *lecuyerstream) {
+
+ //First, we initialize our rng stream.
+ rng * stream = MCMCpack_get_rng(*lecuyer, seedarray, *lecuyerstream);
+
+ //We create an index for the wrap-around.
+ int count = -1;
+
+ //We fill our array now from our function.
+ for(int i = 0; i < *number; ++i) {
+ count = (count + 1)%(*listLength);
+
+ result[i] = stream->richisq(nu[count]);
+ }
+
+ //Now, we're done, so get out of here.
+ }
+
+ //This next function fills the given list from rng::rigamma() from the
+ //given value of the alpha and beta.
+ void rigammaFromR(const int * number, const double * alpha,
+ const double * beta, const int * listLengths,
+ double * result, const int *lecuyer,
+ const int *seedarray, const int *lecuyerstream) {
+
+ //First, we initialize our rng stream.
+ rng * stream = MCMCpack_get_rng(*lecuyer, seedarray, *lecuyerstream);
+
+ //We create an index for the wrap-around.
+ int count[2] = {-1,-1};
+
+ //We fill our array now from our function.
+ for(int i = 0; i < *number; ++i) {
+
+ for(int j = 0; j < 2; ++j) {
+ count[j] = (count[j] + 1)%(listLengths[j]);
+ }
+
+ result[i] = stream->rigamma(alpha[count[0]], beta[count[1]]);
+ }
+
+ //Now, we're done, so get out of here.
+ }
+
+ //This next function fills the given list from rng::rbern() from the
+ //given value of the p parameter.
+ void rbernFromR(const int * number, const double * p,
+ const int * listLength, double * result,
+ const int *lecuyer, const int *seedarray,
+ const int *lecuyerstream) {
+
+ //First, we initialize our rng stream.
+ rng * stream = MCMCpack_get_rng(*lecuyer, seedarray, *lecuyerstream);
+
+ //We create an index for the wrap-around.
+ int count = -1;
+
+ //We fill our array now from our function.
+ for(int i = 0; i < *number; ++i) {
+ count = (count + 1)%(*listLength);
+
+ result[i] = stream->rbern(p[count]);
+ }
+
+ //Now, we're done, so get out of here.
+ }
+
+ //This next function implements the algorithm for rtnormcombo for a
+ //list of data by calling double rng::rtnorm_combo() for each element.
+ void rtnormcomboFromR(const int * number,
+ const double* m, const double* v,
+ const double* below, const double* above,
+ const int * listLengths,
+ double* result, const int* lecuyer,
+ const int* seedarray, const int* lecuyerstream) {
+
+ //First, we initialize our rng stream.
+ rng * stream = MCMCpack_get_rng(*lecuyer, seedarray, *lecuyerstream);
+
+ //Second, we create some indices for wrap-around. Note that we will
+ //be assuming that the listLenths array has four elements, which it
+ //should if it was called from the R code.
+ int count[4]={-1,-1,-1,-1};
+
+ //Third, we use the function to calculate the return values, while
+ //being sure to use R-like vector wraparound for all of our parameters.
+ for(int i=0; i < *number; ++i) {
+
+ for(int j=0; j < 4; ++j) {
+ count[j] = (count[j] + 1) % listLengths[j];
+ }
+
+ result[i] = stream->rtnorm_combo(m[count[0]], v[count[1]], below[count[2]],
+ above[count[3]]);
+ }
+
+ //Finally, we're done, so we just end.
+ }
+
+ //This next function implements the algorithm for rtbnormslice for a
+ //list of data by calling double rng::rtbnorm_slice() for each element.
+ void rtbnormsliceFromR(const int * number,
+ const double* m, const double* v,
+ const double* below, const int* iter,
+ const int * listLengths,
+ double* result, const int* lecuyer,
+ const int* seedarray, const int* lecuyerstream) {
+
+ //First, we initialize our rng stream.
+ rng * stream = MCMCpack_get_rng(*lecuyer, seedarray, *lecuyerstream);
+
+ //Second, we create some indices for wrap-around. Note that we will
+ //be assuming that the listLenths array has four elements, which it
+ //should if it was called from the R code.
+ int count[4]={-1,-1,-1,-1};
+
+ //Third, we use the function to calculate the return values, while
+ //being sure to use R-like vector wraparound for all of our parameters.
+ for(int i=0; i < *number; ++i) {
+
+ for(int j=0; j < 4; ++j) {
+ count[j] = (count[j] + 1) % listLengths[j];
+ }
+
+ result[i] = stream->rtbnorm_slice(m[count[0]], v[count[1]], below[count[2]],
+ iter[count[3]]);
+ }
+
+ //Finally, we're done, so we just end.
+ }
+
+ //This next function implements the algorithm for rtanormslice for a
+ //list of data by calling double rng::rtanorm_slice() for each element.
+ void rtanormsliceFromR(const int * number,
+ const double* m, const double* v,
+ const double* above, const int* iter,
+ const int * listLengths,
+ double* result, const int* lecuyer,
+ const int* seedarray, const int* lecuyerstream) {
+
+ //First, we initialize our rng stream.
+ rng * stream = MCMCpack_get_rng(*lecuyer, seedarray, *lecuyerstream);
+
+ //Second, we create some indices for wrap-around. Note that we will
+ //be assuming that the listLenths array has four elements, which it
+ //should if it was called from the R code.
+ int count[4]={-1,-1,-1,-1};
+
+ //Third, we use the function to calculate the return values, while
+ //being sure to use R-like vector wraparound for all of our parameters.
+ for(int i=0; i < *number; ++i) {
+
+ for(int j=0; j < 4; ++j) {
+ count[j] = (count[j] + 1) % listLengths[j];
+ }
+
+ result[i] = stream->rtanorm_slice(m[count[0]], v[count[1]], above[count[2]],
+ iter[count[3]]);
+ }
+
+ //Finally, we're done, so we just end.
+ }
+
+ //This next function implements the algorithm for rtbnormcombo for a
+ //list of data by calling double rng::rtbnorm_combo() for each element.
+ void rtbnormcomboFromR(const int * number,
+ const double* m, const double* v,
+ const double* below, const int* iter,
+ const int * listLengths,
+ double* result, const int* lecuyer,
+ const int* seedarray, const int* lecuyerstream) {
+
+ //First, we initialize our rng stream.
+ rng * stream = MCMCpack_get_rng(*lecuyer, seedarray, *lecuyerstream);
+
+ //Second, we create some indices for wrap-around. Note that we will
+ //be assuming that the listLenths array has four elements, which it
+ //should if it was called from the R code.
+ int count[4]={-1,-1,-1,-1};
+
+ //Third, we use the function to calculate the return values, while
+ //being sure to use R-like vector wraparound for all of our parameters.
+ for(int i=0; i < *number; ++i) {
+
+ for(int j=0; j < 4; ++j) {
+ count[j] = (count[j] + 1) % listLengths[j];
+ }
+
+ result[i] = stream->rtbnorm_combo(m[count[0]], v[count[1]], below[count[2]],
+ iter[count[3]]);
+ }
+
+ //Finally, we're done, so we just end.
+ }
+
+ //This next function implements the algorithm for rtanormcombo for a
+ //list of data by calling double rng::rtanorm_combo() for each element.
+ void rtanormcomboFromR(const int * number,
+ const double* m, const double* v,
+ const double* above, const int* iter,
+ const int * listLengths,
+ double* result, const int* lecuyer,
+ const int* seedarray, const int* lecuyerstream) {
+
+ //First, we initialize our rng stream.
+ rng * stream = MCMCpack_get_rng(*lecuyer, seedarray, *lecuyerstream);
+
+ //Second, we create some indices for wrap-around. Note that we will
+ //be assuming that the listLenths array has four elements, which it
+ //should if it was called from the R code.
+ int count[4]={-1,-1,-1,-1};
+
+ //Third, we use the function to calculate the return values, while
+ //being sure to use R-like vector wraparound for all of our parameters.
+ for(int i=0; i < *number; ++i) {
+
+ for(int j=0; j < 4; ++j) {
+ count[j] = (count[j] + 1) % listLengths[j];
+ }
+
+ result[i] = stream->rtanorm_combo(m[count[0]], v[count[1]], above[count[2]],
+ iter[count[3]]);
+ }
+
+ //Finally, we're done, so we just end.
+ }
+
+ //This next function implements the algorithm for rwish for a
+ //by calling double rng::rwish() for the whole matrix.
+ void rwishFromR(const int * v, const double* s,
+ const int * dimension,
+ double* result, const int* lecuyer,
+ const int* seedarray, const int* lecuyerstream) {
+
+ //First, we initialize our rng stream, being sure to store the parameters
+ //so that we can check for them the next time this is run.
+ static rng * stream = MCMCpack_get_rng(*lecuyer, seedarray, *lecuyerstream);
+ static int store_lecuyer = *lecuyer;
+ static int store_seedarray[6] = {0,0,0,0,0,0};
+ static int store_lecuyerstream = *lecuyerstream;
+ if((*lecuyer!=store_lecuyer)||(*lecuyerstream!=store_lecuyerstream)||
+ (store_seedarray[0]!=seedarray[0])||(store_seedarray[1]!=seedarray[1])||
+ (store_seedarray[2]!=seedarray[2])||(store_seedarray[3]!=seedarray[3])||
+ (store_seedarray[4]!=seedarray[4])||(store_seedarray[5]!=seedarray[5])) {
+ for(int i = 0; i < 6; ++i) store_seedarray[i] = seedarray[i];
+ store_lecuyer = *lecuyer;
+ store_lecuyerstream = *lecuyerstream;
+ delete stream;
+ stream = MCMCpack_get_rng(*lecuyer, seedarray, *lecuyerstream);
+ }
+
+ //Second, we run our rng::rwish() function to get the result
+ Matrix<double> * myMat = new Matrix<double>(*dimension,*dimension,s);
+ Matrix<double> myMatrix2 = stream->rwish(*v,*myMat);
+
+ for(int i = 0; i < (*dimension)*(*dimension); ++i) {
+ result[i] = myMatrix2[i];
+ }
+
+ //So that we don't have any memory leaks:
+ delete myMat;
+
+ //Finally, we're done, so we just end.
+ }
+
+ //This next function implements the algorithm for rdirich for a
+ //by calling double rng::rdirich() for each element.
+ void rdirichFromR(const int * n, const double* alpha,
+ const int * vsize, const int * hsize,
+ double* result, const int* lecuyer,
+ const int* seedarray, const int* lecuyerstream) {
+
+ //First, we initialize our rng stream.
+ rng * stream = MCMCpack_get_rng(*lecuyer, seedarray, *lecuyerstream);
+
+ //Second, we run our rng::rdirich() function to get the result
+ Matrix<double> alphaMatrix(*hsize,*vsize,alpha);
+ Matrix<double> temp, tempVector;
+
+ for(int i=0; i < *n; ++i) {
+ tempVector = alphaMatrix(i%(*hsize), _);
+ tempVector.resize(*vsize, 1);
+ temp = stream->rdirich(tempVector);
+ for(int j = 0; j < *vsize; ++j)
+ result[i*(*vsize)+j] = temp[j];
+ }
+
+ //Finally, we're done, so we just end.
+ }
+
+ //This next function implements the algorithm for rmvnorm for a
+ //by calling double rng::rmvnorm() for each element.
+ void rmvnormFromR(const int * n, const double* mu,
+ const int * vsize, const int * hsize,
+ const double * sigma,
+ double* result, const int* lecuyer,
+ const int* seedarray, const int* lecuyerstream) {
+
+ //First, we initialize our rng stream.
+ rng * stream = MCMCpack_get_rng(*lecuyer, seedarray, *lecuyerstream);
+
+ //Second, we run our rng::rdirich() function to get the result
+ Matrix<double> muMatrix(*hsize,*vsize,mu);
+ Matrix<double> sigmaMatrix(*vsize,*vsize,sigma);
+ Matrix<double> temp, tempVector;
+
+ for(int i=0; i < *n; ++i) {
+ tempVector = muMatrix(i%(*hsize), _);
+ tempVector.resize(*vsize, 1);
+ temp = stream->rmvnorm(tempVector, sigmaMatrix);
+ for(int j = 0; j < *vsize; ++j)
+ result[i*(*vsize)+j] = temp[j];
+ }
+
+ //Finally, we're done, so we just end.
+ }
+
+ //This next function implements the algorithm for rmvt for a
+ //by calling double rng::rmvt() for each element.
+ void rmvtFromR(const int * n, const double* sigma,
+ const int * size, const double * nu, const int * sizeNu,
+ double* result, const int* lecuyer,
+ const int* seedarray, const int* lecuyerstream) {
+
+ //First, we initialize our rng stream.
+ rng * stream = MCMCpack_get_rng(*lecuyer, seedarray, *lecuyerstream);
+
+ //Second, we run our rng::rdirich() function to get the result
+ Matrix<double> sigmaMatrix(*size,*size,sigma);
+ Matrix<double> temp;
+
+ for(int i=0; i < *n; ++i) {
+ temp = stream->rmvt(sigmaMatrix, nu[(i%(*sizeNu))]);
+ for(int j = 0; j < *size; ++j)
+ result[i*(*size)+j] = temp[j];
+ }
+
+ //Finally, we're done, so we just end.
+ }
+
+ // other *FromR functions should be added here
+
+
+
+} // end extern "C"
diff --git a/src/MCMChierEI.cc b/src/MCMChierEI.cc
index 07a6f2c..2022d7b 100644
--- a/src/MCMChierEI.cc
+++ b/src/MCMChierEI.cc
@@ -144,13 +144,13 @@ static double shrinkage(double (*logfun)(double[], const double&,
double Lbar = L;
double Rbar = R;
- int ind0;
- if (index==0){
- ind0 = 1;
- }
- else {
- ind0 = 0;
- }
+ // int ind0;
+ //if (index==0){
+ // ind0 = 1;
+ // }
+ //else {
+ // ind0 = 0;
+ // }
double theta_x1[2];
theta_x1[0] = theta[0];
theta_x1[1] = theta[1];
@@ -266,7 +266,7 @@ extern "C"{
const int p_init = 50;
const Matrix<double> widthmat(warmup_iter - warmup_burnin, 2);
- // warm up sampling to chose slice sampling parameters adaptively
+ // warm up sampling to choose slice sampling parameters adaptively
for (int iter=0; iter<warmup_iter; ++iter){
// loop over tables
diff --git a/src/MCMCirt1d.cc b/src/MCMCirt1d.cc
index 7e753bc..5d139e2 100644
--- a/src/MCMCirt1d.cc
+++ b/src/MCMCirt1d.cc
@@ -59,6 +59,8 @@ extern "C"{
const int* thetaineqcol,
const int* store
) {
+
+ using namespace SCYTHE; //Added by Matthew S. Fasman on 11/04/2004
// put together matrices
const Matrix<int> X = r2scythe(*Xrow, *Xcol, Xdata);
diff --git a/src/MCMCmixfactanal.cc b/src/MCMCmixfactanal.cc
index e5a70be..dcb7136 100644
--- a/src/MCMCmixfactanal.cc
+++ b/src/MCMCmixfactanal.cc
@@ -20,7 +20,10 @@
// revised version of older MCMCordfactanal
// 7/20/2004 KQ
// updated to new version of Scythe 7/25/2004
+// fixed a bug pointed out by Alexander Raach 1/16/2005 KQ
+//
+#include<iostream>
#include "matrix.h"
#include "distributions.h"
@@ -66,7 +69,8 @@ mixfactanalpost (double* sampledata, const int* samplerow,
const double* b0data, const int* b0row, const int* b0col,
const int* storelambda,
const int* storescores,
- int* accepts
+ int* acceptsdata, const int* acceptsrow,
+ const int* acceptscol
) {
// put together matrices
@@ -94,36 +98,8 @@ mixfactanalpost (double* sampledata, const int* samplerow,
Lampprecdata);
const Matrix <double> a0 = r2scythe(*a0row, *a0col, a0data);
const Matrix <double> b0 = r2scythe(*b0row, *b0col, b0data);
-
-
- /*
- const Matrix<double> Xstar = r2scythe(*Xrow, *Xcol, Xdata);
- const Matrix<int> X = Matrix<int>(*Xrow, *Xcol);
- for (int i=0; i<(*Xrow * *Xcol); ++i)
- X[i] = static_cast<int>(Xstar[i]);
-
- Matrix<double> Lambda = r2scythe(*Lamstartrow, *Lamstartcol,
- Lamstartdata);
- const Matrix<double> gamma = r2scythe(*gamrow, *gamcol, gamdata);
- const Matrix<double> Psi = r2scythe(*Psistartrow, *Psistartcol,
- Psistartdata);
- const Matrix<double> Psi_inv = invpd(Psi);
- const Matrix<int> ncateg = Matrix<int>(*ncatrow, *ncatcol);
- for (int i=0; i<(*ncatrow * *ncatcol); ++i)
- X[i] = static_cast<int>(ncatdata[i]);
- const Matrix<double> Lambda_eq = r2scythe(*Lameqrow, *Lameqcol,
- Lameqdata);
- const Matrix<double> Lambda_ineq = r2scythe(*Lamineqrow, *Lamineqcol,
- Lamineqdata);
- const Matrix<double> Lambda_prior_mean = r2scythe(*Lampmeanrow,
- *Lampmeancol,
- Lampmeandata);
- const Matrix<double> Lambda_prior_prec = r2scythe(*Lampprecrow,
- *Lamppreccol,
- Lampprecdata);
- const Matrix <double> a0 = r2scythe(*a0row, *a0col, a0data);
- const Matrix <double> b0 = r2scythe(*b0row, *b0col, b0data);
- */
+ Matrix<int> accepts = r2scythe(*acceptsrow, *acceptscol, acceptsdata);
+
// initialize rng stream
rng *stream = MCMCpack_get_rng(*lecuyer, seedarray, *lecuyerstream);
@@ -146,7 +122,6 @@ mixfactanalpost (double* sampledata, const int* samplerow,
// starting values for phi and gamma_p
Matrix<double> phi = Matrix<double>(N,D-1);
phi = cbind(ones<double>(N,1), phi);
- Matrix<double> gamma_p = gamma(_,0);
// storage matrices (row major order)
Matrix<double> Lambda_store;
@@ -231,17 +206,22 @@ mixfactanalpost (double* sampledata, const int* samplerow,
// sample gamma
for (int j=0; j<K; ++j){ // do the sampling for each categ. var
+ Matrix<double> gamma_p = gamma(_,j);
+ if (ncateg[j] <= 2){
+ ++accepts[j];
+ }
if (ncateg[j] > 2){
const Matrix<double> X_mean = phi * t(Lambda(j,_));
for (int i=2; i<(ncateg[j]); ++i){
if (i==(ncateg[j]-1)){
- gamma_p[i] = stream->rtbnorm_combo(gamma(i,j), ::pow(tune[j], 2.0),
- gamma_p[i-1]);
+ gamma_p[i] = stream->rtbnorm_combo(gamma(i,j),
+ ::pow(tune[j], 2.0),
+ gamma_p[i-1]);
}
else {
gamma_p[i] = stream->rtnorm_combo(gamma(i,j), ::pow(tune[j], 2.0),
- gamma_p[i-1],
- gamma(i+1, j));
+ gamma_p[i-1],
+ gamma(i+1, j));
}
}
double loglikerat = 0.0;
@@ -271,24 +251,25 @@ mixfactanalpost (double* sampledata, const int* samplerow,
}
}
}
- for (int k=2; k<(ncateg[j]-1); ++k){
+ for (int k=2; k<ncateg[j]; ++k){
loggendenrat = loggendenrat
+ log(pnorm(gamma(k+1,j), gamma(k,j), tune[j]) -
- pnorm(gamma(k-1,j), gamma(k,j), tune[j]) )
+ pnorm(gamma_p[k-1], gamma(k,j), tune[j]) )
- log(pnorm(gamma_p[k+1], gamma_p[k], tune[j]) -
- pnorm(gamma_p[k-1], gamma_p[k], tune[j]) );
+ pnorm(gamma(k-1,j), gamma_p[k], tune[j]) );
}
double logacceptrat = loglikerat + loggendenrat;
+
if (stream->runif() <= exp(logacceptrat)){
for (int i=0; i<*gamrow; ++i){
if (gamma(i,j) == 300) break;
gamma(i,j) = gamma_p[i];
}
- ++accepts[0];
+ ++accepts[j];
}
}
}
-
+
// print results to screen
if (*verbose == 1 && iter % 500 == 0){
Rprintf("\n\nMCMCmixfactanal iteration %i of %i \n", (iter+1), tot_iter);
@@ -304,9 +285,12 @@ mixfactanalpost (double* sampledata, const int* samplerow,
Rprintf("%10.5f", Psi(i,i));
}
Rprintf("\n");
- Rprintf("\nMetropolis-Hastings acceptance rate = %10.5f\n",
- static_cast<double>(*accepts)/(static_cast<double>((iter+1) *
- n_ord_ge3)));
+ Rprintf("\nMetropolis-Hastings acceptance rates = \n");
+ for (int j = 0; j<K; ++j){
+ Rprintf("%6.2f",
+ static_cast<double>(accepts[j]) /
+ static_cast<double>((iter+1)));
+ }
}
// store results
@@ -358,6 +342,10 @@ mixfactanalpost (double* sampledata, const int* samplerow,
const int size = *samplerow * *samplecol;
for (int i=0; i<size; ++i)
sampledata[i] = output[i];
+
+ for (int j=0; j<K; ++j)
+ acceptsdata[j] = accepts[j];
+
}
diff --git a/src/MCMCmnlMH.cc b/src/MCMCmnlMH.cc
new file mode 100644
index 0000000..7def301
--- /dev/null
+++ b/src/MCMCmnlMH.cc
@@ -0,0 +1,171 @@
+// MCMCmnlMH.cc samples from the posterior distribution of a multinomial
+// logit model using a random walk Metropolis algorithm.
+//
+// The initial version of this file was generated by the
+// auto.Scythe.call() function in the MCMCpack R package
+// written by:
+//
+// 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
+//
+// This file was initially generated on Wed Dec 29 15:27:08 2004
+// 12/31/2004 filled out template and got it initial version working (KQ)
+
+#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;
+
+static double mnl_logpost(const Matrix<double>& Y, const Matrix<double>& X,
+ const Matrix<double>& beta,
+ const Matrix<double>& beta_prior_mean,
+ const Matrix<double>& beta_prior_prec){
+
+ // likelihood
+ double loglike = 0.0;
+ Matrix<double> numer = exp(X * beta);
+ numer = reshape(numer, Y.rows(), Y.cols());
+ double *denom = new double[Y.rows()];
+ for (int i=0; i<Y.rows(); ++i){
+ denom[i] = 0.0;
+ for (int j=0; j<Y.cols(); ++j){
+ if (Y(i,j) != -999){
+ denom[i] += numer(i,j);
+ }
+ }
+ for (int j=0; j<Y.cols(); ++j){
+ if (Y(i,j) == 1.0){
+ loglike += std::log(numer(i,j) / denom[i]);
+ }
+ }
+ }
+
+ delete [] denom;
+
+ // prior
+ double logprior = 0.0;
+ if (beta_prior_prec(0,0) != 0){
+ logprior = lndmvn(beta, beta_prior_mean, invpd(beta_prior_prec));
+ }
+
+ return (loglike + logprior);
+}
+
+
+extern "C" {
+
+ // MCMC sampling for multinomial logit via Metropolis-Hastings
+ void MCMCmnlMH(double *sampledata, const int *samplerow,
+ const int *samplecol, const double *Ydata,
+ const int *Yrow, const int *Ycol,
+ const double *Xdata, const int *Xrow, const int *Xcol,
+ const int *burnin, const int *mcmc, const int *thin,
+ const double *tunedata, const int *tunerow,
+ const int *tunecol, const int *lecuyer,
+ const int *seedarray, const int *lecuyerstream,
+ const int *verbose, const double *betastartdata,
+ const int *betastartrow, const int *betastartcol,
+ const double *b0data, const int *b0row,
+ const int *b0col, const double *B0data,
+ const int *B0row, const int *B0col,
+ const double *Vdata, const int *Vrow,
+ const int *Vcol) {
+
+ // pull together Matrix objects
+ // REMEMBER TO ACCESS PASSED ints AND doubles PROPERLY
+ const Matrix <double> Y = r2scythe(*Yrow, *Ycol, Ydata);
+ const Matrix <double> X = r2scythe(*Xrow, *Xcol, Xdata);
+ const Matrix <double> tune = r2scythe(*tunerow, *tunecol, tunedata);
+ Matrix <double> beta = r2scythe(*betastartrow, *betastartcol,
+ betastartdata);
+ const Matrix <double> b0 = r2scythe(*b0row, *b0col, b0data);
+ const Matrix <double> B0 = r2scythe(*B0row, *B0col, B0data);
+ const Matrix <double> V = r2scythe(*Vrow, *Vcol, Vdata);
+
+ // 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 = X.cols();
+
+ // storage matrix or matrices
+ Matrix<double> storemat(nstore, k);
+
+ // initialize rng stream
+ rng *stream = MCMCpack_get_rng(*lecuyer, seedarray, *lecuyerstream);
+
+ // proposal parameters
+ const Matrix<double> propV = tune * invpd(B0 + invpd(V)) * tune;
+ const Matrix<double> propC = cholesky(propV) ;
+
+ double logpost_cur = mnl_logpost(Y, X, beta, b0, B0);
+
+ int count = 0;
+ int accepts = 0;
+ ///// MCMC SAMPLING OCCURS IN THIS FOR LOOP
+ for(int iter = 0; iter < tot_iter; ++iter){
+
+ // sample beta
+ const Matrix<double> beta_can = gaxpy(propC, stream->rnorm(k,1), beta);
+
+ const double logpost_can = mnl_logpost(Y,X,beta_can, b0, B0);
+ const double ratio = ::exp(logpost_can - logpost_cur);
+
+ if (stream->runif() < ratio){
+ beta = beta_can;
+ logpost_cur = logpost_can;
+ ++accepts;
+ }
+
+ // store values in matrices
+ if (iter >= *burnin && ((iter % *thin)==0)){
+ for (int j = 0; j < k; j++)
+ storemat(count, j) = beta[j];
+ ++count;
+ }
+
+ // print output to stdout
+ if(*verbose == 1 && iter % 500 == 0){
+ Rprintf("\n\nMCMCmnl Metropolis iteration %i of %i \n", (iter+1), tot_iter);
+ Rprintf("beta = \n");
+ for (int j=0; j<k; ++j)
+ Rprintf("%10.5f\n", beta[j]);
+ Rprintf("Metropolis acceptance rate for beta = %3.5f\n\n",
+ static_cast<double>(accepts) /
+ static_cast<double>(iter+1));
+ }
+
+ void R_CheckUserInterrupt(void); // allow user interrupts
+ } // end MCMC loop
+
+ delete stream; // clean up random number stream
+
+ // load draws into sample array
+ const int size = *samplerow * *samplecol;
+ for(int i = 0; i < size; ++i)
+ sampledata[i] = storemat[i];
+
+ } // end MCMCmnlMH
+} // end extern "C"
diff --git a/src/MCMCmnlslice.cc b/src/MCMCmnlslice.cc
new file mode 100644
index 0000000..8b97e3a
--- /dev/null
+++ b/src/MCMCmnlslice.cc
@@ -0,0 +1,335 @@
+// MCMCmnlslice.cc DESCRIPTION HERE
+//
+// The initial version of this file was generated by the
+// auto.Scythe.call() function in the MCMCpack R package
+// written by:
+//
+// 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
+//
+// This file was initially generated on Wed Dec 29 15:27:40 2004
+// REVISION HISTORY
+
+#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;
+
+static double mnl_logpost(const Matrix<double>& Y, const Matrix<double>& X,
+ const Matrix<double>& beta,
+ const Matrix<double>& beta_prior_mean,
+ const Matrix<double>& beta_prior_prec){
+
+ // likelihood
+ double loglike = 0.0;
+ Matrix<double> numer = exp(X * beta);
+ numer = reshape(numer, Y.rows(), Y.cols());
+ double *denom = new double[Y.rows()];
+ for (int i=0; i<Y.rows(); ++i){
+ denom[i] = 0.0;
+ for (int j=0; j<Y.cols(); ++j){
+ if (Y(i,j) != -999){
+ denom[i] += numer(i,j);
+ }
+ }
+ for (int j=0; j<Y.cols(); ++j){
+ if (Y(i,j) == 1.0){
+ loglike += std::log(numer(i,j) / denom[i]);
+ }
+ }
+ }
+
+ delete [] denom;
+
+ // prior
+ double logprior = 0.0;
+ if (beta_prior_prec(0,0) != 0){
+ logprior = lndmvn(beta, beta_prior_mean, invpd(beta_prior_prec));
+ }
+
+ return (loglike + logprior);
+}
+
+
+
+
+// eventually all of the slice sampling functions should be made more
+// general and put in MCMCfcds.{h cc}
+//
+// Radford Neal's (2000) doubling procedure coded for a logdensity
+ static void doubling(double (*logfun)(const Matrix<double>&,
+ const Matrix<double>&,
+ const Matrix<double>&,
+ const Matrix<double>&,
+ const Matrix<double>&),
+ const Matrix<double>& beta, const int& index,
+ const double& z,
+ const double& w, const int& p,
+ const Matrix<double>& Y,
+ const Matrix<double>& X,
+ const Matrix<double>& beta_prior_mean,
+ const Matrix<double>& beta_prior_prec,
+ rng* stream, double& L, double& R){
+
+ const double U = stream->runif();
+ const double x0 = beta[index];
+ Matrix<double> beta_L = beta;
+ Matrix<double> beta_R = beta;
+ L = x0 - w*U;
+ beta_L[index] = L;
+ R = L + w;
+ beta_R[index] = R;
+ int K = p;
+ while (K > 0 &&
+ (z < logfun(Y, X, beta_L, beta_prior_mean, beta_prior_prec) |
+ z < logfun(Y, X, beta_R, beta_prior_mean, beta_prior_prec))){
+ double V = stream->runif();
+ if (V < 0.5){
+ L = L - (R - L);
+ beta_L[index] = L;
+ }
+ else {
+ R = R + (R - L);
+ beta_R[index] = R;
+ }
+ --K;
+ }
+}
+
+// Radford Neal's (2000) Accept procedure coded for a logdensity
+static const bool Accept(double (*logfun)(const Matrix<double>&,
+ const Matrix<double>&,
+ const Matrix<double>&,
+ const Matrix<double>&,
+ const Matrix<double>&),
+ const Matrix<double>& beta, const int& index,
+ const double x0,
+ const double& z, const double& w,
+ const Matrix<double>& Y,
+ const Matrix<double>& X,
+ const Matrix<double>& beta_prior_mean,
+ const Matrix<double>& beta_prior_prec,
+ 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 && beta[index] >= M) || (x0 >= M && beta[index] < M)){
+ D = true;
+ }
+ if (beta[index] < M){
+ Rhat = M;
+ }
+ else {
+ Lhat = M;
+ }
+
+ Matrix<double> beta_L = beta;
+ Matrix<double> beta_R = beta;
+ beta_L[index] = Lhat;
+ beta_R[index] = Rhat;
+
+ if (D && z >= logfun(Y, X, beta_L, beta_prior_mean, beta_prior_prec) &&
+ z >= logfun(Y, X, beta_R, beta_prior_mean, beta_prior_prec)){
+ return(false);
+ }
+ }
+ return(true);
+}
+
+
+// Radford Neal's (2000) shrinkage procedure coded for a log density
+static double shrinkage(double (*logfun)(const Matrix<double>&,
+ const Matrix<double>&,
+ const Matrix<double>&,
+ const Matrix<double>&,
+ const Matrix<double>&),
+ const Matrix<double>& beta, const int& index,
+ const double& z, const double& w,
+ const Matrix<double>& Y,
+ const Matrix<double>& X,
+ const Matrix<double>& beta_prior_mean,
+ const Matrix<double>& beta_prior_prec,
+ rng* stream, const double& L, const double& R){
+
+ double Lbar = L;
+ double Rbar = R;
+ Matrix<double> beta_x1 = beta;
+ const double x0 = beta[index];
+ for (;;){
+ const double U = stream->runif();
+ const double x1 = Lbar + U*(Rbar - Lbar);
+ beta_x1[index] = x1;
+ if (z < logfun(Y, X, beta_x1, beta_prior_mean, beta_prior_prec) &&
+ Accept(logfun, beta_x1, index, x0, z, w,
+ Y, X, beta_prior_mean, beta_prior_prec, L, R)){
+ return(x1);
+ }
+ if (x1 < x0){
+ Lbar = x1;
+ }
+ else {
+ Rbar = x1;
+ }
+ } // end infinite loop
+}
+
+
+
+
+
+
+extern "C" {
+
+ // MCMC sampling for MNL model via slice sampling
+ void MCMCmnlslice(double *sampledata, const int *samplerow,
+ const int *samplecol, const double *Ydata,
+ const int *Yrow, const int *Ycol,
+ const double *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 double *betastartdata,
+ const int *betastartrow, const int *betastartcol,
+ const double *b0data, const int *b0row,
+ const int *b0col, const double *B0data,
+ const int *B0row, const int *B0col,
+ const double *Vdata, const int *Vrow,
+ const int *Vcol) {
+
+ // pull together Matrix objects
+ // REMEMBER TO ACCESS PASSED ints AND doubles PROPERLY
+ const Matrix <double> Y = r2scythe(*Yrow, *Ycol, Ydata);
+ const Matrix <double> X = r2scythe(*Xrow, *Xcol, Xdata);
+ Matrix <double> beta = r2scythe(*betastartrow, *betastartcol,
+ betastartdata);
+ const Matrix <double> b0 = r2scythe(*b0row, *b0col, b0data);
+ const Matrix <double> B0 = r2scythe(*B0row, *B0col, B0data);
+ const Matrix <double> V = r2scythe(*Vrow, *Vcol, Vdata);
+
+ // 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 = X.cols();
+
+ // storage matrix or matrices
+ Matrix<double> storemat(nstore, k);
+
+ // initialize rng stream
+ rng *stream = MCMCpack_get_rng(*lecuyer, seedarray, *lecuyerstream);
+
+ // proposal parameters
+ const Matrix<double> propV = invpd(B0 + invpd(V));
+ const Matrix<double> w_init = ones<double>(k, 1);
+ for (int i=0; i<k; ++i)
+ w_init[i] = sqrt(propV(i,i))*0.05;
+
+ // starting values
+ double L = -1.0;
+ double R = 1.0;
+
+
+ const int warmup_iter = 100;
+ const int warmup_burnin = 10;
+ const int p_init = 15;
+ const Matrix<double> widthmat(warmup_iter - warmup_burnin, k);
+ // warm up sampling to choose the slice sampling parameters
+ for (int iter=0; iter<warmup_iter; ++iter){
+ for (int index=0; index<k; ++index){
+
+ double funval = mnl_logpost(Y, X, beta, b0,
+ B0);
+ double z = funval - stream->rexp(1.0);
+ doubling(&mnl_logpost, beta, index, z, w_init[index], p_init, Y, X,
+ b0, B0, stream, L, R);
+
+ beta[index] = shrinkage(&mnl_logpost, beta, index, z,
+ w_init[index], Y, X, b0,
+ B0, stream, L, R);
+ if (iter >= warmup_burnin){
+ widthmat(iter-warmup_burnin, index) = R - L;
+ }
+ }
+ }
+ const Matrix<double> w = meanc(widthmat);
+ Matrix<int> p = ones<int>(k,1);
+ for (int index=0; index<k; ++index){
+ int p_temp = 2;
+ while ((w[index] * pow(2.0, p_temp) ) < max(widthmat(_,index))){
+ ++p_temp;
+ }
+ p[index] = p_temp + 1;
+ }
+
+
+
+ int count = 0;
+ ///// REAL MCMC SAMPLING OCCURS IN THIS FOR LOOP
+ for(int iter = 0; iter < tot_iter; ++iter){
+ for (int index=0; index<k; ++index){
+
+ double funval = mnl_logpost(Y, X, beta, b0,
+ B0);
+ double z = funval - stream->rexp(1.0);
+ doubling(&mnl_logpost, beta, index, z, w[index], p[index],
+ Y, X, b0, B0, stream, L, R);
+ beta[index] = shrinkage(&mnl_logpost, beta, index, z,
+ w[index], Y, X, b0,
+ B0, stream, L, R);
+ }
+
+
+ // store draws in storage matrix (or matrices)
+ if(iter >= *burnin && (iter % *thin == 0)){
+ for (int j = 0; j < k; j++)
+ storemat(count, j) = beta[j];
+ ++count;
+ }
+
+ // print output to stdout
+ if(*verbose == 1 && iter % 50 == 0){
+ Rprintf("\n\nMCMCmnl slice iteration %i of %i \n", (iter+1),
+ tot_iter);
+ Rprintf("beta = \n");
+ for (int j=0; j<k; ++j)
+ Rprintf("%10.5f\n", beta[j]);
+ }
+
+ void R_CheckUserInterrupt(void); // allow user interrupts
+ } // end MCMC loop
+
+ delete stream; // clean up random number stream
+
+ // load draws into sample array
+ const int size = *samplerow * *samplecol;
+ for(int i = 0; i < size; ++i)
+ sampledata[i] = storemat[i];
+
+ } // end MCMCmnlslice
+} // end extern "C"
diff --git a/src/MCMCoprobit.cc b/src/MCMCoprobit.cc
index 5a02404..31e82a2 100644
--- a/src/MCMCoprobit.cc
+++ b/src/MCMCoprobit.cc
@@ -1,190 +1,194 @@
-// MCMCoprobit.cc is C++ code to estimate a ordinalprobit regression
-// model with a multivariate normal prior
-//
-// 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
-//
-// updated to the new version of Scythe 7/26/2004 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;
-
-extern "C"{
-
-
- void MCMCoprobit(double *sampledata, const int *samplerow,
- const int *samplecol, const int *Y,
- const double *Xdata,
- const int *Xrow, const int *Xcol, const int *burnin,
- const int *mcmc, const int *thin, const double* tune,
- const int *lecuyer, const int *seedarray,
- const int *lecuyerstream, const int *verbose,
- const double *betadata, const int *betarow,
- const int *betacol, const double* gammadata,
- const int* gammarow, const int* gammacol,
- const double *b0data, const int *b0row, const int *b0col,
- const double *B0data, const int *B0row, const int *B0col) {
-
- // pull together Matrix objects
- const Matrix <double> X = r2scythe(*Xrow, *Xcol, Xdata);
- Matrix <double> beta = r2scythe(*betarow, *betacol,
- betadata);
- Matrix <double> gamma = r2scythe(*gammarow, *gammacol,
- gammadata);
- const Matrix <double> b0 = r2scythe(*b0row, *b0col, b0data);
- const Matrix <double> B0 = r2scythe(*B0row, *B0col, B0data);
-
- // define constants and from cross-product matrices
- const int tot_iter = *burnin + *mcmc; // total number of mcmc iterations
- const int nstore = *mcmc / *thin; // number of draws to store
- const int k = X.cols();
- const int N = X.rows();
- const int ncat = gamma.rows() - 1;
- const Matrix<double> XpX = crossprod(X);
-
- // storage matrix or matrices
- Matrix<double> storemat(nstore, k+ncat+1);
-
- // initialize rng stream
- rng *stream = MCMCpack_get_rng(*lecuyer, seedarray, *lecuyerstream);
-
- // initialize Z
- Matrix<double> gamma_p = gamma;
- Matrix<double> Z(N,1);
- Matrix<double> Xbeta = X * beta;
-
-
- // Gibbs loop
- int count = 0;
- int accepts = 0;
- for (int iter = 0; iter < tot_iter; ++iter){
-
- // [gamma | Z, beta]
- for (int i=2; i<(ncat); ++i){
- if (i==(ncat-1)){
- gamma_p[i] = stream->rtbnorm_combo(gamma[i], ::pow(tune[0], 2.0),
- gamma_p[i-1]);
- }
- else {
- gamma_p[i] = stream->rtnorm_combo(gamma[i], ::pow(tune[0], 2.0),
- gamma_p[i-1],
- gamma[i+1]);
- }
- }
- double loglikerat = 0.0;
- double loggendenrat = 0.0;
-
- // loop over observations and construct the acceptance ratio
- for (int i=0; i<N; ++i){
- if (Y[i] == ncat){
- loglikerat = loglikerat
- + log(1.0 -
- pnorm(gamma_p[Y[i]-1] - Xbeta[i]) )
- - log(1.0 -
- pnorm(gamma[Y[i]-1] - Xbeta[i]) );
- }
- else if (Y[i] == 1){
- loglikerat = loglikerat
- + log(pnorm(gamma_p[Y[i]] - Xbeta[i]) )
- - log(pnorm(gamma[Y[i]] - Xbeta[i]) );
- }
- else{
- loglikerat = loglikerat
- + log(pnorm(gamma_p[Y[i]] - Xbeta[i]) -
- pnorm(gamma_p[Y[i]-1] - Xbeta[i]) )
- - log(pnorm(gamma[Y[i]] - Xbeta[i]) -
- pnorm(gamma[Y[i]-1] - Xbeta[i]) );
- }
- }
- for (int j=2; j<(ncat-1); ++j){
- loggendenrat = loggendenrat
- + log(pnorm(gamma[j+1], gamma[j], tune[0]) -
- pnorm(gamma[j-1], gamma[j], tune[0]) )
- - log(pnorm(gamma_p[j+1], gamma_p[j], tune[0]) -
- pnorm(gamma_p[j-1], gamma_p[j], tune[0]) );
- }
- double logacceptrat = loglikerat + loggendenrat;
- if (stream->runif() <= exp(logacceptrat)){
- gamma = gamma_p;
- ++accepts;
- }
-
-
- // [Z| gamma, beta, y]
- //Matrix<double> Z_mean = X * beta;
- for (int i=0; i<N; ++i){
- Z[i] = stream->rtnorm_combo(Xbeta[i], 1.0, gamma[Y[i]-1], gamma[Y[i]]);
- }
-
-
- // [beta|Z, gamma]
- const Matrix<double> XpZ = t(X) * Z;
- beta = NormNormregress_beta_draw(XpX, XpZ, b0, B0, 1.0, stream);
- Xbeta = X * beta;
-
-
- // store values in matrices
- if (iter >= *burnin && ((iter % *thin)==0)){
- for (int j=0; j<k; ++j)
- storemat(count, j) = beta[j];
- for (int j=0; j<(ncat+1); ++j)
- storemat(count, j+k) = gamma[j];
- ++count;
- }
-
- // print output to stdout
- if(*verbose == 1 && iter % 500 == 0){
- Rprintf("\n\nMCMCprobit iteration %i of %i \n", (iter+1), tot_iter);
- Rprintf("beta = \n");
- for (int j=0; j<k; ++j)
- Rprintf("%10.5f\n", beta[j]);
- Rprintf("Metropolis acceptance rate for beta = %3.5f\n\n",
- static_cast<double>(accepts) /
- static_cast<double>(iter+1));
- }
-
-
- void R_CheckUserInterrupt(void); // allow user interrupts
- }
+// MCMCoprobit.cc is C++ code to estimate a ordinalprobit regression
+// model with a multivariate normal prior
+//
+// 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
+//
+// updated to the new version of Scythe 7/26/2004 KQ
+// fixed a bug pointed out by Alexander Raach 1/16/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;
+
+extern "C"{
+
+
+ void MCMCoprobit(double *sampledata, const int *samplerow,
+ const int *samplecol, const int *Y,
+ const double *Xdata,
+ const int *Xrow, const int *Xcol, const int *burnin,
+ const int *mcmc, const int *thin, const double* tune,
+ const int *lecuyer, const int *seedarray,
+ const int *lecuyerstream, const int *verbose,
+ const double *betadata, const int *betarow,
+ const int *betacol, const double* gammadata,
+ const int* gammarow, const int* gammacol,
+ const double *b0data, const int *b0row, const int *b0col,
+ const double *B0data, const int *B0row, const int *B0col) {
+
+ // pull together Matrix objects
+ const Matrix <double> X = r2scythe(*Xrow, *Xcol, Xdata);
+ Matrix <double> beta = r2scythe(*betarow, *betacol,
+ betadata);
+ Matrix <double> gamma = r2scythe(*gammarow, *gammacol,
+ gammadata);
+ const Matrix <double> b0 = r2scythe(*b0row, *b0col, b0data);
+ const Matrix <double> B0 = r2scythe(*B0row, *B0col, B0data);
+
+ // define constants and from cross-product matrices
+ const int tot_iter = *burnin + *mcmc; // total number of mcmc iterations
+ const int nstore = *mcmc / *thin; // number of draws to store
+ const int k = X.cols();
+ const int N = X.rows();
+ const int ncat = gamma.rows() - 1;
+ const Matrix<double> XpX = crossprod(X);
+
+ // storage matrix or matrices
+ Matrix<double> storemat(nstore, k+ncat+1);
+
+ // initialize rng stream
+ rng *stream = MCMCpack_get_rng(*lecuyer, seedarray, *lecuyerstream);
+
+ // initialize Z
+ Matrix<double> gamma_p = gamma;
+ Matrix<double> Z(N,1);
+ Matrix<double> Xbeta = X * beta;
+
+
+ // Gibbs loop
+ int count = 0;
+ int accepts = 0;
+ for (int iter = 0; iter < tot_iter; ++iter){
+
+ // [gamma | Z, beta]
+ for (int i=2; i<(ncat); ++i){
+ if (i==(ncat-1)){
+ gamma_p[i] = stream->rtbnorm_combo(gamma[i], ::pow(tune[0], 2.0),
+ gamma_p[i-1]);
+ }
+ else {
+ gamma_p[i] = stream->rtnorm_combo(gamma[i], ::pow(tune[0], 2.0),
+ gamma_p[i-1],
+ gamma[i+1]);
+ }
+ }
+ double loglikerat = 0.0;
+ double loggendenrat = 0.0;
+
+ // loop over observations and construct the acceptance ratio
+ for (int i=0; i<N; ++i){
+ if (Y[i] == ncat){
+ loglikerat = loglikerat
+ + log(1.0 -
+ pnorm(gamma_p[Y[i]-1] - Xbeta[i]) )
+ - log(1.0 -
+ pnorm(gamma[Y[i]-1] - Xbeta[i]) );
+ }
+ else if (Y[i] == 1){
+ loglikerat = loglikerat
+ + log(pnorm(gamma_p[Y[i]] - Xbeta[i]) )
+ - log(pnorm(gamma[Y[i]] - Xbeta[i]) );
+ }
+ else{
+ loglikerat = loglikerat
+ + log(pnorm(gamma_p[Y[i]] - Xbeta[i]) -
+ pnorm(gamma_p[Y[i]-1] - Xbeta[i]) )
+ - log(pnorm(gamma[Y[i]] - Xbeta[i]) -
+ pnorm(gamma[Y[i]-1] - Xbeta[i]) );
+ }
+ }
+ for (int j=2; j<ncat; ++j){
+ loggendenrat = loggendenrat
+ + log(pnorm(gamma[j+1], gamma[j], tune[0]) -
+ pnorm(gamma_p[j-1], gamma[j], tune[0]) )
+ - log(pnorm(gamma_p[j+1], gamma_p[j], tune[0]) -
+ pnorm(gamma[j-1], gamma_p[j], tune[0]) );
+
+ }
+ double logacceptrat = loglikerat + loggendenrat;
+ if (stream->runif() <= exp(logacceptrat)){
+ gamma = gamma_p;
+ ++accepts;
+ }
+
+
+ // [Z| gamma, beta, y]
+ //Matrix<double> Z_mean = X * beta;
+ for (int i=0; i<N; ++i){
+ Z[i] = stream->rtnorm_combo(Xbeta[i], 1.0, gamma[Y[i]-1], gamma[Y[i]]);
+ }
+
+
+ // [beta|Z, gamma]
+ const Matrix<double> XpZ = t(X) * Z;
+ beta = NormNormregress_beta_draw(XpX, XpZ, b0, B0, 1.0, stream);
+ Xbeta = X * beta;
+
+
+ // store values in matrices
+ if (iter >= *burnin && ((iter % *thin)==0)){
+ for (int j=0; j<k; ++j)
+ storemat(count, j) = beta[j];
+ for (int j=0; j<(ncat+1); ++j)
+ storemat(count, j+k) = gamma[j];
+ ++count;
+ }
+
+ // print output to stdout
+ if(*verbose == 1 && iter % 500 == 0){
+ Rprintf("\n\nMCMCoprobit iteration %i of %i \n", (iter+1), tot_iter);
+ Rprintf("beta = \n");
+ for (int j=0; j<k; ++j)
+ Rprintf("%10.5f\n", beta[j]);
+ Rprintf("Metropolis acceptance rate for gamma = %3.5f\n\n",
+ static_cast<double>(accepts) /
+ static_cast<double>(iter+1));
+ }
+
+
+ void R_CheckUserInterrupt(void); // allow user interrupts
+ }
delete stream; // clean up random number stream
-
- // return output
- const int size = *samplerow * *samplecol;
- for (int i=0; i<size; ++i)
- sampledata[i] = storemat[i];
-
- Rprintf("\n\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n");
- Rprintf("The Metropolis acceptance rate for beta was %3.5f",
- static_cast<double>(accepts) / static_cast<double>(tot_iter));
- Rprintf("\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n");
-
- }
-
-}
+
+ // return output
+ const int size = *samplerow * *samplecol;
+ for (int i=0; i<size; ++i)
+ sampledata[i] = storemat[i];
+
+ Rprintf("\n\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n");
+ Rprintf("The Metropolis acceptance rate for beta was %3.5f",
+ static_cast<double>(accepts) / static_cast<double>(tot_iter));
+ Rprintf("\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n");
+
+ }
+
+}
diff --git a/src/MCMCordfactanal.cc b/src/MCMCordfactanal.cc
index 100927e..f531143 100644
--- a/src/MCMCordfactanal.cc
+++ b/src/MCMCordfactanal.cc
@@ -20,6 +20,10 @@
// revised version of older MCMCordfactanal
// 7/16/2004 KQ
// updated to new version of Scythe ADM 7/24/2004
+// fixed a bug pointed out by Alexander Raach 1/16/2005 KQ
+//
+
+#include <iostream>
#include "matrix.h"
#include "distributions.h"
@@ -60,7 +64,8 @@ ordfactanalpost (double* sampledata, const int* samplerow,
const double* Lampprecdata, const int* Lampprecrow,
const int* Lamppreccol, const int* storelambda,
const int* storescores,
- int* accepts, const int* outswitch
+ int* acceptsdata, const int* acceptsrow,
+ const int* acceptscol, const int* outswitch
) {
// put together matrices
@@ -77,6 +82,7 @@ ordfactanalpost (double* sampledata, const int* samplerow,
const Matrix<double> Lambda_prior_prec = r2scythe(*Lampprecrow,
*Lamppreccol,
Lampprecdata);
+ Matrix<int> accepts = r2scythe(*acceptsrow, *acceptscol, acceptsdata);
// initialize rng stream
@@ -99,9 +105,9 @@ ordfactanalpost (double* sampledata, const int* samplerow,
// starting values for phi, Xstar, and gamma_p
Matrix<double> phi = Matrix<double>(N, D-1);
+ //Matrix<double> phi = stream->rnorm(N, D-1);
phi = cbind(ones<double>(N,1), phi);
Matrix<double> Xstar = Matrix<double>(N, K);
- Matrix<double> gamma_p = gamma(_,0);
// storage matrices (row major order)
Matrix<double> Lambda_store;
@@ -158,6 +164,7 @@ ordfactanalpost (double* sampledata, const int* samplerow,
// sample gamma
for (int j=0; j<K; ++j){ // do the sampling for each manifest var
+ Matrix<double> gamma_p = gamma(_,j);
Matrix<double> X_mean = phi * t(Lambda(j,_));
for (int i=2; i<(ncateg[j]); ++i){
if (i==(ncateg[j]-1)){
@@ -198,12 +205,12 @@ ordfactanalpost (double* sampledata, const int* samplerow,
}
}
}
- for (int k=2; k<(ncateg[j]-1); ++k){
+ for (int k=2; k<ncateg[j]; ++k){
loggendenrat = loggendenrat
- + log(pnorm(gamma(k+1,j), gamma(k,j), tune[j]) -
- pnorm(gamma(k-1,j), gamma(k,j), tune[j]) )
- - log(pnorm(gamma_p[k+1], gamma_p[k], tune[j]) -
- pnorm(gamma_p[k-1], gamma_p[k], tune[j]) );
+ + log(pnorm(gamma(k+1,j), gamma(k,j), tune[j]) -
+ pnorm(gamma_p[k-1], gamma(k,j), tune[j]) )
+ - log(pnorm(gamma_p[k+1], gamma_p[k], tune[j]) -
+ pnorm(gamma(k-1,j), gamma_p[k], tune[j]) );
}
double logacceptrat = loglikerat + loggendenrat;
if (stream->runif() <= exp(logacceptrat)){
@@ -211,7 +218,7 @@ ordfactanalpost (double* sampledata, const int* samplerow,
if (gamma(i,j) == 300) break;
gamma(i,j) = gamma_p[i];
}
- ++accepts[0];
+ ++accepts[j];
}
}
@@ -226,8 +233,12 @@ ordfactanalpost (double* sampledata, const int* samplerow,
}
Rprintf("\n");
}
- Rprintf("\nMetropolis-Hastings acceptance rate = %10.5f\n",
- static_cast<double>(*accepts)/(static_cast<double>((iter+1)*K)));
+ Rprintf("\nMetropolis-Hastings acceptance rates = \n");
+ for (int j = 0; j<K; ++j){
+ Rprintf("%6.2f",
+ static_cast<double>(accepts[j]) /
+ static_cast<double>((iter+1)));
+ }
}
if (verbose[0] == 1 && iter % 500 == 0 && *outswitch == 2){
Rprintf("\n\nMCMCirtKd iteration %i of %i \n", (iter+1), tot_iter);
@@ -278,6 +289,9 @@ ordfactanalpost (double* sampledata, const int* samplerow,
int size = *samplerow * *samplecol;
for (int i=0; i<size; ++i)
sampledata[i] = output[i];
+
+ for (int j=0; j<K; ++j)
+ acceptsdata[j] = accepts[j];
}
diff --git a/src/MCMCregress.cc b/src/MCMCregress.cc
index f0ca317..5d732f6 100644
--- a/src/MCMCregress.cc
+++ b/src/MCMCregress.cc
@@ -48,16 +48,17 @@ extern "C" {
// simulate from posterior density and return an mcmc by parameters
// matrix of the posterior density sample
void MCMCregress(double *sampledata, const int *samplerow,
- const int *samplecol, const double *Ydata, const int *Yrow,
- const int *Ycol, const double *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 double *betastartdata, const int *betastartrow,
- 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 *samplecol, const double *Ydata, const int *Yrow,
+ const int *Ycol, const double *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 double *betastartdata, const int *betastartrow,
+ 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) {
+
// pull together Matrix objects
Matrix <double> Y = r2scythe(*Yrow, *Ycol, Ydata);
Matrix <double> X = r2scythe(*Xrow, *Xcol, Xdata);
@@ -122,4 +123,3 @@ extern "C" {
} // end MCMCregress
} // end extern "C"
-
diff --git a/src/MCMCrng.cc b/src/MCMCrng.cc
index f3118da..e800a3b 100644
--- a/src/MCMCrng.cc
+++ b/src/MCMCrng.cc
@@ -26,10 +26,10 @@
#include "mersenne.h"
#include "lecuyer.h"
-#include <R.h> // needed to use Rprintf()
+//#include <R.h> // needed to use Rprintf()
using namespace std;
-using namespace SCYTHE;
+//using namespace SCYTHE;
namespace SCYTHE {
diff --git a/src/MCMCrng.h b/src/MCMCrng.h
index b84fc1b..8acc11a 100644
--- a/src/MCMCrng.h
+++ b/src/MCMCrng.h
@@ -26,10 +26,10 @@
#include "mersenne.h"
#include "lecuyer.h"
-#include <R.h> // needed to use Rprintf()
+//#include <R.h> // needed to use Rprintf()
using namespace std;
-using namespace SCYTHE;
+//using namespace SCYTHE;
namespace SCYTHE {
diff --git a/src/MCMCregress.cc b/src/MCMCtobit.cc
similarity index 50%
copy from src/MCMCregress.cc
copy to src/MCMCtobit.cc
index f0ca317..4fca55b 100644
--- a/src/MCMCregress.cc
+++ b/src/MCMCtobit.cc
@@ -1,5 +1,7 @@
-// MCMCregress.cc is a program that simualates draws from the posterior
-// density of a linear regression model with Gaussian errors.
+// MCMCtobit.cc is a program that simualates draws from the posterior
+// density of a linear regression model with Gaussian errors when the
+// dependent variable is censored from below and/or above.
+//
//
// The initial version of this file was generated by the
// auto.Scythe.call() function in the MCMCpack R package
@@ -21,12 +23,9 @@
//
// Copyright (C) 2004 Andrew D. Martin and Kevin M. Quinn
//
-// This file was initially generated on Fri Jul 23 15:07:21 2004
-//
+// This file was initially generated on Tue Sep 14 00:50:08 2004
// ADM and KQ 10/10/2002 [ported to Scythe0.3]
-// ADM 6/2/04 [re-written using template]
-// KQ 6/18/04 [modified to meet new developer specification]
-// ADM 7/22/04 [modified to work with new Scythe and rngs]
+// BG 09/18/2004 [updated to new specification, added above censoring]
#include "matrix.h"
#include "distributions.h"
@@ -37,6 +36,7 @@
#include "MCMCrng.h"
#include "MCMCfcds.h"
+
#include <R.h> // needed to use Rprintf()
#include <R_ext/Utils.h> // needed to allow user interrupts
@@ -45,54 +45,73 @@ using namespace std;
extern "C" {
- // simulate from posterior density and return an mcmc by parameters
- // matrix of the posterior density sample
- void MCMCregress(double *sampledata, const int *samplerow,
- const int *samplecol, const double *Ydata, const int *Yrow,
- const int *Ycol, const double *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 double *betastartdata, const int *betastartrow,
- 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) {
+ // MCMCtobit is a linear regression model with a censored dependent variable
+ void MCMCtobit(double *sampledata, const int *samplerow,
+ const int *samplecol,
+ const double *Ydata, const int *Yrow, const int *Ycol,
+ const double *Xdata,
+ const int *Xrow, const int *Xcol,
+ const double *below, const double *above,
+ const int *burnin, const int *mcmc, const int *thin,
+ const int *lecuyer,
+ const int *seedarray, const int *lecuyerstream,
+ const int *verbose,
+ const double *betastartdata, const int *betastartrow,
+ 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) {
// pull together Matrix objects
+ // REMEMBER TO ACCESS PASSED ints AND doubles PROPERLY
Matrix <double> Y = r2scythe(*Yrow, *Ycol, Ydata);
Matrix <double> X = r2scythe(*Xrow, *Xcol, Xdata);
- Matrix <double> betastart = r2scythe(*betastartrow,
- *betastartcol, betastartdata);
+ Matrix <double> betastart = r2scythe(*betastartrow, *betastartcol,
+ betastartdata);
Matrix <double> b0 = r2scythe(*b0row, *b0col, b0data);
Matrix <double> B0 = r2scythe(*B0row, *B0col, B0data);
-
- // define constants and form cross-product matrices
- const int tot_iter = *burnin + *mcmc; // total number of mcmc iterations
- const int nstore = *mcmc / *thin; // number of draws to store
+ Matrix <double> t_X = t(X);
+
+
+ // 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 = X.cols ();
+ const int N = X.rows ();
const Matrix <double> XpX = crossprod(X);
- const Matrix <double> XpY = t(X) * Y;
-
- // storage matrices
+
+ // storage matrix or matrices
Matrix <double> betamatrix (k, nstore);
Matrix <double> sigmamatrix (1, nstore);
-
+
// initialize rng stream
rng *stream = MCMCpack_get_rng(*lecuyer, seedarray, *lecuyerstream);
// set starting values
Matrix <double> beta = betastart;
-
- // Gibbs sampler
+
+ ///// MCMC SAMPLING OCCURS IN THIS FOR LOOP
int count = 0;
- 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);
-
+ Matrix <double> Z = Y;
+ for(int iter = 0; iter < tot_iter; ++iter){
+ double sigma2 = NormIGregress_sigma2_draw (X, Z, beta, *c0, *d0,
+ stream);
+ Matrix <double> Z_mean = X * beta;
+
+ for (int i=0; i<N; ++i) {
+ if (Y[i] <= *below)
+ Z[i] = stream->rtanorm_combo(Z_mean[i], sigma2, *below);
+ if (Y[i] >= *above)
+ Z[i] = stream->rtbnorm_combo(Z_mean[i], sigma2, *above);
+ }
+
+ Matrix <double> XpZ = t (X) * Z;
+
+ beta = NormNormregress_beta_draw (XpX, XpZ, b0, B0, sigma2, stream);
+
// store draws in storage matrix (or matrices)
- if (iter >= *burnin && (iter % *thin == 0)) {
+ if (iter >= *burnin && (iter % 1 == 0)) {
sigmamatrix (0, count) = sigma2;
for (int j = 0; j < k; j++)
betamatrix (j, count) = beta[j];
@@ -100,26 +119,25 @@ extern "C" {
}
// print output to stdout
- if(*verbose == 1 && iter % 500 == 0) {
- Rprintf("\n\nMCMCregress iteration %i of %i \n",
- (iter+1), tot_iter);
+ if(*verbose == 1 && iter % 1000 == 0) {
+ Rprintf("\n\nMCMCtobit iteration %i of %i \n",
+ (iter+1), tot_iter);
Rprintf("beta = \n");
for (int j=0; j<k; ++j)
Rprintf("%10.5f\n", beta[j]);
Rprintf("sigma2 = %10.5f\n", sigma2);
}
-
+
void R_CheckUserInterrupt(void); // allow user interrupts
} // end MCMC loop
-
+
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];
-
- } // end MCMCregress
+
+ } // end MCMCtobit
} // end extern "C"
-
diff --git a/src/Makevars b/src/Makevars
index 98c8c6d..ae8cde5 100644
--- a/src/Makevars
+++ b/src/Makevars
@@ -1,2 +1,2 @@
-PKG_CXXFLAGS = -Wall -pedantic -DSCYTHE_COMPILE_DIRECT -DSCYTHE_NO_RANGE -O3 -fno-gcse -funroll-loops -DHAVE_TRUNC
+PKG_CXXFLAGS = -DSCYTHE_COMPILE_DIRECT -DSCYTHE_NO_RANGE
diff --git a/src/Makevars.in b/src/Makevars.in
index 37e0a7e..f71ce14 100644
--- a/src/Makevars.in
+++ b/src/Makevars.in
@@ -1,2 +1,2 @@
-PKG_CXXFLAGS = -Wall -pedantic -DSCYTHE_COMPILE_DIRECT -DSCYTHE_NO_RANGE -O3 -fno-gcse -funroll-loops @MV_HAVE_IEEEFP_H@ @MV_HAVE_TRUNC@
+PKG_CXXFLAGS = -DSCYTHE_COMPILE_DIRECT -DSCYTHE_NO_RANGE @MV_HAVE_IEEEFP_H@ @MV_HAVE_TRUNC@
diff --git a/src/optimize.cc b/src/optimize.cc
index 4472e65..fd6528c 100644
--- a/src/optimize.cc
+++ b/src/optimize.cc
@@ -1,642 +1,642 @@
-/*
- * Scythe Statistical Library
- * Copyright (C) 2000-2002 Andrew D. Martin and Kevin M. Quinn;
- * 2002-2004 Andrew D. Martin, Kevin M. Quinn, and Daniel
- * Pemstein. All Rights Reserved.
- *
- * This program is free software; you can redistribute it and/or modify
- * under the terms of the GNU General Public License as published by
- * Free Software Foundation; either version 2 of the License, or (at
- * your option) any later version. See the text files COPYING
- * and LICENSE, distributed with this source code, for further
- * information.
- * --------------------------------------------------------------------
- * scythestat/optimize.cc
- *
- * Provides implementations of various numerical optimization
- * routines.
- *
- */
-
-#ifndef SCYTHE_OPTIMIZE_CC
-#define SCYTHE_OPTIMIZE_CC
-
-#include <cmath>
-#include <iostream>
-
-#ifdef SCYTHE_COMPILE_DIRECT
-#include "optimize.h"
-#include "error.h"
-#include "util.h"
-#include "distributions.h"
-#include "la.h"
-#include "ide.h"
-#include "smath.h"
-#include "stat.h"
-#else
-#include "scythestat/optimize.h"
-#include "scythestat/error.h"
-#include "scythestat/util.h"
-#include "scythestat/distributions.h"
-#include "scythestat/la.h"
-#include "scythestat/ide.h"
-#include "scythestat/smath.h"
-#include "scythestat/stat.h"
-#endif
-
-// Avoid NameSpace Pollution
-namespace SCYTHE {
-
- /* Functions (private to this file) that do very little... */
-#ifdef __MINGW32__
- template <class T>
- static
- T
- donothing(const Matrix<T> &x) {return 0.0;}
-
- template <class T>
- static
- T
- donothing(const T &x) { return 0.0; }
-#else
- namespace {
- template <class T>
- T
- donothing(const Matrix<T> &x) {return 0.0;}
-
- template <class T>
- T
- donothing(const T &x) { return 0.0; }
- }
-#endif
-
-
- /* Return the machine epsilon
- * Notes: Algorithm taken from Sedgewick, Robert. 1992. Algorithms
- * in C++. Addison Wesley. pg. 561
- */
- template <class T>
- T
- epsilon()
- {
- T eps, del, neweps;
- del = (T) 0.5;
- eps = (T) 0.0;
- neweps = (T) 1.0;
-
- while ( del > 0 ) {
- if ( 1 + neweps > 1 ) { /* Then the value might be too large */
- eps = neweps; /* ...save the current value... */
- neweps -= del; /* ...and decrement a bit */
- } else { /* Then the value is too small */
- neweps += del; /* ...so increment it */
- }
- del *= 0.5; /* Reduce the adjustment by half */
- }
-
- return eps;
- }
-
- /* Calculate the definite integral of a function from a to b */
- template <class T>
- T
- intsimp(T (*fun)(const T &), const T &a, const T &b, const int &N)
- {
- if (a > b)
- throw scythe_invalid_arg (__FILE__, __PRETTY_FUNCTION__, __LINE__,
- "Lower limit larger than upper");
-
- if (N <= 0)
- throw scythe_invalid_arg (__FILE__, __PRETTY_FUNCTION__, __LINE__,
- "Number of subintervals negative");
-
- T I = (T) 0;
- T w = (b - a) / N;
- for (int i = 1; i <= N; i++)
- I += w * (fun(a +(i - 1) *w) + 4 * fun(a - w / 2 + i * w) +
- fun(a + i * w)) / 6;
-
- return I;
- }
-
- /* Calculate the definite integral of a function from a to b
- * Notes: Algorithm taken from Sedgewick, Robert. 1992. Algorithms
- * in C++. Addison Wesley. pg. 562
- */
- template <class T>
- T
- adaptsimp(T (*fun)(const T &), const T &a, const T &b, const int &N,
- const T &tol)
- {
- if (a > b)
- throw scythe_invalid_arg (__FILE__, __PRETTY_FUNCTION__, __LINE__,
- "Lower limit larger than upper");
-
- if (N <= 0)
- throw scythe_invalid_arg (__FILE__, __PRETTY_FUNCTION__, __LINE__,
- "Number of subintervals negative");
-
- T I = intsimp(fun, a, b, N);
- if (::fabs(I - intsimp(fun, a, b, N / 2)) > tol)
- return adaptsimp(fun, a, (a + b) / 2, N, tol)
- + adaptsimp(fun, (a + b) / 2, b, N, tol);
-
- return I;
- }
-
-
- /* Numerically calculates the first derivative of a function
- * Notes: Algorithm taken from Nocedal and Wright. 1999. section 7.1
- * with some additional tricks from Press et al. NRC
- */
- template <class T>
- Matrix<T>
- gradfdif (T (*fun)(const Matrix<T> &, const Matrix<T> &,
- const Matrix<T> &), const Matrix<T> &theta,
- const Matrix<T> &y, const Matrix<T> &X)
- {
- if (! theta.isColVector())
- throw scythe_dimension_error (__FILE__, __PRETTY_FUNCTION__,
- __LINE__, "Theta not column vector");
-
- int k = theta.size();
- // stepsize CAREFUL-- THIS IS MACHINE-SPECIFIC!!!
- T h = std::sqrt(epsilon<T>()); // 2.25e-16
- //T h = std::sqrt(2.25e-16);
-
-
- Matrix<T> grad(k,1);
-
- for (int i = 0; i < k; ++i){
- Matrix<T> e(k,1);
- e[i] = h;
- Matrix<T> temp = theta + e;
- donothing(temp);
- e = temp - theta;
- grad[i] = (fun(theta + e, y, X) - fun(theta, y, X)) / e[i];
- }
-
- return grad;
- }
-
-
- /* Numerically calculates the first derivative of a function
- * Notes: Algorithm taken from Nocedal and Wright. 1999. section 7.1
- * with some additional tricks from Press et al. NRC
- */
- template <class T>
- T
- gradfdifls (T (*fun)(const Matrix<T> &, const Matrix<T> &,
- const Matrix<T> &), const T &alpha,
- const Matrix<T> &theta, const Matrix<T> &p,
- const Matrix<T> &y, const Matrix<T> &X)
- {
- if (! theta.isColVector())
- throw scythe_dimension_error (__FILE__, __PRETTY_FUNCTION__,
- __LINE__, "Theta not column vector");
- if (! p.isColVector())
- throw scythe_dimension_error (__FILE__, __PRETTY_FUNCTION__,
- __LINE__, "p not column vector");
-
- int k = theta.size();
- // stepsize CAREFUL-- THIS IS MACHINE-SPECIFIC!!!
- T h = std::sqrt(epsilon<T>()); //2.2e-16
- //T h = std::sqrt(2.2e-16);
-
- T deriv;
-
- for (int i = 0; i < k; ++i) {
- T temp = alpha + h;
- donothing(temp);
- T e = temp - alpha;
- deriv = (fun(theta + (alpha + e) * p, y, X)
- - fun(theta + alpha * p, y, X)) / e;
- }
-
- return deriv;
- }
-
-
-
- /* Numerically calculates the gradient of a function
- * Notes: Algorithm taken from Nocedal and Wright. 1999. section 7.1
- * with some additional tricks from Press et al. NRC
- */
- template <class T>
- Matrix<T>
- jacfdif(Matrix<T> (*fun)(const Matrix<T> &), const Matrix<T> &theta)
- {
- if (! theta.isColVector())
- throw scythe_dimension_error (__FILE__, __PRETTY_FUNCTION__,
- __LINE__, "Theta not column vector");
-
- Matrix<T> fval = fun(theta);
-
- int k = theta.rows();
- int n = fval.rows();
- // stepsize CAREFUL -- THIS IS MACHINE-SPECIFIC!!!
- T h = std::sqrt(epsilon<T>()); //2.2e-16
- //T h = std::sqrt(2.2e-16);
- Matrix<T> J(n,k);
-
- for (int i = 0; i < k; ++i) {
- Matrix<T> e(k,1);
- e[i] = h;
- Matrix<T> temp = theta + e;
- donothing(temp);
- e = temp - theta;
- Matrix<T> fthetae = fun(theta + e);
- Matrix<T> ftheta = fun(theta);
- for (int j = 0; j < n; ++j) {
- J(j,i) = (fthetae[j] - ftheta[j]) / e[i];
- }
- }
-
- return J;
- }
-
- /* Numerically calculates the gradient of a function
- * Notes: Algorithm taken from Nocedal and Wright. 1999. section 7.1
- * with some additional tricks from Press et al. NRC
- */
- template <class T>
- Matrix<T>
- jacfdif(Matrix<T> (*fun)(const Matrix<T> &, const Matrix<T> &),
- const Matrix<T> &theta, const Matrix<T> &psi)
- {
- if (! theta.isColVector())
- throw scythe_dimension_error (__FILE__, __PRETTY_FUNCTION__,
- __LINE__, "Theta not column vector");
-
- Matrix<T> fval = fun(theta, psi);
-
- int k = theta.rows();
- int n = fval.rows();
- // stepsize CAREFUL -- THIS IS MACHINE-SPECIFIC!!!
- //T h = std::sqrt(epsilon<T>()); //2.2e-16
- T h = std::sqrt(2.2e-16);
- Matrix<T> J(n, k);
-
- for (int i = 0; i < k; ++i) {
- Matrix<T> e(k,1);
- e[i] = h;
- Matrix<T> temp = theta + e;
- donothing(temp);
- e = temp - theta;
- Matrix<T> fthetae = fun(theta + e, psi);
- Matrix<T> ftheta = fun(theta, psi);
- for (int j = 0; j < n; ++j) {
- J(j,i) = (fthetae[j] - ftheta[j]) / e[i];
- }
- }
-
- return J;
- }
-
- /* Numerically calculates the Hessian of a function */
- template <class T>
- Matrix<T>
- hesscdif (T (*fun)(const Matrix<T> &, const Matrix<T> &,
- const Matrix<T> &), const Matrix<T> &theta,
- const Matrix<T> &y, const Matrix<T> &X)
- {
- if (! theta.isColVector())
- throw scythe_dimension_error (__FILE__, __PRETTY_FUNCTION__,
- __LINE__, "Theta not column vector");
-
- T fval = fun(theta,y,X);
-
- int k = theta.rows();
-
- // stepsize CAREFUL -- THIS IS MACHINE-SPECIFIC!!!
- T h2 = (T) 1e-10;
- T h = ::sqrt(h2);
- Matrix<T> H(k,k);
-
- for (int i=0; i<k; ++i) {
- Matrix<T> ei(k, 1);
- ei[i] = h;
- Matrix<T> temp = theta + ei;
- donothing(temp);
- ei = temp - theta;
- for (int j = 0; j < k; ++j){
- Matrix<T> ej(k,1);
- ej[j] = h;
- temp = theta + ej;
- donothing(temp);
- ej = temp - theta;
-
- if (i==j){
- H(i,i) = ( -fun(theta + 2.0 * ei, y, X) + 16.0 *
- fun(theta+ei,y,X) - 30.0 * fval + 16.0 *
- fun(theta-ei,y,X) -
- fun(theta-2.0 * ei, y, X)) / (12.0 * h2);
- } else {
- H(i,j) = ( fun(theta + ei + ej, y, X) - fun(theta+ei-ej, y, X)
- - fun(theta - ei + ej, y, X) + fun(theta-ei-ej, y, X))
- / (4.0 * h2);
- }
- }
- }
-
- return H;
- }
-
- /* Performs a linesearch to find the step length (\a alpha)
- * Notes: Algorithm taken from Nocedal and Wright. 1999. Procedure
- * 3.1
- */
- template <class T>
- T
- linesearch1(T (*fun)(const Matrix<T> &, const Matrix<T> &,
- const Matrix<T> &), const Matrix<T> &theta,
- const Matrix<T> &p, const Matrix<T> &y,
- const Matrix<T> &X)
- {
- if (! theta.isColVector())
- throw scythe_dimension_error (__FILE__, __PRETTY_FUNCTION__,
- __LINE__, "Theta not column vector");
- if (! p.isColVector())
- throw scythe_dimension_error (__FILE__, __PRETTY_FUNCTION__,
- __LINE__, "p not column vector");
-
- T alpha_bar = (T) 1.0;
- T rho = (T) 0.9;
- T c = (T) 0.5;
- T alpha = alpha_bar;
- Matrix<T> fgrad = gradfdif(fun, theta, y, X);
-
- while (fun(theta + alpha * p, y, X) > (fun(theta, y, X) + c
- * alpha * t(fgrad) * p)[0]) {
- alpha = rho * alpha;
- }
-
- return alpha;
- }
-
- /* Performs a linesearch to find the step length (\a alpha)
- * Notes: Algorithm taken from Nocedal and Wright. 1999. Algorithm
- * 3.2
- */
- template <class T>
- T
- linesearch2(T (*fun)(const Matrix<T> &, const Matrix<T> &,
- const Matrix<T> &), const Matrix<T> &theta,
- const Matrix<T> &p, const Matrix<T> &y,
- const Matrix<T> &X, rng *myrng)
- {
- if (! theta.isColVector())
- throw scythe_dimension_error (__FILE__, __PRETTY_FUNCTION__,
- __LINE__, "Theta not column vector");
- if (! p.isColVector())
- throw scythe_dimension_error (__FILE__, __PRETTY_FUNCTION__,
- __LINE__, "p not column vector");
-
- T alpha_last = (T) 0.0;
- T alpha_cur = (T) 1.0;
- T alpha_max = (T) 10.0;
- T c1 = (T) 1e-4;
- T c2 = (T) 0.5;
- int max_iter = 50;
- T fgradalpha0 = gradfdifls(fun, (T) 0, theta, p, y, X);
-
- for (int i = 0; i < max_iter; ++i) {
- T phi_cur = fun(theta + alpha_cur * p, y, X);
- T phi_last = fun(theta + alpha_last * p, y, X);
-
- if ((phi_cur > (fun(theta, y, X) + c1 * alpha_cur * fgradalpha0))
- ||
- ((phi_cur >= phi_last) && (i > 0))) {
- T alphastar = zoom(fun, alpha_last, alpha_cur, theta, p, y, X);
- return alphastar;
- }
-
- T fgradalpha_cur = gradfdifls(fun, alpha_cur, theta, p, y, X);
- if ( ::fabs(fgradalpha_cur) <= -1 * c2 * fgradalpha0)
- return alpha_cur;
-
- if ( fgradalpha_cur >= (T) 0.0) {
- T alphastar = zoom(fun, alpha_cur, alpha_last, theta, p, y, X);
- return alphastar;
- }
-
- alpha_last = alpha_cur;
- alpha_cur = myrng->runif() * (alpha_max - alpha_cur) + alpha_cur;
- }
-
- return 0.001;
- }
-
-
- /* Finds the minimum of a function once bracketed (i.e. over a
- * closed interval).
- * Notes: Algorithm taken from Nocedal and Wright. 1999. Algorithm
- * 3.3
- */
- template <class T>
- T
- zoom (T (*fun)(const Matrix<T> &, const Matrix<T> &,
- const Matrix<T> &), const T &alo, const T &ahi,
- const Matrix<T> &theta, const Matrix<T> &p,
- const Matrix<T> &y, const Matrix<T> &X )
- {
- if (! theta.isColVector())
- throw scythe_dimension_error (__FILE__, __PRETTY_FUNCTION__,
- __LINE__, "Theta not column vector");
- if (! p.isColVector())
- throw scythe_dimension_error (__FILE__, __PRETTY_FUNCTION__,
- __LINE__, "p not column vector");
-
- T alpha_lo = alo;
- T alpha_hi = ahi;
- T alpha_j = (alo + ahi) / 2.0;
- T phi_0 = fun(theta, y, X);
- T c1 = (T) 1e-4;
- T c2 = (T) 0.5;
- T fgrad0 = gradfdifls(fun,(T) 0, theta, p, y, X);
-
- int count = 0;
- int maxit = 20;
- while(count < maxit) {
- T phi_j = fun(theta + alpha_j * p, y, X);
- T phi_lo = fun(theta + alpha_lo * p, y, X);
-
- if ((phi_j > (phi_0 + c1 * alpha_j * fgrad0))
- || (phi_j >= phi_lo)){
- alpha_hi = alpha_j;
- } else {
- T fgradj = gradfdifls(fun, alpha_j, theta, p, y, X);
- if (::fabs(fgradj) <= -1 * c2 * fgrad0){
- return alpha_j;
- }
- if ( fgradj * (alpha_hi - alpha_lo) >= 0){
- alpha_hi = alpha_lo;
- }
- alpha_lo = alpha_j;
- }
- ++count;
- }
-
- return alpha_j;
- }
-
-
-
- /* Find the minimum of a function using the BFGS algorithm
- * Notes: Algorithm taken from Nocedal and Wright. 1999. algorithm
- * 8.1
- */
- template <class T>
- Matrix<T>
- BFGS (T (*fun)(const Matrix<T> &, const Matrix<T> &,
- const Matrix<T> &), rng *myrng,
- const Matrix<T> &theta,
- const Matrix<T> &y, const Matrix<T> &X,
- const int &maxit, const T &tolerance)
- {
- if (! theta.isColVector())
- throw scythe_dimension_error (__FILE__, __PRETTY_FUNCTION__,
- __LINE__, "Theta not column vector");
-
- int n = theta.size();
-
- // H is initial inverse hessian
- Matrix<T> H = inv(hesscdif(fun, theta, y, X));
-
- // gradient at starting values
- Matrix<T> fgrad = gradfdif(fun, theta, y, X);
- Matrix<T> thetamin = theta;
- Matrix<T> fgrad_new = fgrad;
- Matrix<T> I = eye<T>(n);
-
- int count = 0;
- while( (t(fgrad_new)*fgrad_new)[0] > tolerance) {
- Matrix<T> p = -1 * H * fgrad;
- T alpha = linesearch2(fun, thetamin, p, y, X, myrng);
- Matrix<T> thetamin_new = thetamin + alpha*p;
- fgrad_new = gradfdif(fun, thetamin_new, y, X);
- Matrix<T> s = thetamin_new - thetamin;
- Matrix<T> y = fgrad_new - fgrad;
- T rho = 1.0 / (t(y) * s)[0];
- H = (I - rho * s * t(y)) * H *(I - rho * y * t(s))
- + rho * s * (!s);
-
- thetamin = thetamin_new;
- fgrad = fgrad_new;
- ++count;
-
- std::cout << "BFGS iteration = " << count << std::endl;
- std::cout << "thetamin = " << (!thetamin).toString() << std::endl;
- std::cout << "gradient = " << (!fgrad).toString() << std::endl;
- std::cout << "t(gradient) * gradient = " <<
- ((!fgrad) * fgrad).toString() << std::endl;
-
- if (count > maxit)
- throw scythe_convergence_error (__FILE__, __PRETTY_FUNCTION__,
- __LINE__, "Failed to converge. Try better starting values");
- }
-
- return thetamin;
- }
-
-
- /* Zero a function using Broyen's Method
- * Notes: Algorithm taken from Nocedal and Wright. 1999. algorithm 11
- * line search is not used to determine alpha (this should probably
- * be changed at some point.
- */
- template <class T>
- Matrix<T>
- nls_broyden(Matrix<T> (*fun)(const Matrix<T> &),
- const Matrix<T> &theta, const int &maxit = 5000,
- const T &tolerance = 1e-6)
- {
- if (! theta.isColVector())
- throw scythe_dimension_error (__FILE__, __PRETTY_FUNCTION__,
- __LINE__, "Theta not column vector");
-
- Matrix<T> thetastar = theta;
- Matrix<T> B = jacfdif(fun, thetastar);
-
- Matrix<T> fthetastar;
- Matrix<T> p;
- Matrix<T> thetastar_new;
- Matrix<T> fthetastar_new;
- Matrix<T> s;
- Matrix<T> y;
-
- for (int i = 0; i < maxit; ++i) {
- fthetastar = fun(thetastar);
- p = lu_solve(B, -1 * fthetastar);
- T alpha = (T) 1.0;
- thetastar_new = thetastar + alpha*p;
- fthetastar_new = fun(thetastar_new);
- s = thetastar_new - thetastar;
- y = fthetastar_new - fthetastar;
- B = B + ((y - B * s) * (!s)) / ((!s) * s);
- thetastar = thetastar_new;
- if (max(fabs(fthetastar_new)) < tolerance)
- return thetastar;
- }
-
- throw scythe_convergence_error (__FILE__, __PRETTY_FUNCTION__,
- __LINE__,std::string("Failed to converge. Try better starting") &
- " values or increase maxit");
- }
-
-
- /* Zero a function using Broyen's Method
- * Notes: Algorithm taken from Nocedal and Wright. 1999. algorithm
- * 11.3
- * line search is not used to determine alpha (this should probably
- * be changed at some point.
- */
- template <class T>
- Matrix<T>
- nls_broyden(Matrix<T> (*fun)(const Matrix<T> &, const Matrix<T> &),
- const Matrix<T> &theta, const Matrix<T> &psi,
- const int &maxit=5000, const T &tolerance=1e-6)
- {
- if (! theta.isColVector())
- throw scythe_dimension_error (__FILE__, __PRETTY_FUNCTION__,
- __LINE__, "Theta not column vector");
- if (! psi.isColVector())
- throw scythe_dimension_error (__FILE__, __PRETTY_FUNCTION__,
- __LINE__, "Psi not column vector");
-
- Matrix<T> thetastar = theta;
- Matrix<T> B = jacfdif(fun, thetastar, psi);
-
- Matrix<T> fthetastar;
- Matrix<T> p;
- Matrix<T> thetastar_new;
- Matrix<T> fthetastar_new;
- Matrix<T> s;
- Matrix<T> y;
-
- for (int i = 0; i < maxit; ++i) {
- fthetastar = fun(thetastar, psi);
- p = lu_solve(B, -1 * fthetastar);
- T alpha = (T) 1.0;
- thetastar_new = thetastar + alpha*p;
- fthetastar_new = fun(thetastar_new, psi);
- s = thetastar_new - thetastar;
- y = fthetastar_new - fthetastar;
- B = B + ((y - B * s) * (!s)) / ((!s) * s);
- thetastar = thetastar_new;
- if (max(fabs(fthetastar_new)) < tolerance)
- return thetastar;
- }
-
- throw scythe_convergence_error (__FILE__, __PRETTY_FUNCTION__,
- __LINE__,std::string("Failed to converge. Try better starting") &
- " values or increase maxit");
- }
-
-} // namespace dec
-
-#ifndef SCYTHE_COMPILE_DIRECT
-#include "scythestat/eti/optimize.t"
-#endif
-
-#endif /* SCYTHE_OPTIMIZE_CC */
+/*
+ * Scythe Statistical Library
+ * Copyright (C) 2000-2002 Andrew D. Martin and Kevin M. Quinn;
+ * 2002-2004 Andrew D. Martin, Kevin M. Quinn, and Daniel
+ * Pemstein. All Rights Reserved.
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * under the terms of the GNU General Public License as published by
+ * Free Software Foundation; either version 2 of the License, or (at
+ * your option) any later version. See the text files COPYING
+ * and LICENSE, distributed with this source code, for further
+ * information.
+ * --------------------------------------------------------------------
+ * scythestat/optimize.cc
+ *
+ * Provides implementations of various numerical optimization
+ * routines.
+ *
+ */
+
+#ifndef SCYTHE_OPTIMIZE_CC
+#define SCYTHE_OPTIMIZE_CC
+
+#include <cmath>
+#include <iostream>
+
+#ifdef SCYTHE_COMPILE_DIRECT
+#include "optimize.h"
+#include "error.h"
+#include "util.h"
+#include "distributions.h"
+#include "la.h"
+#include "ide.h"
+#include "smath.h"
+#include "stat.h"
+#else
+#include "scythestat/optimize.h"
+#include "scythestat/error.h"
+#include "scythestat/util.h"
+#include "scythestat/distributions.h"
+#include "scythestat/la.h"
+#include "scythestat/ide.h"
+#include "scythestat/smath.h"
+#include "scythestat/stat.h"
+#endif
+
+// Avoid NameSpace Pollution
+namespace SCYTHE {
+
+ /* Functions (private to this file) that do very little... */
+#ifdef __MINGW32__
+ template <class T>
+ static
+ T
+ donothing(const Matrix<T> &x) {return 0.0;}
+
+ template <class T>
+ static
+ T
+ donothing(const T &x) { return 0.0; }
+#else
+ namespace {
+ template <class T>
+ T
+ donothing(const Matrix<T> &x) {return 0.0;}
+
+ template <class T>
+ T
+ donothing(const T &x) { return 0.0; }
+ }
+#endif
+
+
+ /* Return the machine epsilon
+ * Notes: Algorithm taken from Sedgewick, Robert. 1992. Algorithms
+ * in C++. Addison Wesley. pg. 561
+ */
+ template <class T>
+ T
+ epsilon()
+ {
+ T eps, del, neweps;
+ del = (T) 0.5;
+ eps = (T) 0.0;
+ neweps = (T) 1.0;
+
+ while ( del > 0 ) {
+ if ( 1 + neweps > 1 ) { /* Then the value might be too large */
+ eps = neweps; /* ...save the current value... */
+ neweps -= del; /* ...and decrement a bit */
+ } else { /* Then the value is too small */
+ neweps += del; /* ...so increment it */
+ }
+ del *= 0.5; /* Reduce the adjustment by half */
+ }
+
+ return eps;
+ }
+
+ /* Calculate the definite integral of a function from a to b */
+ template <class T>
+ T
+ intsimp(T (*fun)(const T &), const T &a, const T &b, const int &N)
+ {
+ if (a > b)
+ throw scythe_invalid_arg (__FILE__, __PRETTY_FUNCTION__, __LINE__,
+ "Lower limit larger than upper");
+
+ if (N <= 0)
+ throw scythe_invalid_arg (__FILE__, __PRETTY_FUNCTION__, __LINE__,
+ "Number of subintervals negative");
+
+ T I = (T) 0;
+ T w = (b - a) / N;
+ for (int i = 1; i <= N; i++)
+ I += w * (fun(a +(i - 1) *w) + 4 * fun(a - w / 2 + i * w) +
+ fun(a + i * w)) / 6;
+
+ return I;
+ }
+
+ /* Calculate the definite integral of a function from a to b
+ * Notes: Algorithm taken from Sedgewick, Robert. 1992. Algorithms
+ * in C++. Addison Wesley. pg. 562
+ */
+ template <class T>
+ T
+ adaptsimp(T (*fun)(const T &), const T &a, const T &b, const int &N,
+ const T &tol)
+ {
+ if (a > b)
+ throw scythe_invalid_arg (__FILE__, __PRETTY_FUNCTION__, __LINE__,
+ "Lower limit larger than upper");
+
+ if (N <= 0)
+ throw scythe_invalid_arg (__FILE__, __PRETTY_FUNCTION__, __LINE__,
+ "Number of subintervals negative");
+
+ T I = intsimp(fun, a, b, N);
+ if (::fabs(I - intsimp(fun, a, b, N / 2)) > tol)
+ return adaptsimp(fun, a, (a + b) / 2, N, tol)
+ + adaptsimp(fun, (a + b) / 2, b, N, tol);
+
+ return I;
+ }
+
+
+ /* Numerically calculates the first derivative of a function
+ * Notes: Algorithm taken from Nocedal and Wright. 1999. section 7.1
+ * with some additional tricks from Press et al. NRC
+ */
+ template <class T>
+ Matrix<T>
+ gradfdif (T (*fun)(const Matrix<T> &, const Matrix<T> &,
+ const Matrix<T> &), const Matrix<T> &theta,
+ const Matrix<T> &y, const Matrix<T> &X)
+ {
+ if (! theta.isColVector())
+ throw scythe_dimension_error (__FILE__, __PRETTY_FUNCTION__,
+ __LINE__, "Theta not column vector");
+
+ int k = theta.size();
+ // stepsize CAREFUL-- THIS IS MACHINE-SPECIFIC!!!
+ T h = std::sqrt(epsilon<T>()); // 2.25e-16
+ //T h = std::sqrt(2.25e-16);
+
+
+ Matrix<T> grad(k,1);
+
+ for (int i = 0; i < k; ++i){
+ Matrix<T> e(k,1);
+ e[i] = h;
+ Matrix<T> temp = theta + e;
+ donothing(temp);
+ e = temp - theta;
+ grad[i] = (fun(theta + e, y, X) - fun(theta, y, X)) / e[i];
+ }
+
+ return grad;
+ }
+
+
+ /* Numerically calculates the first derivative of a function
+ * Notes: Algorithm taken from Nocedal and Wright. 1999. section 7.1
+ * with some additional tricks from Press et al. NRC
+ */
+ template <class T>
+ T
+ gradfdifls (T (*fun)(const Matrix<T> &, const Matrix<T> &,
+ const Matrix<T> &), const T &alpha,
+ const Matrix<T> &theta, const Matrix<T> &p,
+ const Matrix<T> &y, const Matrix<T> &X)
+ {
+ if (! theta.isColVector())
+ throw scythe_dimension_error (__FILE__, __PRETTY_FUNCTION__,
+ __LINE__, "Theta not column vector");
+ if (! p.isColVector())
+ throw scythe_dimension_error (__FILE__, __PRETTY_FUNCTION__,
+ __LINE__, "p not column vector");
+
+ int k = theta.size();
+ // stepsize CAREFUL-- THIS IS MACHINE-SPECIFIC!!!
+ T h = std::sqrt(epsilon<T>()); //2.2e-16
+ //T h = std::sqrt(2.2e-16);
+
+ T deriv;
+
+ for (int i = 0; i < k; ++i) {
+ T temp = alpha + h;
+ donothing(temp);
+ T e = temp - alpha;
+ deriv = (fun(theta + (alpha + e) * p, y, X)
+ - fun(theta + alpha * p, y, X)) / e;
+ }
+
+ return deriv;
+ }
+
+
+
+ /* Numerically calculates the gradient of a function
+ * Notes: Algorithm taken from Nocedal and Wright. 1999. section 7.1
+ * with some additional tricks from Press et al. NRC
+ */
+ template <class T>
+ Matrix<T>
+ jacfdif(Matrix<T> (*fun)(const Matrix<T> &), const Matrix<T> &theta)
+ {
+ if (! theta.isColVector())
+ throw scythe_dimension_error (__FILE__, __PRETTY_FUNCTION__,
+ __LINE__, "Theta not column vector");
+
+ Matrix<T> fval = fun(theta);
+
+ int k = theta.rows();
+ int n = fval.rows();
+ // stepsize CAREFUL -- THIS IS MACHINE-SPECIFIC!!!
+ T h = std::sqrt(epsilon<T>()); //2.2e-16
+ //T h = std::sqrt(2.2e-16);
+ Matrix<T> J(n,k);
+
+ for (int i = 0; i < k; ++i) {
+ Matrix<T> e(k,1);
+ e[i] = h;
+ Matrix<T> temp = theta + e;
+ donothing(temp);
+ e = temp - theta;
+ Matrix<T> fthetae = fun(theta + e);
+ Matrix<T> ftheta = fun(theta);
+ for (int j = 0; j < n; ++j) {
+ J(j,i) = (fthetae[j] - ftheta[j]) / e[i];
+ }
+ }
+
+ return J;
+ }
+
+ /* Numerically calculates the gradient of a function
+ * Notes: Algorithm taken from Nocedal and Wright. 1999. section 7.1
+ * with some additional tricks from Press et al. NRC
+ */
+ template <class T>
+ Matrix<T>
+ jacfdif(Matrix<T> (*fun)(const Matrix<T> &, const Matrix<T> &),
+ const Matrix<T> &theta, const Matrix<T> &psi)
+ {
+ if (! theta.isColVector())
+ throw scythe_dimension_error (__FILE__, __PRETTY_FUNCTION__,
+ __LINE__, "Theta not column vector");
+
+ Matrix<T> fval = fun(theta, psi);
+
+ int k = theta.rows();
+ int n = fval.rows();
+ // stepsize CAREFUL -- THIS IS MACHINE-SPECIFIC!!!
+ //T h = std::sqrt(epsilon<T>()); //2.2e-16
+ T h = std::sqrt(2.2e-16);
+ Matrix<T> J(n, k);
+
+ for (int i = 0; i < k; ++i) {
+ Matrix<T> e(k,1);
+ e[i] = h;
+ Matrix<T> temp = theta + e;
+ donothing(temp);
+ e = temp - theta;
+ Matrix<T> fthetae = fun(theta + e, psi);
+ Matrix<T> ftheta = fun(theta, psi);
+ for (int j = 0; j < n; ++j) {
+ J(j,i) = (fthetae[j] - ftheta[j]) / e[i];
+ }
+ }
+
+ return J;
+ }
+
+ /* Numerically calculates the Hessian of a function */
+ template <class T>
+ Matrix<T>
+ hesscdif (T (*fun)(const Matrix<T> &, const Matrix<T> &,
+ const Matrix<T> &), const Matrix<T> &theta,
+ const Matrix<T> &y, const Matrix<T> &X)
+ {
+ if (! theta.isColVector())
+ throw scythe_dimension_error (__FILE__, __PRETTY_FUNCTION__,
+ __LINE__, "Theta not column vector");
+
+ T fval = fun(theta,y,X);
+
+ int k = theta.rows();
+
+ // stepsize CAREFUL -- THIS IS MACHINE-SPECIFIC!!!
+ T h2 = (T) 1e-10;
+ T h = ::sqrt(h2);
+ Matrix<T> H(k,k);
+
+ for (int i=0; i<k; ++i) {
+ Matrix<T> ei(k, 1);
+ ei[i] = h;
+ Matrix<T> temp = theta + ei;
+ donothing(temp);
+ ei = temp - theta;
+ for (int j = 0; j < k; ++j){
+ Matrix<T> ej(k,1);
+ ej[j] = h;
+ temp = theta + ej;
+ donothing(temp);
+ ej = temp - theta;
+
+ if (i==j){
+ H(i,i) = ( -fun(theta + 2.0 * ei, y, X) + 16.0 *
+ fun(theta+ei,y,X) - 30.0 * fval + 16.0 *
+ fun(theta-ei,y,X) -
+ fun(theta-2.0 * ei, y, X)) / (12.0 * h2);
+ } else {
+ H(i,j) = ( fun(theta + ei + ej, y, X) - fun(theta+ei-ej, y, X)
+ - fun(theta - ei + ej, y, X) + fun(theta-ei-ej, y, X))
+ / (4.0 * h2);
+ }
+ }
+ }
+
+ return H;
+ }
+
+ /* Performs a linesearch to find the step length (\a alpha)
+ * Notes: Algorithm taken from Nocedal and Wright. 1999. Procedure
+ * 3.1
+ */
+ template <class T>
+ T
+ linesearch1(T (*fun)(const Matrix<T> &, const Matrix<T> &,
+ const Matrix<T> &), const Matrix<T> &theta,
+ const Matrix<T> &p, const Matrix<T> &y,
+ const Matrix<T> &X)
+ {
+ if (! theta.isColVector())
+ throw scythe_dimension_error (__FILE__, __PRETTY_FUNCTION__,
+ __LINE__, "Theta not column vector");
+ if (! p.isColVector())
+ throw scythe_dimension_error (__FILE__, __PRETTY_FUNCTION__,
+ __LINE__, "p not column vector");
+
+ T alpha_bar = (T) 1.0;
+ T rho = (T) 0.9;
+ T c = (T) 0.5;
+ T alpha = alpha_bar;
+ Matrix<T> fgrad = gradfdif(fun, theta, y, X);
+
+ while (fun(theta + alpha * p, y, X) > (fun(theta, y, X) + c
+ * alpha * t(fgrad) * p)[0]) {
+ alpha = rho * alpha;
+ }
+
+ return alpha;
+ }
+
+ /* Performs a linesearch to find the step length (\a alpha)
+ * Notes: Algorithm taken from Nocedal and Wright. 1999. Algorithm
+ * 3.2
+ */
+ template <class T>
+ T
+ linesearch2(T (*fun)(const Matrix<T> &, const Matrix<T> &,
+ const Matrix<T> &), const Matrix<T> &theta,
+ const Matrix<T> &p, const Matrix<T> &y,
+ const Matrix<T> &X, rng *myrng)
+ {
+ if (! theta.isColVector())
+ throw scythe_dimension_error (__FILE__, __PRETTY_FUNCTION__,
+ __LINE__, "Theta not column vector");
+ if (! p.isColVector())
+ throw scythe_dimension_error (__FILE__, __PRETTY_FUNCTION__,
+ __LINE__, "p not column vector");
+
+ T alpha_last = (T) 0.0;
+ T alpha_cur = (T) 1.0;
+ T alpha_max = (T) 10.0;
+ T c1 = (T) 1e-4;
+ T c2 = (T) 0.5;
+ int max_iter = 50;
+ T fgradalpha0 = gradfdifls(fun, (T) 0, theta, p, y, X);
+
+ for (int i = 0; i < max_iter; ++i) {
+ T phi_cur = fun(theta + alpha_cur * p, y, X);
+ T phi_last = fun(theta + alpha_last * p, y, X);
+
+ if ((phi_cur > (fun(theta, y, X) + c1 * alpha_cur * fgradalpha0))
+ ||
+ ((phi_cur >= phi_last) && (i > 0))) {
+ T alphastar = zoom(fun, alpha_last, alpha_cur, theta, p, y, X);
+ return alphastar;
+ }
+
+ T fgradalpha_cur = gradfdifls(fun, alpha_cur, theta, p, y, X);
+ if ( ::fabs(fgradalpha_cur) <= -1 * c2 * fgradalpha0)
+ return alpha_cur;
+
+ if ( fgradalpha_cur >= (T) 0.0) {
+ T alphastar = zoom(fun, alpha_cur, alpha_last, theta, p, y, X);
+ return alphastar;
+ }
+
+ alpha_last = alpha_cur;
+ alpha_cur = myrng->runif() * (alpha_max - alpha_cur) + alpha_cur;
+ }
+
+ return 0.001;
+ }
+
+
+ /* Finds the minimum of a function once bracketed (i.e. over a
+ * closed interval).
+ * Notes: Algorithm taken from Nocedal and Wright. 1999. Algorithm
+ * 3.3
+ */
+ template <class T>
+ T
+ zoom (T (*fun)(const Matrix<T> &, const Matrix<T> &,
+ const Matrix<T> &), const T &alo, const T &ahi,
+ const Matrix<T> &theta, const Matrix<T> &p,
+ const Matrix<T> &y, const Matrix<T> &X )
+ {
+ if (! theta.isColVector())
+ throw scythe_dimension_error (__FILE__, __PRETTY_FUNCTION__,
+ __LINE__, "Theta not column vector");
+ if (! p.isColVector())
+ throw scythe_dimension_error (__FILE__, __PRETTY_FUNCTION__,
+ __LINE__, "p not column vector");
+
+ T alpha_lo = alo;
+ T alpha_hi = ahi;
+ T alpha_j = (alo + ahi) / 2.0;
+ T phi_0 = fun(theta, y, X);
+ T c1 = (T) 1e-4;
+ T c2 = (T) 0.5;
+ T fgrad0 = gradfdifls(fun,(T) 0, theta, p, y, X);
+
+ int count = 0;
+ int maxit = 20;
+ while(count < maxit) {
+ T phi_j = fun(theta + alpha_j * p, y, X);
+ T phi_lo = fun(theta + alpha_lo * p, y, X);
+
+ if ((phi_j > (phi_0 + c1 * alpha_j * fgrad0))
+ || (phi_j >= phi_lo)){
+ alpha_hi = alpha_j;
+ } else {
+ T fgradj = gradfdifls(fun, alpha_j, theta, p, y, X);
+ if (::fabs(fgradj) <= -1 * c2 * fgrad0){
+ return alpha_j;
+ }
+ if ( fgradj * (alpha_hi - alpha_lo) >= 0){
+ alpha_hi = alpha_lo;
+ }
+ alpha_lo = alpha_j;
+ }
+ ++count;
+ }
+
+ return alpha_j;
+ }
+
+
+
+ /* Find the minimum of a function using the BFGS algorithm
+ * Notes: Algorithm taken from Nocedal and Wright. 1999. algorithm
+ * 8.1
+ */
+ template <class T>
+ Matrix<T>
+ BFGS (T (*fun)(const Matrix<T> &, const Matrix<T> &,
+ const Matrix<T> &), rng *myrng,
+ const Matrix<T> &theta,
+ const Matrix<T> &y, const Matrix<T> &X,
+ const int &maxit, const T &tolerance)
+ {
+ if (! theta.isColVector())
+ throw scythe_dimension_error (__FILE__, __PRETTY_FUNCTION__,
+ __LINE__, "Theta not column vector");
+
+ int n = theta.size();
+
+ // H is initial inverse hessian
+ Matrix<T> H = inv(hesscdif(fun, theta, y, X));
+
+ // gradient at starting values
+ Matrix<T> fgrad = gradfdif(fun, theta, y, X);
+ Matrix<T> thetamin = theta;
+ Matrix<T> fgrad_new = fgrad;
+ Matrix<T> I = eye<T>(n);
+
+ int count = 0;
+ while( (t(fgrad_new)*fgrad_new)[0] > tolerance) {
+ Matrix<T> p = -1 * H * fgrad;
+ T alpha = linesearch2(fun, thetamin, p, y, X, myrng);
+ Matrix<T> thetamin_new = thetamin + alpha*p;
+ fgrad_new = gradfdif(fun, thetamin_new, y, X);
+ Matrix<T> s = thetamin_new - thetamin;
+ Matrix<T> y = fgrad_new - fgrad;
+ T rho = 1.0 / (t(y) * s)[0];
+ H = (I - rho * s * t(y)) * H *(I - rho * y * t(s))
+ + rho * s * (!s);
+
+ thetamin = thetamin_new;
+ fgrad = fgrad_new;
+ ++count;
+
+ std::cout << "BFGS iteration = " << count << std::endl;
+ std::cout << "thetamin = " << (!thetamin).toString() << std::endl;
+ std::cout << "gradient = " << (!fgrad).toString() << std::endl;
+ std::cout << "t(gradient) * gradient = " <<
+ ((!fgrad) * fgrad).toString() << std::endl;
+
+ if (count > maxit)
+ throw scythe_convergence_error (__FILE__, __PRETTY_FUNCTION__,
+ __LINE__, "Failed to converge. Try better starting values");
+ }
+
+ return thetamin;
+ }
+
+
+ /* Zero a function using Broyen's Method
+ * Notes: Algorithm taken from Nocedal and Wright. 1999. algorithm 11
+ * line search is not used to determine alpha (this should probably
+ * be changed at some point.
+ */
+ template <class T>
+ Matrix<T>
+ nls_broyden(Matrix<T> (*fun)(const Matrix<T> &),
+ const Matrix<T> &theta, const int &maxit = 5000,
+ const T &tolerance = 1e-6)
+ {
+ if (! theta.isColVector())
+ throw scythe_dimension_error (__FILE__, __PRETTY_FUNCTION__,
+ __LINE__, "Theta not column vector");
+
+ Matrix<T> thetastar = theta;
+ Matrix<T> B = jacfdif(fun, thetastar);
+
+ Matrix<T> fthetastar;
+ Matrix<T> p;
+ Matrix<T> thetastar_new;
+ Matrix<T> fthetastar_new;
+ Matrix<T> s;
+ Matrix<T> y;
+
+ for (int i = 0; i < maxit; ++i) {
+ fthetastar = fun(thetastar);
+ p = lu_solve(B, -1 * fthetastar);
+ T alpha = (T) 1.0;
+ thetastar_new = thetastar + alpha*p;
+ fthetastar_new = fun(thetastar_new);
+ s = thetastar_new - thetastar;
+ y = fthetastar_new - fthetastar;
+ B = B + ((y - B * s) * (!s)) / ((!s) * s);
+ thetastar = thetastar_new;
+ if (max(fabs(fthetastar_new)) < tolerance)
+ return thetastar;
+ }
+
+ throw scythe_convergence_error (__FILE__, __PRETTY_FUNCTION__,
+ __LINE__,std::string("Failed to converge. Try better starting") &
+ " values or increase maxit");
+ }
+
+
+ /* Zero a function using Broyen's Method
+ * Notes: Algorithm taken from Nocedal and Wright. 1999. algorithm
+ * 11.3
+ * line search is not used to determine alpha (this should probably
+ * be changed at some point.
+ */
+ template <class T>
+ Matrix<T>
+ nls_broyden(Matrix<T> (*fun)(const Matrix<T> &, const Matrix<T> &),
+ const Matrix<T> &theta, const Matrix<T> &psi,
+ const int &maxit=5000, const T &tolerance=1e-6)
+ {
+ if (! theta.isColVector())
+ throw scythe_dimension_error (__FILE__, __PRETTY_FUNCTION__,
+ __LINE__, "Theta not column vector");
+ if (! psi.isColVector())
+ throw scythe_dimension_error (__FILE__, __PRETTY_FUNCTION__,
+ __LINE__, "Psi not column vector");
+
+ Matrix<T> thetastar = theta;
+ Matrix<T> B = jacfdif(fun, thetastar, psi);
+
+ Matrix<T> fthetastar;
+ Matrix<T> p;
+ Matrix<T> thetastar_new;
+ Matrix<T> fthetastar_new;
+ Matrix<T> s;
+ Matrix<T> y;
+
+ for (int i = 0; i < maxit; ++i) {
+ fthetastar = fun(thetastar, psi);
+ p = lu_solve(B, -1 * fthetastar);
+ T alpha = (T) 1.0;
+ thetastar_new = thetastar + alpha*p;
+ fthetastar_new = fun(thetastar_new, psi);
+ s = thetastar_new - thetastar;
+ y = fthetastar_new - fthetastar;
+ B = B + ((y - B * s) * (!s)) / ((!s) * s);
+ thetastar = thetastar_new;
+ if (max(fabs(fthetastar_new)) < tolerance)
+ return thetastar;
+ }
+
+ throw scythe_convergence_error (__FILE__, __PRETTY_FUNCTION__,
+ __LINE__,std::string("Failed to converge. Try better starting") &
+ " values or increase maxit");
+ }
+
+} // namespace dec
+
+#ifndef SCYTHE_COMPILE_DIRECT
+#include "scythestat/eti/optimize.t"
+#endif
+
+#endif /* SCYTHE_OPTIMIZE_CC */
diff --git a/src/rng.cc b/src/rng.cc
index 4d0fc72..4a2aa97 100644
--- a/src/rng.cc
+++ b/src/rng.cc
@@ -43,81 +43,81 @@
#endif
namespace SCYTHE {
-
- /* Default constructor */
- rng::rng ()
- {
- }
-
- rng::~rng()
- {
- }
-
- /* Random Numbers */
- Matrix<double>
- rng::runif (const int &rows, const int &cols)
- {
- int size = rows * cols;
- if (size <= 0) {
- throw scythe_invalid_arg (__FILE__, __PRETTY_FUNCTION__,
+
+ /* Default constructor */
+ rng::rng ()
+ {
+ }
+
+ rng::~rng()
+ {
+ }
+
+ /* Random Numbers */
+ Matrix<double>
+ rng::runif (const int &rows, const int &cols)
+ {
+ int size = rows * cols;
+ if (size <= 0) {
+ throw scythe_invalid_arg (__FILE__, __PRETTY_FUNCTION__,
__LINE__, "Attempted to create Matrix of size <= 0");
- }
-
- Matrix<double> temp(rows, cols, false);
- for (int i = 0; i < size; ++i)
- temp[i] = runif();
-
- return temp;
- }
-
- double
- rng::rbeta (const double& alpha, const double& beta)
- {
- static double report;
- double xalpha, xbeta;
-
- // Check for allowable parameters
- if (alpha <= 0) {
- throw scythe_invalid_arg(__FILE__, __PRETTY_FUNCTION__,
- __LINE__, "alpha <= 0");
- }
- if (beta <= 0) {
- throw scythe_invalid_arg(__FILE__, __PRETTY_FUNCTION__,
- __LINE__, "beta <= 0");
- }
-
- xalpha = rchisq (2 * alpha);
- xbeta = rchisq (2 * beta);
- report = xalpha / (xalpha + xbeta);
-
- return (report);
- }
-
- Matrix<double>
- rng::rbeta (const int& rows, const int& cols, const double& alpha,
- const double& beta)
- {
- int size = rows * cols;
- if (size <= 0) {
- throw scythe_invalid_arg (__FILE__, __PRETTY_FUNCTION__,
+ }
+
+ Matrix<double> temp(rows, cols, false);
+ for (int i = 0; i < size; ++i)
+ temp[i] = runif();
+
+ return temp;
+ }
+
+ double
+ rng::rbeta (const double& alpha, const double& beta)
+ {
+ static double report;
+ double xalpha, xbeta;
+
+ // Check for allowable parameters
+ if (alpha <= 0) {
+ throw scythe_invalid_arg(__FILE__, __PRETTY_FUNCTION__,
+ __LINE__, "alpha <= 0");
+ }
+ if (beta <= 0) {
+ throw scythe_invalid_arg(__FILE__, __PRETTY_FUNCTION__,
+ __LINE__, "beta <= 0");
+ }
+
+ xalpha = rchisq (2 * alpha);
+ xbeta = rchisq (2 * beta);
+ report = xalpha / (xalpha + xbeta);
+
+ return (report);
+ }
+
+ Matrix<double>
+ rng::rbeta (const int& rows, const int& cols, const double& alpha,
+ const double& beta)
+ {
+ int size = rows * cols;
+ if (size <= 0) {
+ throw scythe_invalid_arg (__FILE__, __PRETTY_FUNCTION__,
__LINE__, "Attempted to create Matrix of size <= 0");
- }
-
- Matrix<double> temp(rows, cols, false);
- for (int i = 0; i < size; ++i)
- temp[i] = rbeta (alpha, beta);
-
- return temp;
- }
-
+ }
+
+ Matrix<double> temp(rows, cols, false);
+ for (int i = 0; i < size; ++i)
+ temp[i] = rbeta (alpha, beta);
+
+ return temp;
+ }
+
/* Return a pseudo-random deviate from a non-cental hypergeometric
* distribution
*/
double
- rng::rnchypgeom(const double& m1, const double& n1,
- const double& n2, const double& psi,
- const double& delta)
- {
+ rng::rnchypgeom(const double& m1, const double& n1,
+ const double& n2, const double& psi,
+ const double& delta)
+ {
// Calculate mode of mass function
double a = psi - 1;
double b = -1 * ((n1+m1+2)*psi + n2 - m1);
@@ -226,24 +226,24 @@ namespace SCYTHE {
exit(500000);
}
- Matrix<double>
- rng::rnchypgeom(const int &rows, const int &cols, const double &m1,
- const double &n1, const double &n2,
- const double &psi, const double &delta)
- {
- int size = rows * cols;
- if (size <= 0) {
- throw scythe_invalid_arg (__FILE__, __PRETTY_FUNCTION__,
+ Matrix<double>
+ rng::rnchypgeom(const int &rows, const int &cols, const double &m1,
+ const double &n1, const double &n2,
+ const double &psi, const double &delta)
+ {
+ int size = rows * cols;
+ if (size <= 0) {
+ throw scythe_invalid_arg (__FILE__, __PRETTY_FUNCTION__,
__LINE__, "Attempted to create Matrix of size <= 0");
- }
-
- Matrix<double> temp(rows, cols, false);
- for (int i = 0; i < size; ++i)
- temp[i] = rnchypgeom (m1, n1, n2, psi, delta);
-
- return temp;
- }
-
+ }
+
+ Matrix<double> temp(rows, cols, false);
+ for (int i = 0; i < size; ++i)
+ temp[i] = rnchypgeom (m1, n1, n2, psi, delta);
+
+ return temp;
+ }
+
/* Random Numbers */
int
rng::rbinom (const int& n, const double& p)
@@ -853,9 +853,11 @@ namespace SCYTHE {
double
rng::rtnorm(const double& m, const double& v, const double& below,
- const double& above)
+ const double& above)
{
if (below > above) {
+ std::cout << "mean = " << m << " and var = " << v << std::endl <<
+ "below = " << below << " and above = " << above << std::endl;
throw scythe_invalid_arg (__FILE__, __PRETTY_FUNCTION__,
__LINE__, \
"Truncation bound not logically consistent");
@@ -896,9 +898,9 @@ namespace SCYTHE {
double
rng::rtnorm_combo(const double& m, const double& v,
- const double& below, const double& above)
+ const double& below, const double& above)
{
- if (below > above) {
+ if (below >= above) {
throw scythe_invalid_arg (__FILE__, __PRETTY_FUNCTION__,
__LINE__, "Truncation bound not logically consistent");
}
@@ -947,7 +949,7 @@ namespace SCYTHE {
double
rng::rtbnorm_slice (const double& m, const double& v,
- const double& below, const int& iter)
+ const double& below, const int& iter)
{
if (below < m) {
throw scythe_invalid_arg (__FILE__, __PRETTY_FUNCTION__,
@@ -977,7 +979,7 @@ namespace SCYTHE {
double
rng::rtanorm_slice (const double& m, const double& v,
- const double& above, const int& iter)
+ const double& above, const int& iter)
{
if (above > m) {
throw scythe_invalid_arg (__FILE__, __PRETTY_FUNCTION__,
@@ -1010,13 +1012,13 @@ namespace SCYTHE {
double
rng::rtbnorm_combo (const double& m, const double& v,
- const double& below, const int& iter)
+ const double& below, const int& iter)
{
if (v <= 0){
throw scythe_invalid_arg (__FILE__, __PRETTY_FUNCTION__,
__LINE__, "Variance non-positive");
}
-
+
double s = ::sqrt(v);
// do rejection sampling and return value
//if (m >= below){
@@ -1027,7 +1029,7 @@ namespace SCYTHE {
return x;
} else if ((m/s - below/s ) > -5.0 ){
// use the inverse cdf method
- double above = (m+30.0)*s;
+ double above = std::numeric_limits<double>::infinity();
double x = rtnorm(m, v, below, above);
return x;
} else {
@@ -1051,7 +1053,7 @@ namespace SCYTHE {
double
rng::rtanorm_combo (const double& m, const double& v,
- const double& above, const int& iter)
+ const double& above, const int& iter)
{
if (v <= 0){
throw scythe_invalid_arg (__FILE__, __PRETTY_FUNCTION__,
@@ -1066,7 +1068,7 @@ namespace SCYTHE {
return x;
} else if ((m/s - above/s ) < 5.0 ){
// use the inverse cdf method
- double below = (m-30.0)*s;
+ double below = -std::numeric_limits<double>::infinity();
double x = rtnorm(m, v, below, above);
return x;
} else {
diff --git a/src/smath.cc b/src/smath.cc
index 7628485..19cb34a 100644
--- a/src/smath.cc
+++ b/src/smath.cc
@@ -47,16 +47,16 @@ namespace SCYTHE {
return A;
}
- /* calc the inverse hyperbolic cosine of each element of a Matrix */
- template <class T>
- Matrix<T>
- acosh (Matrix<T> A)
- {
- for (int i = 0; i < A.size(); ++i)
- A[i] = ::acosh(A[i]);
-
- return A;
- }
+// /* calc the inverse hyperbolic cosine of each element of a Matrix */
+// template <class T>
+// Matrix<T>
+// acosh (Matrix<T> A)
+// {
+// for (int i = 0; i < A.size(); ++i)
+// A[i] = ::acosh(A[i]);
+//
+// return A;
+// }
/* calc the inverse sine of each element of a Matrix */
template <class T>
@@ -69,16 +69,16 @@ namespace SCYTHE {
return A;
}
- /* calc the inverse hyperbolic sine of each element of a Matrix */
- template <class T>
- Matrix<T>
- asinh (Matrix<T> A)
- {
- for (int i = 0; i < A.size(); ++i)
- A[i] = ::asinh(A[i]);
-
- return A;
- }
+// /* calc the inverse hyperbolic sine of each element of a Matrix */
+// template <class T>
+// Matrix<T>
+// asinh (Matrix<T> A)
+// {
+// for (int i = 0; i < A.size(); ++i)
+// A[i] = ::asinh(A[i]);
+//
+// return A;
+// }
/* calc the inverse tangent of each element of a Matrix */
template <class T>
@@ -91,16 +91,16 @@ namespace SCYTHE {
return A;
}
- /* calc the inverse hyperbolic tangent of each element of a Matrix */
- template <class T>
- Matrix<T>
- atanh (Matrix<T> A)
- {
- for (int i = 0; i < A.size(); ++i)
- A[i] = ::atanh(A[i]);
-
- return A;
- }
+// /* calc the inverse hyperbolic tangent of each element of a Matrix */
+// template <class T>
+// Matrix<T>
+// atanh (Matrix<T> A)
+// {
+// for (int i = 0; i < A.size(); ++i)
+// A[i] = ::atanh(A[i]);
+//
+// return A;
+// }
/* calc the angle whose tangent is y/x */
template <class T>
@@ -235,16 +235,16 @@ namespace SCYTHE {
return A;
}
- /* calc the exponent - 1 of each element of a Matrix */
- template <class T>
- Matrix<T>
- expm1 (Matrix<T> A)
- {
- for (int i = 0; i < A.size(); ++i)
- A[i] = ::expm1(A[i]);
-
- return A;
- }
+// /* calc the exponent - 1 of each element of a Matrix */
+// template <class T>
+// Matrix<T>
+// expm1 (Matrix<T> A)
+// {
+// for (int i = 0; i < A.size(); ++i)
+// A[i] = ::expm1(A[i]);
+//
+// return A;
+// }
/* calc the absval of each element of a Matrix */
template <class T>
diff --git a/src/smath.h b/src/smath.h
index 6703456..1e67050 100644
--- a/src/smath.h
+++ b/src/smath.h
@@ -47,25 +47,25 @@ namespace SCYTHE {
template <class T>
Matrix<T> acos (Matrix<T>);
- /* acosh - inverse hyperbolic cosine function */
- template <class T>
- Matrix<T> acosh (Matrix<T>);
+// /* acosh - inverse hyperbolic cosine function */
+// template <class T>
+// Matrix<T> acosh (Matrix<T>);
/* asin - inverse sine function */
template <class T>
Matrix <T> asin (Matrix<T>);
- /* asinh - inverse hyperbolic sine function */
- template <class T>
- Matrix<T> asinh (Matrix<T>);
+// /* asinh - inverse hyperbolic sine function */
+// template <class T>
+// Matrix<T> asinh (Matrix<T>);
/* atan - inverse tangent function */
template <class T>
Matrix<T> atan (Matrix<T>);
- /* atanh - inverse hyperbolic tangent function */
- template <class T>
- Matrix<T> atanh (Matrix<T>);
+// /* atanh - inverse hyperbolic tangent function */
+// template <class T>
+// Matrix<T> atanh (Matrix<T>);
/* atan2 - returns angle whose tangent is y/x in the full angular
* range [-pit,+pi]. Domain error if both x and y zero
@@ -110,9 +110,9 @@ namespace SCYTHE {
template <class T>
Matrix<T> exp (Matrix<T>);
- /* expm1 - exponent minus 1 */
- template <class T>
- Matrix<T> expm1 (Matrix<T>);
+// /* expm1 - exponent minus 1 */
+// template <class T>
+// Matrix<T> expm1 (Matrix<T>);
/* fabs - Calculate the absolute value of each Matrix element */
template <class T>
--
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