[r-cran-mcmcpack] 32/90: Imported Upstream version 0.9-6
Andreas Tille
tille at debian.org
Fri Dec 16 09:07:37 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 1fa77cdba8e4a55f948df0e660017bc3a7fd5ff7
Author: Andreas Tille <tille at debian.org>
Date: Fri Dec 16 08:07:13 2016 +0100
Imported Upstream version 0.9-6
---
DESCRIPTION | 8 +-
HISTORY | 15 ++
LICENSE | 14 +-
NAMESPACE | 3 +-
R/BayesFactors.R | 13 ++
R/MCMCSVDreg.R | 11 +-
R/MCMCdynamicEI.R | 19 ++-
R/MCMCfactanal.R | 9 +-
R/MCMChierEI.R | 19 ++-
R/MCMCirt1d.R | 10 ++
R/MCMCirtHier1d.R | 270 ++++++++++++++++++++++++++++++++
R/MCMCirtKd.R | 7 +
R/MCMCirtKdRob.R | 7 +
R/MCMClogit.R | 14 +-
R/MCMCmetrop1R.R | 10 ++
R/MCMCmixfactanal.R | 7 +
R/MCMCmnl.R | 12 +-
R/MCMCoprobit.R | 10 ++
R/MCMCordfactanal.R | 7 +
R/MCMCpoisson.R | 23 ++-
R/MCMCpoissonChangepoint.R | 10 ++
R/MCMCprobit.R | 14 +-
R/MCMCregress.R | 14 +-
R/MCMCtobit.R | 18 ++-
R/MCmodels.R | 14 ++
R/automate.R | 11 ++
R/btsutil.R | 20 ++-
R/distn.R | 21 ++-
R/hidden.R | 11 +-
R/procrust.R | 10 ++
R/scythe.R | 12 +-
R/tomog.R | 15 +-
R/utility.R | 13 +-
R/zzz.R | 12 ++
README | 5 +-
man/MCMCSVDreg.Rd | 2 +-
man/MCMCfactanal.Rd | 2 +-
man/MCMCirt1d.Rd | 2 +-
man/MCMCirtHier1d.Rd | 305 ++++++++++++++++++++++++++++++++++++
man/MCMClogit.Rd | 2 +-
man/MCMCmnl.Rd | 10 +-
man/MCMCpoisson.Rd | 2 +-
man/MCMCprobit.Rd | 2 +-
man/MCMCregress.Rd | 2 +-
man/MCMCtobit.Rd | 2 +-
man/PErisk.Rd | 2 +-
src/MCMCSVDreg.cc | 8 +-
src/MCMCdynamicEI.cc | 8 +
src/MCMCfactanal.cc | 8 +-
src/MCMCfcds.h | 93 ++++++++++-
src/MCMChierEI.cc | 6 +
src/MCMCirt1d.cc | 6 +
src/MCMCirtHier1d.cc | 350 ++++++++++++++++++++++++++++++++++++++++++
src/MCMCirtKdRob.cc | 8 +-
src/MCMClogit.cc | 10 +-
src/MCMClogituserprior.cc | 9 +-
src/MCMCmetrop1R.cc | 43 ++++--
src/MCMCmixfactanal.cc | 8 +-
src/MCMCmnl.h | 9 +-
src/MCMCmnlMH.cc | 9 +-
src/MCMCmnlslice.cc | 9 +-
src/MCMCoprobit.cc | 9 +-
src/MCMCordfactanal.cc | 8 +-
src/MCMCpoisson.cc | 9 +-
src/MCMCpoissonChangepoint.cc | 6 +-
src/MCMCprobit.cc | 8 +-
src/MCMCprobitres.cc | 9 +-
src/MCMCregress.cc | 8 +-
src/MCMCrng.h | 8 +-
src/MCMCtobit.cc | 9 +-
70 files changed, 1548 insertions(+), 131 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index 527a495..7953237 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,12 +1,12 @@
Package: MCMCpack
-Version: 0.9-4
-Date: 2008-3-3
+Version: 0.9-6
+Date: 2009-2-17
Title: Markov chain Monte Carlo (MCMC) Package
Author: Andrew D. Martin <admartin at wustl.edu>,
Kevin M. Quinn <kevin_quinn at harvard.edu>,
Jong Hee Park <jhp at uchicago.edu>
Maintainer: Andrew D. Martin <admartin at wustl.edu>
-Depends: R (>= 2.6.0), coda (>= 0.11-3), MASS, stats
+Depends: R (>= 2.8.0), coda (>= 0.11-3), MASS, stats
Description: This package contains functions to perform Bayesian
inference using posterior simulation for a number of statistical
models. Most simulation is done in compiled C++ written in the Scythe
@@ -19,4 +19,4 @@ Description: This package contains functions to perform Bayesian
License: GPL version 2
SystemRequirements: gcc (>= 4.0)
URL: http://mcmcpack.wustl.edu
-Packaged: Mon Mar 3 11:38:07 2008; adm
+Packaged: Tue Feb 17 10:39:42 2009; adm
diff --git a/HISTORY b/HISTORY
index 8de0702..d88800d 100644
--- a/HISTORY
+++ b/HISTORY
@@ -2,6 +2,21 @@
// Changes and Bug Fixes
//
+0.9-6 to 1.0-1
+ * added one-dimensional dynamic IRT model
+
+0.9-5 to 0.9-6
+ * fixed bug in MCMCmetrop1R.R [thanks to many for spotting memory issue]
+ * fixed by in Lecuyer seeds [thanks to Eduardo Leoni]
+
+0.9-4 to 0.9-5
+ * changed how missing values are handled when calculating agreement scores
+ for starting values in factor analysis and IRT models
+ * added warning about incomparable models to BayesFactor
+ * fixed bug with MCMCmnl example [thanks to Duncan Murdoch]
+ * fixed formula interface bug for all models [thanks to Duncan Murdoch]
+ * added hierarchical IRT code [written by Mike Malecki]
+
0.9-3 to 0.9-4
* fixed issue with MCMCmetrop1R.R [thanks to Jim Albert]
* fixed Linux compilation issue [thanks to Chris Lawrence]
diff --git a/LICENSE b/LICENSE
index cd38f59..6f9c100 100644
--- a/LICENSE
+++ b/LICENSE
@@ -1,5 +1,7 @@
Markov chain Monte Carlo Package (MCMCpack)
Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+ and Jong Hee Park
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
@@ -18,16 +20,14 @@ Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
Please contact:
Andrew D. Martin, Ph.D.
-Department of Political Science
-Washington University
-Campus Box 1063
-One Brookings Drive
-St. Louis, MO 63130
+Professor and Chair, Political Science, Arts & Sciences
+Professor and CERL Director, School of Law
+Washington University in St. Louis
+
(314) 935-5863 (Office)
-(314) 753-8377 (Cell)
(314) 935-5856 (Fax)
Email: admartin at wustl.edu
-WWW: http://adm.wustl.edu
+WWW: http://adm.wustl.edu
With any problems or questions.
diff --git a/NAMESPACE b/NAMESPACE
index aecee80..57c881e 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -17,10 +17,10 @@ export(
MCnormalnormal,
MCmultinomdirichlet,
MCMCdynamicEI,
-## MCMCdynamicIRT1d,
MCMCfactanal,
MCMChierEI,
MCMCirt1d,
+ MCMCirtHier1d,
MCMCirtKd,
MCMCirtKdRob,
MCMClogit,
@@ -51,6 +51,5 @@ export(
xpnd
)
-
S3method(print, BayesFactor)
S3method(summary, BayesFactor)
diff --git a/R/BayesFactors.R b/R/BayesFactors.R
index 7dc5b20..3720fca 100644
--- a/R/BayesFactors.R
+++ b/R/BayesFactors.R
@@ -1,8 +1,18 @@
+##########################################################################
## BayesFactor.R contains functions useful for calculating and comparing
## marginal likelihoods
##
+## 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.
+##
## Originally written by KQ 1/26/2006
##
+## Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+## Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+## and Jong Hee Park
+##########################################################################
+
## log densities
"logdinvgamma" <- function(sigma2, a, b){
@@ -152,6 +162,9 @@
attr(model.list[[j]], "logmarglike")
BF.mat[i,j] <- exp(BF.log.mat[i,j])
}
+ else{
+ warning(paste(model.names[i], " and ", model.names[j], " do not have exactly identical y data.\nBayes factors are not defined.\n", sep=""))
+ }
}
}
diff --git a/R/MCMCSVDreg.R b/R/MCMCSVDreg.R
index 61bfab2..e237970 100644
--- a/R/MCMCSVDreg.R
+++ b/R/MCMCSVDreg.R
@@ -1,3 +1,4 @@
+##########################################################################
## MCMCSVDreg.R samples from the posterior distribution of a Gaussian
## linear regression model in which the X matrix has been decomposed
## with an SVD. Useful for prediction when number of columns of X
@@ -6,8 +7,16 @@
## See West, Mike. 2000. "Bayesian Regression in the 'Large p, Small n'
## Paradigm." Duke ISDS Discussion Paper 2000-22.
##
+## 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.
+##
## KQ 9/9/2005
##
+## Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+## Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+## and Jong Hee Park
+##########################################################################
parse.formula.SVDreg <- function(formula, data, intercept){
@@ -47,7 +56,7 @@ parse.formula.SVDreg <- function(formula, data, intercept){
"MCMCSVDreg" <-
- function(formula, data=parent.frame(), burnin = 1000, mcmc = 10000,
+ function(formula, data=NULL, burnin = 1000, mcmc = 10000,
thin=1, verbose = 0, seed = NA, tau2.start = 1,
g0 = 0, a0 = 0.001, b0 = 0.001, c0=2, d0=2, w0=1,
beta.samp=FALSE, intercept=TRUE, ...) {
diff --git a/R/MCMCdynamicEI.R b/R/MCMCdynamicEI.R
index 6e80575..24201a9 100644
--- a/R/MCMCdynamicEI.R
+++ b/R/MCMCdynamicEI.R
@@ -1,7 +1,18 @@
-# sample from the posterior of Quinn's dynamic ecological inference model
-# in R using linked C++ code in Scythe
-#
-# KQ 10/25/2002
+##########################################################################
+## sample from the posterior of Quinn's dynamic ecological inference model
+## in R using linked C++ code in Scythe
+##
+## 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.
+##
+## KQ 10/25/2002
+##
+## Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+## Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+## and Jong Hee Park
+##########################################################################
+
"MCMCdynamicEI" <-
function(r0, r1, c0, c1, burnin=5000, mcmc=50000,
diff --git a/R/MCMCfactanal.R b/R/MCMCfactanal.R
index 27de0b5..b4cd928 100644
--- a/R/MCMCfactanal.R
+++ b/R/MCMCfactanal.R
@@ -19,14 +19,21 @@
## Kevin M. Quinn
## Harvard University
##
+## 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.
+##
## May 7, 2003
## revised to accomodate new spec 7/5/2004 KQ
##
+## Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+## Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+## and Jong Hee Park
##########################################################################
"MCMCfactanal" <-
function(x, factors, lambda.constraints=list(),
- data=parent.frame(), burnin = 1000, mcmc = 20000,
+ data=NULL, burnin = 1000, mcmc = 20000,
thin=1, verbose = 0, seed = NA,
lambda.start = NA, psi.start = NA,
l0=0, L0=0, a0=0.001, b0=0.001,
diff --git a/R/MCMChierEI.R b/R/MCMChierEI.R
index 688d969..9aa99c1 100644
--- a/R/MCMChierEI.R
+++ b/R/MCMChierEI.R
@@ -1,7 +1,18 @@
-# sample from the posterior distribution of Wakefield's hierarchical model
-# for ecological inference in R using linked C++ code in Scythe
-#
-# KQ 10/22/2002
+##########################################################################
+## sample from the posterior distribution of Wakefield's hierarchical model
+## for ecological inference in R using linked C++ code in Scythe
+##
+## 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.
+##
+## KQ 10/22/2002
+##
+## Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+## Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+## and Jong Hee Park
+##########################################################################
+
"MCMChierEI" <-
function(r0, r1, c0, c1, burnin=5000, mcmc=50000, thin=1,
diff --git a/R/MCMCirt1d.R b/R/MCMCirt1d.R
index 1591060..935aca7 100644
--- a/R/MCMCirt1d.R
+++ b/R/MCMCirt1d.R
@@ -1,9 +1,19 @@
+##########################################################################
## sample from the posterior distribution of a one-dimensional item
## response theory model in R using linked C++ code in Scythe.
##
+## 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.
+##
## ADM and KQ 1/23/2003
## updated extensively ADM & KQ 7/28/2004
## store.ability arg added KQ 1/27/2006
+##
+## Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+## Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+## and Jong Hee Park
+##########################################################################
"MCMCirt1d" <-
function(datamatrix, theta.constraints=list(), burnin = 1000,
diff --git a/R/MCMCirtHier1d.R b/R/MCMCirtHier1d.R
new file mode 100644
index 0000000..05dfd83
--- /dev/null
+++ b/R/MCMCirtHier1d.R
@@ -0,0 +1,270 @@
+## sample from the posterior distribution of a one-dimensional item
+## response theory model in R using linked C++ code in Scythe.
+##
+## ADM and KQ 1/23/2003
+## updated extensively ADM & KQ 7/28/2004
+## store.ability arg added KQ 1/27/2006
+## Hierarchical Subject Parameters MJM 2007-11-07
+## Parameter Expansion 2008-11-18
+## Grouped Subject Parameters (Wdata) started MJM,YS 2008-11
+
+"MCMCirtHier1d" <-
+ function(datamatrix, Xjdata,
+ burnin = 1000, mcmc = 20000, thin=1,
+ verbose = 0, seed = NA,
+ theta.start = NA, a.start = NA, b.start = NA,
+ beta.start=NA, b0=0, B0=.01, c0=.001, d0=.001,
+ ab0=0, AB0=.25, store.item = FALSE, store.ability=TRUE,
+ drop.constant.items=TRUE,
+ marginal.likelihood=c("none","Chib95"), px=TRUE,
+ px_a0 = 10, px_b0=10,
+ ... ) {
+
+ ## checks
+ check.offset(list(...))
+ check.mcmc.parameters(burnin, mcmc, thin)
+
+ ## check vote matrix and convert to work with C++ code
+ if (drop.constant.items==TRUE){
+ x.col.var <- apply(datamatrix, 2, var, na.rm=TRUE)
+ keep.inds <- x.col.var>0
+ keep.inds[is.na(keep.inds)] <- FALSE
+ datamatrix <- datamatrix[,keep.inds]
+ }
+ datamatrix <- as.matrix(datamatrix)
+ K <- ncol(datamatrix) # cases, bills, items, etc
+ J <- nrow(datamatrix) # justices, legislators, subjects, etc
+ L <- ncol(Xjdata) # predictors on theta Xj
+
+ if(sum(datamatrix==1 | datamatrix==0 | is.na(datamatrix)) != (J*K)) {
+ cat("Error: Data matrix contains elements other than 0, 1 or NA.\n")
+ stop("Please check data and call ", calling.function(), " again.\n",
+ call.=FALSE)
+ }
+ datamatrix[is.na(datamatrix)] <- 9
+ item.names <- colnames(as.data.frame(datamatrix))
+ subject.names <- rownames(as.data.frame(datamatrix))
+ beta.names <- c(names(as.data.frame(Xjdata)),"sigmasq")
+
+ ## names
+ item.names <- colnames(datamatrix)
+ if (is.null(item.names)){
+ item.names <- paste("item", 1:K, sep="")
+ }
+
+ ## check Xj matrix and set up betastart
+ Xjdata <- as.matrix(Xjdata)
+ if(nrow(Xjdata) != nrow(datamatrix)) {
+ cat("Error: subject covariates not of same length as datamatrix\n")
+ stop("Please check data and try ",calling.function()," again.\n",call.=FALSE)
+ }
+
+
+ ## prior for (alpha, beta)
+ holder <- form.mvn.prior(ab0, AB0, 2)
+ ab0 <- holder[[1]]
+ AB0 <- holder[[2]]
+
+ ## starting values for theta error checking
+ ## could use factor.score.start.check EXCEPT
+ ## We have done away with eq and ineq constraints.
+ if (max(is.na(theta.start))==1) {
+ theta.start <- factor.score.eigen.start(agree.mat(datamatrix), 1)
+ }
+ else if(is.numeric(theta.start) & length(theta.start) == J ) {
+ theta.start <- theta.start * matrix(1, J, 1)
+ }
+ else {
+ cat("Inappropriate value of theta.start passed.\n")
+ stop("Please respecify and call", calling.function(), " again.\n",
+ call.=FALSE)
+ }
+ ## starting values for (a, b)
+ ab.starts <- matrix(NA, K, 2)
+ for (i in 1:K){
+ local.y <- datamatrix[,i]
+ local.y[local.y==9] <- NA
+ if (var(na.omit(local.y))==0){
+ ab.starts[i,] <- c(0,10)
+ }
+ else {
+ ab.starts[i,] <- coef(suppressWarnings(glm(local.y~theta.start,
+ family=binomial(probit),
+ control=glm.control(
+ maxit=8, epsilon=1e-3)
+ )))
+ }
+ }
+ ab.starts[,1] <- -1 * ab.starts[,1] # make this into a difficulty param
+
+ ## starting values for a and b error checking
+ if (is.na(a.start)) {
+ a.start <- ab.starts[,1]
+ }
+ else if(is.null(dim(a.start))) {
+ a.start <- a.start * matrix(1,K,1)
+ }
+ else if((dim(a.start)[1] != K) || (dim(a.start)[2] != 1)) {
+ cat("Error: Starting value for a not conformable.\n")
+ stop("Please respecify and call ", calling.function(), " again.\n",
+ call.=FALSE)
+ }
+ if (is.na(b.start)) {
+ b.start <- ab.starts[,2]
+ }
+ else if(is.null(dim(b.start))) {
+ b.start <- b.start * matrix(1,K,1)
+ }
+ else if((dim(b.start)[1] != K) || (dim(b.start)[2] != 1)) {
+ cat("Error: Starting value for b not conformable.\n")
+ stop("Please respecify and call ", calling.function(), " again.\n",
+ call.=FALSE)
+ }
+
+ cat("doing betastart ")
+ ## starting values are regression of theta.start on Xj
+ ## or passed vector, or same value all beta
+ if (max(is.na(beta.start))==1) { # beta.start NA
+ beta.start <- coef(suppressWarnings(glm.fit(Xjdata,theta.start)))
+ } else if ( length(beta.start) == L ) { # beta.start vector
+ beta.start <- matrix(beta.start,L,1)
+ } else if ( length(beta.start) == 1 ) { # beta.start scalar
+ beta.start <- beta.start * matrix(1,L,1)
+ } else {
+ cat("Error: Starting value for beta not conformable.\n")
+ stop("Please respecify and call ", calling.function(), " again.\n",
+ call.=FALSE)
+ }
+ print(beta.start)
+ ## prior for beta
+ holder <- form.mvn.prior(b0, B0, L)
+ b0 <- holder[[1]]
+ B0 <- holder[[2]]
+ check.ig.prior(c0, d0)
+
+ ## get marginal likelihood argument
+ marginal.likelihood <- match.arg(marginal.likelihood)
+ B0.eigenvalues <- eigen(B0)$values
+ if (min(B0.eigenvalues) < 0){
+ stop("B0 is not positive semi-definite.\nPlease respecify and call again.\n")
+ }
+ if (isTRUE(all.equal(min(B0.eigenvalues), 0))){
+ if (marginal.likelihood != "none"){
+ warning("Cannot calculate marginal likelihood with improper prior\n")
+ marginal.likelihood <- "none"
+ }
+ }
+ logmarglike <- NULL
+ chib <- 0
+ if (marginal.likelihood == "Chib95"){
+ chib <- 1
+ }
+
+ cat("setting up posterior holder\n" )
+
+ ## define holder for posterior sample
+ if(store.item == FALSE & store.ability == TRUE) {
+ sample <- matrix(data=0, mcmc/thin, J+L+1)
+ }
+ else if (store.item == TRUE & store.ability == FALSE){
+ sample <- matrix(data=0, mcmc/thin, L+1 + 2*K)
+ }
+ else if (store.item == TRUE & store.ability == TRUE){
+ sample <- matrix(data=0, mcmc/thin, L+1 + J + 2 * K)
+ }
+ else{
+ stop("Either store.item or store.ability should be true.\n")
+ }
+
+ ## seeds
+ seeds <- form.seeds(seed)
+ lecuyer <- seeds[[1]]
+ seed.array <- seeds[[2]]
+ lecuyer.stream <- seeds[[3]]
+
+ # call C++ code to draw sample
+ posterior <- .C("MCMCirtHier1d",
+ sampledata = as.double(sample),
+ samplerow = as.integer(nrow(sample)),
+ samplecol = as.integer(ncol(sample)),
+ Xdata = as.integer(datamatrix),
+ Xrow = as.integer(nrow(datamatrix)),
+ Xcol = as.integer(ncol(datamatrix)),
+ 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),
+ thetastartdata = as.double(theta.start),
+ thetastartrow = as.integer(nrow(theta.start)),
+ thetastartcol = as.integer(ncol(theta.start)),
+ astartdata = as.double(a.start),
+ astartrow = as.integer(length(a.start)),
+ astartcol = as.integer(1),
+ bstartdata = as.double(b.start),
+ bstartrow = as.integer(length(b.start)),
+ bstartcol = as.integer(1),
+ ab0data = as.double(ab0),
+ ab0row = as.integer(nrow(ab0)),
+ ab0col = as.integer(ncol(ab0)),
+ AB0data = as.double(AB0),
+ AB0row = as.integer(nrow(AB0)),
+ AB0col = as.integer(ncol(AB0)),
+ Xjdata = as.double(Xjdata),
+ Xjrow = as.integer(nrow(Xjdata)),
+ Xjcol = as.integer(ncol(Xjdata)),
+ betastartdata = as.double(beta.start),
+ betastartrow = as.integer(length(beta.start)),
+ betastartcol = as.integer(1),
+ b0data = as.double(b0),
+ b0row = as.integer(length(b0)),
+ b0col = as.integer(1),
+ B0data = as.double(B0),
+ B0row = as.integer(nrow(B0)),
+ B0col = as.integer(ncol(B0)),
+ c0 = as.double(c0),
+ d0 = as.double(d0),
+ storei = as.integer(store.item),
+ storea = as.integer(store.ability),
+ logmarglikeholder = as.double(0.0),
+ chib = as.integer(chib),
+ px= as.integer(px),
+ px_a0 = as.double(px_a0),
+ px_b0 = as.double(px_b0),
+ PACKAGE="MCMCpack"
+ )
+
+ beta.names <- paste("beta.",beta.names,sep="")
+ theta.names <- paste("theta.", subject.names, sep = "")
+ alpha.beta.names <- paste(rep(c("alpha.","beta."), K),
+ rep(item.names, each = 2),
+ sep = "")
+
+ # put together matrix and build MCMC object to return
+ sample <- matrix(posterior$sampledata, posterior$samplerow,
+ posterior$samplecol,
+ byrow=FALSE)
+
+ output <- mcmc(data=sample, start=burnin+1, end=burnin+mcmc, thin=thin)
+ if (marginal.likelihood == "Chib95"){
+ logmarglike <- posterior$logmarglikeholder
+ }
+
+ names <- NULL
+ if(store.ability == TRUE) {
+ names <- c(names, theta.names)
+ }
+ if (store.item == TRUE){
+ names <- c(names, alpha.beta.names)
+ }
+ names <- c(names,beta.names)
+
+ try( varnames(output) <- names)
+ attr(output,"title") <-
+ "MCMCirtHier1d Posterior Sample"
+ attr(output,"logmarglike") <- posterior$logmarglikeholder
+ return(output)
+
+ }
diff --git a/R/MCMCirtKd.R b/R/MCMCirtKd.R
index 5869eef..f5c7513 100644
--- a/R/MCMCirtKd.R
+++ b/R/MCMCirtKd.R
@@ -9,8 +9,15 @@
## Kevin M. Quinn
## Harvard University
##
+## 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.
+##
## June 8, 2003
##
+## Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+## Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+## and Jong Hee Park
##########################################################################
"MCMCirtKd" <-
diff --git a/R/MCMCirtKdRob.R b/R/MCMCirtKdRob.R
index 2a1c879..903c7db 100644
--- a/R/MCMCirtKdRob.R
+++ b/R/MCMCirtKdRob.R
@@ -17,8 +17,15 @@
## Kevin M. Quinn
## Harvard University
##
+## 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.
+##
## Feb. 17, 2005
##
+## Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+## Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+## and Jong Hee Park
##########################################################################
diff --git a/R/MCMClogit.R b/R/MCMClogit.R
index 84f3c25..cac921a 100644
--- a/R/MCMClogit.R
+++ b/R/MCMClogit.R
@@ -1,17 +1,27 @@
+##########################################################################
## sample from the posterior distribution of a logistic regression
## model in R using linked C++ code in Scythe
##
## KQ 1/23/2003
##
+## 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.
+##
## Modified to meet new developer specification 7/15/2004 KQ
## Modified for new Scythe and rngs 7/25/2004 KQ
## note: B0 is now a precision
## Modified to allow user-specified prior density 8/17/2005 KQ
## Modified to handle marginal likelihood calculation 1/27/2006 KQ
##
+## Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+## Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+## and Jong Hee Park
+##########################################################################
+
"MCMClogit" <-
- function(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
+ function(formula, data=NULL, burnin = 1000, mcmc = 10000,
thin=1, tune=1.1, verbose = 0, seed = NA, beta.start = NA,
b0 = 0, B0 = 0, user.prior.density=NULL, logfun=TRUE,
marginal.likelihood = c("none", "Laplace"), ...) {
@@ -30,7 +40,7 @@
lecuyer.stream <- seeds[[3]]
## form response and model matrices
- holder <- parse.formula(formula, data)
+ holder <- parse.formula(formula, data=data)
Y <- holder[[1]]
X <- holder[[2]]
xnames <- holder[[3]]
diff --git a/R/MCMCmetrop1R.R b/R/MCMCmetrop1R.R
index 7106227..2073536 100644
--- a/R/MCMCmetrop1R.R
+++ b/R/MCMCmetrop1R.R
@@ -1,8 +1,13 @@
+##########################################################################
## samples from a user-written posterior coded in R using a
## random walk Metropolis algorithm
##
## KQ 6/24/2004
##
+## 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.
+##
## modified to work with non-invertible Hessian KQ 6/28/2005
##
## changed the method used to pass additional arguments to the user-defined
@@ -10,6 +15,11 @@
##
## changed to allow more user control of optim KQ 6/18/2006
##
+## Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+## Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+## and Jong Hee Park
+##########################################################################
+
"MCMCmetrop1R" <- function(fun, theta.init,
burnin=500, mcmc=20000, thin=1,
diff --git a/R/MCMCmixfactanal.R b/R/MCMCmixfactanal.R
index 25f662c..f686c77 100644
--- a/R/MCMCmixfactanal.R
+++ b/R/MCMCmixfactanal.R
@@ -20,10 +20,17 @@
## Kevin M. Quinn
## Harvard University
##
+## 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.
+##
## 12/2/2003
## Revised to accommodate new spec 7/20/2004
## Minor bug fix regarding std.mean 6/25/2004
##
+## Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+## Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+## and Jong Hee Park
##########################################################################
"MCMCmixfactanal" <-
diff --git a/R/MCMCmnl.R b/R/MCMCmnl.R
index 316f286..357c360 100644
--- a/R/MCMCmnl.R
+++ b/R/MCMCmnl.R
@@ -1,7 +1,17 @@
+##########################################################################
## MCMCmnl.R samples from the posterior distribution of a multinomial
## logit model using Metropolis-Hastings.
##
+## 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.
+##
## KQ 12/22/2004
+##
+## Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+## Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+## and Jong Hee Park
+##########################################################################
## parse formula and return a list that contains the model response
## matrix as element one, the model matrix as element two,
@@ -209,7 +219,7 @@
"MCMCmnl" <-
- function(formula, baseline=NULL, data = parent.frame(),
+ function(formula, baseline=NULL, data=NULL,
burnin = 1000, mcmc = 10000, thin=1,
mcmc.method = c("IndMH", "RWM", "slice"),
tune = 1.0, tdf=6, verbose = 0, seed = NA,
diff --git a/R/MCMCoprobit.R b/R/MCMCoprobit.R
index 33eb794..ec58c9e 100644
--- a/R/MCMCoprobit.R
+++ b/R/MCMCoprobit.R
@@ -1,9 +1,19 @@
+##########################################################################
## sample from the posterior distribution of an ordered probit model
## via the data augmentation approach of Cowles (1996)
##
+## 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.
+##
## KQ 1/25/2003
## Modified to meet new developer specification 7/26/2004 KQ
## Modified for new Scythe and rngs 7/26/2004 KQ
+##
+## Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+## Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+## and Jong Hee Park
+##########################################################################
"MCMCoprobit" <-
diff --git a/R/MCMCordfactanal.R b/R/MCMCordfactanal.R
index 12f94db..721ed63 100644
--- a/R/MCMCordfactanal.R
+++ b/R/MCMCordfactanal.R
@@ -18,9 +18,16 @@
## Kevin M. Quinn
## Harvard University
##
+## 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.
+##
## May 12, 2003
## Revised to accommodate new spec 7/13/2004
##
+## Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+## Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+## and Jong Hee Park
##########################################################################
"MCMCordfactanal" <-
diff --git a/R/MCMCpoisson.R b/R/MCMCpoisson.R
index bef353e..9ec6b55 100644
--- a/R/MCMCpoisson.R
+++ b/R/MCMCpoisson.R
@@ -1,14 +1,25 @@
-### sample from the posterior distribution of a Poisson regression
-### model in R using linked C++ code in Scythe
-###
-### ADM 1/24/2003
+##########################################################################
+## sample from the posterior distribution of a Poisson regression
+## model in R using linked C++ code in Scythe
+##
+## ADM 1/24/2003
## KQ 3/17/2003 [bug fix]
## Modified to meet new developer specification 7/15/2004 KQ
## Modified for new Scythe and rngs 7/26/2004 KQ
## Modified to handle marginal likelihood calculation 1/27/2006 KQ
+##
+## 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) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+## Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+## and Jong Hee Park
+##########################################################################
+
"MCMCpoisson" <-
- function(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
+ function(formula, data=NULL, burnin = 1000, mcmc = 10000,
thin=1, tune=1.1, verbose = 0, seed = NA, beta.start = NA,
b0 = 0, B0 = 0,
marginal.likelihood = c("none", "Laplace"),...) {
@@ -27,7 +38,7 @@
lecuyer.stream <- seeds[[3]]
## form response and model matrices
- holder <- parse.formula(formula, data)
+ holder <- parse.formula(formula, data=data)
Y <- holder[[1]]
X <- holder[[2]]
xnames <- holder[[3]]
diff --git a/R/MCMCpoissonChangepoint.R b/R/MCMCpoissonChangepoint.R
index 348ee34..8c21797 100644
--- a/R/MCMCpoissonChangepoint.R
+++ b/R/MCMCpoissonChangepoint.R
@@ -1,3 +1,4 @@
+##########################################################################
## sample from the posterior distribution
## of a Poisson model with multiple changepoints
## using linked C++ code in Scythe 1.0
@@ -5,6 +6,15 @@
## JHP 07/01/2007
##
## Revised on 09/12/2007 JHP
+##
+## 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) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+## Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+## and Jong Hee Park
+##########################################################################
"MCMCpoissonChangepoint"<-
function(data, m = 1, c0 = NA, d0 = NA, a = NULL, b = NULL,
diff --git a/R/MCMCprobit.R b/R/MCMCprobit.R
index a2edfd8..edb261c 100644
--- a/R/MCMCprobit.R
+++ b/R/MCMCprobit.R
@@ -1,3 +1,4 @@
+##########################################################################
## sample from the posterior distribution of a probit
## model in R using linked C++ code in Scythe
##
@@ -5,10 +6,19 @@
## Modified to meet new developer specification 7/26/2004 KQ
## Modified for new Scythe and rngs 7/26/2004 KQ
## Modified to handle marginal likelihood calculation 1/27/2006 KQ
+##
+## 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) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+## Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+## and Jong Hee Park
+##########################################################################
"MCMCprobit" <-
- function(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
+ function(formula, data=NULL, burnin = 1000, mcmc = 10000,
thin = 1, verbose = 0, seed = NA, beta.start = NA,
b0 = 0, B0 = 0, bayes.resid=FALSE,
marginal.likelihood = c("none", "Laplace"), ...) {
@@ -27,7 +37,7 @@
lecuyer.stream <- seeds[[3]]
## form response and model matrices
- holder <- parse.formula(formula, data)
+ holder <- parse.formula(formula, data=data)
Y <- holder[[1]]
X <- holder[[2]]
xnames <- holder[[3]]
diff --git a/R/MCMCregress.R b/R/MCMCregress.R
index 68bf281..e6537da 100644
--- a/R/MCMCregress.R
+++ b/R/MCMCregress.R
@@ -1,3 +1,4 @@
+##########################################################################
## MCMCregress.R samples from the posterior distribution of a Gaussian
## linear regression model in R using linked C++ code in Scythe
##
@@ -6,9 +7,18 @@
## Modified to meet new developer specification 6/18/2004 KQ
## Modified for new Scythe and rngs 7/22/2004 ADM
## Modified to handle marginal likelihood calculation 1/26/2006 KQ
+##
+## 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) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+## Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+## and Jong Hee Park
+##########################################################################
"MCMCregress" <-
- function(formula, data=parent.frame(), burnin = 1000, mcmc = 10000,
+ function(formula, data=NULL, burnin = 1000, mcmc = 10000,
thin=1, verbose = 0, seed = NA, beta.start = NA,
b0 = 0, B0 = 0, c0 = 0.001, d0 = 0.001,
marginal.likelihood = c("none", "Laplace", "Chib95"),
@@ -29,7 +39,7 @@
lecuyer.stream <- seeds[[3]]
## form response and model matrices
- holder <- parse.formula(formula, data)
+ holder <- parse.formula(formula, data=data)
Y <- holder[[1]]
X <- holder[[2]]
xnames <- holder[[3]]
diff --git a/R/MCMCtobit.R b/R/MCMCtobit.R
index e5735fd..92bf528 100644
--- a/R/MCMCtobit.R
+++ b/R/MCMCtobit.R
@@ -1,5 +1,19 @@
+##########################################################################
+## tobit regression model
+##
+## 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) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+## Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+## and Jong Hee Park
+##########################################################################
+
+
+
"MCMCtobit" <-
-function(formula, data=parent.frame(), below = 0.0, above = Inf,
+function(formula, data=NULL, below = 0.0, above = Inf,
burnin = 1000, mcmc = 10000,
thin=1, verbose = 0, seed = NA, beta.start = NA,
b0 = 0, B0 = 0, c0 = 0.001, d0 = 0.001, ...) {
@@ -29,7 +43,7 @@ function(formula, data=parent.frame(), below = 0.0, above = Inf,
lecuyer.stream <- seeds[[3]]
# form response and model matrices
- holder <- parse.formula(formula, data)
+ holder <- parse.formula(formula, data=data)
Y <- holder[[1]]
X <- holder[[2]]
xnames <- holder[[3]]
diff --git a/R/MCmodels.R b/R/MCmodels.R
index fca6c63..420e361 100644
--- a/R/MCmodels.R
+++ b/R/MCmodels.R
@@ -1,6 +1,20 @@
+##########################################################################
+## simple instructional models using Monte Carlo simulation
+##
+## 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) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+## Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+## and Jong Hee Park
+##########################################################################
+
## Monte Carlo simulation from the likelihood of a
## binomial distribution with a Beta(alpha, beta) prior
## ADM 1/25/2006
+
+
MCbinomialbeta <- function(y, n, alpha=1, beta=1, mc=1000, ...) {
# check data
diff --git a/R/automate.R b/R/automate.R
index d2f397f..ff4cd1f 100644
--- a/R/automate.R
+++ b/R/automate.R
@@ -1,3 +1,4 @@
+##########################################################################
## this function automates the Scythe C++ call making book-keeping
## much easier
##
@@ -42,7 +43,17 @@
## This also build a skeleton C++ template and clean R template
## for MCMCpack if developer=TRUE.
##
+## 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.
+##
## Updated by ADM and KQ 1/25/2006 (to allow for multiple nonconsts)
+##
+## Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+## Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+## and Jong Hee Park
+##########################################################################
+
"auto.Scythe.call" <-
function(output.object, cc.fun.name, package="MCMCpack",
diff --git a/R/btsutil.R b/R/btsutil.R
index deb1722..cd5f857 100644
--- a/R/btsutil.R
+++ b/R/btsutil.R
@@ -1,15 +1,27 @@
-########## Utility Functions for Bayesian Times Series Models ##########
-
+##########################################################################
+## Utility Functions for Bayesian Times Series Models
+##
## written and maintained by:
## Jong Hee Park
## Department of Political Science
## University of Chicago
## jhp at uchicago.edu
-
+##
## Revised on 09/12/2007 JHP
-
+##
## NOTE: only the plot functions are documented and exported in the
## NAMESPACE.
+##
+## 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) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+## Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+## and Jong Hee Park
+##########################################################################
+
+
##############################################################
## Helper functions for MCMCPoissonChangepoint()
diff --git a/R/distn.R b/R/distn.R
index 1d7b592..00ff600 100644
--- a/R/distn.R
+++ b/R/distn.R
@@ -1,11 +1,16 @@
-########## 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).
+##########################################################################
+## Density Functions and Random Number Generators
+##
+## 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) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+## Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+## and Jong Hee Park
+##########################################################################
+
+
##
## Wishart
diff --git a/R/hidden.R b/R/hidden.R
index a2ce184..118d134 100644
--- a/R/hidden.R
+++ b/R/hidden.R
@@ -205,7 +205,7 @@
# or for another model by passing default starting values to
# the function
"coef.start" <-
- function(beta.start, K, formula, family, data, defaults=NA) {
+ function(beta.start, K, formula, family, data=NULL, defaults=NA) {
if (is.na(beta.start)[1] & is.na(defaults)[1]){ # use GLM estimates
beta.start <- matrix(coef(glm(formula, family=family, data=data)), K, 1)
@@ -663,7 +663,7 @@
stop("Please respecify and call ", calling.function(), " again.",
call.=FALSE)
}
- if( all(lec.seed[1:3]) == 0 ){
+ if( all(lec.seed[1:3] == 0 )){
cat("Error: first three L'Ecuyer seeds == 0\n")
stop("Please respecify and call ", calling.function(), " again.",
call.=FALSE)
@@ -674,7 +674,7 @@
stop("Please respecify and call ", calling.function(), " again.",
call.=FALSE)
}
- if( all(lec.seed[4:6]) == 0 ){
+ if( all(lec.seed[4:6] == 0 )){
cat("Error: last three L'Ecuyer seeds == 0\n")
stop("Please respecify and call ", calling.function(), " again.",
call.=FALSE)
@@ -719,16 +719,15 @@
# parse formula and return a list that contains the model response
# matrix as element one, and the model matrix as element two
"parse.formula" <-
- function(formula, data, intercept=TRUE, justX=FALSE) {
+ function(formula, data=NULL, intercept=TRUE, justX=FALSE) {
# 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$justX <- NULL
mf$drop.unused.levels <- TRUE
mf[[1]] <- as.name("model.frame")
mf <- eval(mf, sys.frame(sys.parent()))
+ mt <- attr(mf, "terms")
if (!intercept){
attributes(mt)$intercept <- 0
}
diff --git a/R/procrust.R b/R/procrust.R
index d6bf7d8..3e24c5e 100644
--- a/R/procrust.R
+++ b/R/procrust.R
@@ -1,3 +1,4 @@
+##########################################################################
## function that performs procrustes transformation on X with target Xstar
##
## returns the rotation matrix R, translation vector tt,
@@ -10,10 +11,19 @@
## Based on Borg and Groenen 1997. Modern Multidimensional
## Scaling. New York: Springer. pp. 340-342.
##
+## 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.
+##
## Kevin Quinn
## Harvard University
## 6/13/2005
##
+## Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+## Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+## and Jong Hee Park
+##########################################################################
+
procrustes <- function(X, Xstar, translation=FALSE, dilation=FALSE){
if (nrow(X) != nrow(Xstar)){
cat("X and Xstar do not have same number of rows.\n")
diff --git a/R/scythe.R b/R/scythe.R
index 574989d..0a6d5a3 100644
--- a/R/scythe.R
+++ b/R/scythe.R
@@ -1,4 +1,14 @@
-########## Scythe Inter-Operation Functions ##########
+##########################################################################
+## Scythe Inter-Operation Functions
+##
+## 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) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+## Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+## and Jong Hee Park
+##########################################################################
# writes a matrix out to an ASCII file that can be read by Scythe.
diff --git a/R/tomog.R b/R/tomog.R
index cb4a7d0..f7ee433 100644
--- a/R/tomog.R
+++ b/R/tomog.R
@@ -1,11 +1,22 @@
-########## Tomography Plots for Ecological Inference ##########
-
+##########################################################################
+## Tomography Plots for Ecological Inference
+##
## produces tomography plots (see King, 1997, A Solution to the
## Ecological Inference Problem, Princeton University Press)
##
+## 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.
+##
## KQ 11/9/2002
##
## Modification added suggested by David Hugh-Jones 6/10/2006
+##
+## Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+## Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+## and Jong Hee Park
+##########################################################################
+
"tomogplot" <-
function(r0, r1, c0, c1,
diff --git a/R/utility.R b/R/utility.R
index b2515c6..fe5c1e0 100644
--- a/R/utility.R
+++ b/R/utility.R
@@ -1,4 +1,15 @@
-########## Utility Functions ##########
+##########################################################################
+## Utility Functions
+##
+## 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) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+## Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+## and Jong Hee Park
+##########################################################################
+
# takes a symmetric matrix x and returns lower diagonal
# note: does not check for symmetry
diff --git a/R/zzz.R b/R/zzz.R
index f078b95..8af0a9b 100644
--- a/R/zzz.R
+++ b/R/zzz.R
@@ -1,3 +1,15 @@
+##########################################################################
+## start-up and clean-up functions
+##
+## 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) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+## Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+## and Jong Hee Park
+##########################################################################
+
.onAttach <- function(...) {
# figure out year automatically (probably could be done more elegantly)
diff --git a/README b/README
index d5461ee..f4de977 100644
--- a/README
+++ b/README
@@ -42,13 +42,12 @@ you have any problems or questions.
--
Andrew D. Martin, Ph.D.
+Professor and Chair, Political Science, Arts & Sciences
Professor and CERL Director, School of Law
-Professor, Political Science, Arts & Sciences
Washington University in St. Louis
(314) 935-5863 (Office)
-(314) 935-5150 (Fax)
+(314) 935-5856 (Fax)
-Office: Anheuser-Busch 470
Email: admartin at wustl.edu
WWW: http://adm.wustl.edu
diff --git a/man/MCMCSVDreg.Rd b/man/MCMCSVDreg.Rd
index b181bc1..5b4e576 100644
--- a/man/MCMCSVDreg.Rd
+++ b/man/MCMCSVDreg.Rd
@@ -13,7 +13,7 @@
}
\usage{
-MCMCSVDreg(formula, data=parent.frame(), burnin = 1000, mcmc = 10000,
+MCMCSVDreg(formula, data=NULL, burnin = 1000, mcmc = 10000,
thin=1, verbose = 0, seed = NA, tau2.start = 1,
g0 = 0, a0 = 0.001, b0 = 0.001, c0=2, d0=2, w0=1,
beta.samp=FALSE, intercept=TRUE, ...)}
diff --git a/man/MCMCfactanal.Rd b/man/MCMCfactanal.Rd
index 5a20c1d..6b41536 100644
--- a/man/MCMCfactanal.Rd
+++ b/man/MCMCfactanal.Rd
@@ -13,7 +13,7 @@
\usage{
MCMCfactanal(x, factors, lambda.constraints=list(),
- data=parent.frame(), burnin = 1000, mcmc = 20000,
+ data=NULL, burnin = 1000, mcmc = 20000,
thin=1, verbose = 0, seed = NA,
lambda.start = NA, psi.start = NA,
l0=0, L0=0, a0=0.001, b0=0.001,
diff --git a/man/MCMCirt1d.Rd b/man/MCMCirt1d.Rd
index 96d94ae..2aa0f76 100644
--- a/man/MCMCirt1d.Rd
+++ b/man/MCMCirt1d.Rd
@@ -4,7 +4,7 @@
Model}
\description{
This function generates a sample from the posterior distribution of a one
- dimentional item response theory (IRT) model, with Normal
+ dimensional item response theory (IRT) model, with Normal
priors on the subject abilities (ideal points), and
multivariate Normal priors on the item parameters. The user
supplies data and priors, and a sample from the posterior
diff --git a/man/MCMCirtHier1d.Rd b/man/MCMCirtHier1d.Rd
new file mode 100644
index 0000000..048b30f
--- /dev/null
+++ b/man/MCMCirtHier1d.Rd
@@ -0,0 +1,305 @@
+\name{MCMCirtHier1d}
+\alias{MCMCirtHier1d}
+\title{Markov Chain Monte Carlo for Hierarchical One Dimensional Item Response Theory
+ Model, Covariates Predicting Latent Ideal Point (Ability)}
+\description{
+ This function generates a sample from the posterior distribution of a one
+ dimensional item response theory (IRT) model, with
+ multivariate Normal priors on the item parameters, and a
+ Normal-Inverse Gamma hierarchical prior on subject ideal points (abilities). The user
+ supplies item-response data, subject covariates, and priors. Note that
+ this identification strategy obviates the constraints used on
+ theta in \code{\link[MCMCpack]{MCMCirt1d}}. A sample from the posterior
+ distribution is returned as an mcmc object, which can be
+ subsequently analyzed with functions provided in the coda
+ package.
+
+ If you are interested in fitting K-dimensional item response theory
+ models, or would rather identify the model by placing constraints
+ on the item parameters, please see \code{\link[MCMCpack]{MCMCirtKd}}.
+ }
+
+\usage{
+MCMCirtHier1d(datamatrix, Xjdata,
+ burnin = 1000, mcmc = 20000, thin=1,
+ verbose = 0, seed = NA,
+ theta.start = NA, a.start = NA, b.start = NA,
+ beta.start=NA, b0=0, B0=.01, c0=.001, d0=.001,
+ ab0=0, AB0=.25, store.item = FALSE, store.ability=TRUE,
+ drop.constant.items=TRUE,
+ marginal.likelihood=c("none","Chib95"),
+ px=TRUE,px_a0 = 10, px_b0=10,
+ ... ) }
+
+\arguments{
+ \item{datamatrix}{The matrix of data. Must be 0, 1, or missing values.
+ The rows of \code{datamatrix} correspond to subjects and the
+ columns correspond to items.}
+
+ \item{Xjdata}{A \code{data.frame} containing second-level predictor covariates for
+ ideal points \eqn{\theta}{theta}. Predictors are modeled as a
+ linear regression on the mean vector of
+ \eqn{\theta}{theta}; the posterior sample contains regression
+ coefficients \eqn{\beta}{beta} and common variance
+ \eqn{\sigma^2}{sigma^2}. See Rivers (2003) for a
+ thorough discussion of identification of IRT models.}
+
+ \item{burnin}{The number of burn-in iterations for the sampler.}
+
+ \item{mcmc}{The number of Gibbs iterations for the sampler.}
+
+ \item{thin}{The thinning interval used in the simulation. The number of
+ Gibbs iterations must be divisible by this value.}
+
+ \item{verbose}{A switch which determines whether or not the progress of
+ the sampler is printed to the screen. If \code{verbose} is greater
+ than 0 then every \code{verbose}th iteration will be printed to the
+ screen.}
+
+ \item{seed}{The seed for the random number generator. If NA, the Mersenne
+ Twister generator is used with default seed 12345; if an integer is
+ passed it is used to seed the Mersenne twister. The user can also
+ pass a list of length two to use the L'Ecuyer random number generator,
+ which is suitable for parallel computation. The first element of the
+ list is the L'Ecuyer seed, which is a vector of length six or NA (if NA
+ a default seed of \code{rep(12345,6)} is used). The second element of
+ list is a positive substream number. See the MCMCpack
+ specification for more details.}
+
+ \item{theta.start}{The starting values for the subject
+ abilities (ideal points). This can either be a scalar or a
+ column vector with dimension equal to the number of voters.
+ If this takes a scalar value, then that value will serve as
+ the starting value for all of the thetas. The default value
+ of NA will choose the starting values based on an
+ eigenvalue-eigenvector decomposition of the agreement score matrix
+ formed from the \code{datamatrix}.}
+
+ \item{a.start}{The starting values for the
+ \eqn{a}{a} difficulty parameters. This can either be
+ a scalar or a column vector with dimension equal to the
+ number of items. If this takes a scalar value, then that
+ value will serve as the starting value for all \eqn{a}{a}.
+ The default value of NA will set the starting values based on
+ a series of probit regressions that condition on the starting
+ values of theta.}
+
+ \item{b.start}{The starting values for the
+ \eqn{b}{b} discrimination parameters. This can either
+ be a scalar or a column vector with dimension equal to the
+ number of items. If this takes a scalar value, then that
+ value will serve as the starting value for all \eqn{b}{b}.
+ The default value of NA will set the starting values based on a
+ series of probit regressions that condition on the starting values
+ of theta.}
+
+
+ \item{beta.start}{The starting values for the \eqn{\beta}{beta}
+ regression coefficients that predict the means of ideal points
+ \eqn{\theta}{theta}. This can either
+ be a scalar or a column vector with length equal to the
+ number of covariates. 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 set the starting values based on a
+ linear regression of the covariates on (either provided or
+ generated) \code{theta.start}.
+ }
+
+
+ \item{b0}{The prior mean of \eqn{\beta}{beta}. Can be either a scalar or a
+ vector of length equal to the number of subject covariates. If a
+ scalar all means with be set to the passed value.}
+
+\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. A default proper but diffuse value
+ of .01 ensures finite marginal likelihood for model comparison.
+ A 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 \eqn{\theta}{theta}). 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
+ \eqn{\theta}{theta}). 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{ab0}{The prior mean of \code{(a, b)}. Can be either a
+ scalar or a 2-vector. If a scalar both means will be set to the
+ passed value. The prior mean is assumed to be the same across all
+ items.}
+
+ \item{AB0}{The prior precision of \code{(a, b)}.This can
+ either be ascalar or a 2 by 2 matrix. If this takes a scalar
+ value, then that value times an identity matrix serves as the
+ prior precision. The prior precision is assumed to be the same
+ across all items.}
+
+ \item{store.item}{A switch that determines whether or not to
+ store the item parameters for posterior analysis.
+ \emph{NOTE: In situations with many items storing the item
+ parameters takes an enormous amount of memory, so
+ \code{store.item} should only be \code{TRUE} if the chain is thinned
+ heavily, or for applications with a small number of items}.
+ By default, the item parameters are not stored.}
+
+ \item{store.ability}{A switch that determines whether or not to
+ store the ability parameters for posterior analysis.
+ \emph{NOTE: In situations with many individuals storing the ability
+ parameters takes an enormous amount of memory, so
+ \code{store.ability} should only be \code{TRUE} if the chain is thinned
+ heavily, or for applications with a small number of individuals}.
+ By default, ability parameters are stored.}
+
+ \item{drop.constant.items}{A switch that determines whether or not
+ items that have no variation
+ should be deleted before fitting the model. Default = TRUE.}
+
+ \item{marginal.likelihood}{Should the marginal likelihood of the
+ second-level model on ideal points be calculated using the method of
+ Chib (1995)? It is stored as an attribute of the posterior \code{mcmc}
+ object and suitable for comparison using \code{\link[MCMCpack]{BayesFactor}}.}
+
+ \item{px}{Use Parameter Expansion to reduce autocorrelation in the chain?
+ PX introduces an unidentified parameter \eqn{alpha} for the residual
+ variance in the latent data (Liu and Wu 1999). Default = TRUE
+}
+ \item{px_a0}{Prior shape parameter for the inverse-gamma distribution on \eqn{alpha}, the residual variance of the latent data. Default=10.}
+
+\item{px_b0}{ Prior scale parameter for the inverse-gamma distribution on \eqn{alpha}, the residual variance of the latent data. Default = 10 }
+
+ \item{...}{further arguments to be passed}
+}
+
+\value{
+ An \code{mcmc} object that contains the sample from the posterior
+ distribution. This object can be summarized by functions provided by
+ the coda package. If \code{marginal.likelihood = "Chib95"} the object will have attribute \code{logmarglike}.
+}
+
+\details{
+ \code{MCMCirtHier1d} simulates from the posterior distribution using
+ standard Gibbs sampling using data augmentation (a Normal draw
+ for the subject abilities, a multivariate Normal draw for
+ (second-level) subject
+ ability predictors, an Inverse-Gamma draw for the (second-level)
+ variance of subject abilities, a multivariate Normal
+ draw for the item parameters, and a truncated Normal draw for
+ the latent utilities). The simulation proper is done in
+ compiled C++ code to maximize efficiency. Please consult the
+ coda documentation for a comprehensive list of functions that
+ can be used to analyze the posterior sample.
+
+ The model takes the following form. We assume that each
+ subject has an subject ability (ideal point) denoted
+ \eqn{\theta_j}{theta_j} and that each item has a difficulty
+ parameter \eqn{a_i}{a_i} and discrimination parameter
+ \eqn{b_i}{b_i}. The observed choice by subject
+ \eqn{j}{j} on item \eqn{i}{i} is the observed data matrix which
+ is \eqn{(I \times J)}{(I * J)}. We assume that the choice is
+ dictated by an unobserved utility: \deqn{z_{i,j} = -\alpha_i +
+ \beta_i \theta_j + \varepsilon_{i,j}}{z_ij = -a_i +
+ b_i*theta_j + epsilon_ij} Where the errors are assumed to be
+distributed standard Normal. This constitutes the measurement or level-1
+ model. The subject abilities (ideal points) are modeled by a second
+ level Normal linear predictor for subject covariates \code{Xjdata}, with
+ common variance \eqn{\sigma^2}{sigma^2}.
+The parameters of interest are
+ the subject abilities (ideal points), item parameters, and
+ second-level coefficients.
+
+ We assume the following priors. For the subject abilities (ideal points):
+ \deqn{\theta_j \sim \mathcal{N}(\mu_{\theta} ,T_{0}^{-1})}{theta_j ~ N(t0, T0^{-1})}
+ For the item parameters, the prior is:
+ \deqn{\left[a_i, b_i \right]' \sim \mathcal{N}_2
+ (ab_{0},AB_{0}^{-1})}{[alpha_i beta_i]' ~ N_2 (ab0, AB0^{-1})}
+
+ 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.
+}
+ \author{Michael Malecki, \email{malecki at wustl.edu}, \url{http://malecki.wustl.edu}.}
+
+\references{
+ James H. Albert. 1992. ``Bayesian Estimation of Normal Ogive Item Response
+ Curves Using Gibbs Sampling." \emph{Journal of Educational Statistics}.
+ 17: 251--269.
+
+ Joshua Clinton, Simon Jackman, and Douglas Rivers. 2004. ``The Statistical
+ Analysis of Roll Call Data." \emph{American Political Science Review}
+ 98: 355--370.
+
+ Valen E. Johnson and James H. Albert. 1999. ``Ordinal Data Modeling."
+ Springer: New York.
+
+ Liu, Jun S. and Ying Nian Wu. 1999. ``Parameter Expansion for Data Augmentation.'' \emph{Journal of the American Statistical Association} 94: 1264--1274.
+
+ Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007.
+ \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/}.
+
+ Douglas Rivers. 2004. ``Identification of Multidimensional Item-Response
+ Models." Stanford University, typescript.
+}
+
+
+\examples{
+ \dontrun{
+data(SupremeCourt)
+
+Xjdata <- data.frame(presparty= c(1,1,0,1,1,1,1,0,0),
+ sex= c(0,0,1,0,0,0,0,1,0))
+
+## Parameter Expansion reduces autocorrelation.
+ posterior1 <- MCMCirtHier1d(t(SupremeCourt),
+ burnin=50000, mcmc=10000, thin=20,
+ verbose=10000,
+ Xjdata=Xjdata,
+ marginal.likelihood="Chib95",
+ px=TRUE)
+
+## But, you can always turn it off.
+ posterior2 <- MCMCirtHier1d(t(SupremeCourt),
+ burnin=50000, mcmc=10000, thin=20,
+ verbose=10000,
+ Xjdata=Xjdata,
+ #marginal.likelihood="Chib95",
+ px=FALSE)
+## Note that the hierarchical model has greater autocorrelation than
+## the naive IRT model.
+ posterior0 <- MCMCirt1d(t(SupremeCourt),
+ theta.constraints=list(Scalia="+", Ginsburg="-"),
+ B0.alpha=.2, B0.beta=.2,
+ burnin=50000, mcmc=100000, thin=100, verbose=10000,
+ store.item=FALSE)
+
+## Randomly 10% Missing -- this affects the expansion parameter, increasing
+## the variance of the (unidentified) latent parameter alpha.
+
+scMiss <- SupremeCourt
+scMiss[matrix(as.logical(rbinom( nrow(SupremeCourt)*ncol(SupremeCourt), 1, .1)), dim(SupremeCourt))] <- NA
+
+ posterior1.miss <- MCMCirtHier1d(t(scMiss),
+ burnin=80000, mcmc=10000, thin=20,
+ verbose=10000,
+ Xjdata=Xjdata,
+ marginal.likelihood="Chib95",
+ px=TRUE)
+
+ }
+}
+
+\keyword{models}
+
+\seealso{\code{\link[coda]{plot.mcmc}},\code{\link[coda]{summary.mcmc}},
+\code{\link[MCMCpack]{MCMCirtKd}}}
+
diff --git a/man/MCMClogit.Rd b/man/MCMClogit.Rd
index 81a21b7..64d9d36 100644
--- a/man/MCMClogit.Rd
+++ b/man/MCMClogit.Rd
@@ -11,7 +11,7 @@
}
\usage{
-MCMClogit(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
+MCMClogit(formula, data=NULL, burnin = 1000, mcmc = 10000,
thin=1, tune=1.1, verbose = 0, seed = NA, beta.start = NA,
b0 = 0, B0 = 0, user.prior.density=NULL, logfun=TRUE,
marginal.likelihood=c("none", "Laplace"), ...) }
diff --git a/man/MCMCmnl.Rd b/man/MCMCmnl.Rd
index 13eaee7..9a8eec8 100644
--- a/man/MCMCmnl.Rd
+++ b/man/MCMCmnl.Rd
@@ -12,7 +12,7 @@
}
\usage{
-MCMCmnl(formula, baseline = NULL, data = parent.frame(),
+MCMCmnl(formula, baseline=NULL, data=NULL,
burnin = 1000, mcmc = 10000, thin = 1,
mcmc.method = c("IndMH", "RWM", "slice"), tune = 1, tdf=6,
verbose = 0, seed = NA, beta.start = NA, b0 = 0, B0 = 0, ...) }
@@ -175,7 +175,7 @@ that can be used to analyze the posterior sample.
choicevar(distPvdA, "sqdist", "PvdA") +
choicevar(distVVD, "sqdist", "VVD") +
choicevar(distCDA, "sqdist", "CDA"),
- baseline="D66", mcmc.method="MH", B0=0,
+ baseline="D66", mcmc.method="IndMH", B0=0,
verbose=500, mcmc=100000, thin=10, tune=1.0,
data=Nethvote)
@@ -187,7 +187,7 @@ that can be used to analyze the posterior sample.
## just individual-specific X vars
post2<- MCMCmnl(vote ~
relig + class + income + educ + age + urban,
- baseline="D66", mcmc.method="MH", B0=0,
+ baseline="D66", mcmc.method="IndMH", B0=0,
verbose=500, mcmc=100000, thin=10, tune=0.5,
data=Nethvote)
@@ -203,7 +203,7 @@ that can be used to analyze the posterior sample.
choicevar(distVVD, "sqdist", "VVD") +
choicevar(distCDA, "sqdist", "CDA") +
relig + class + income + educ + age + urban,
- baseline="D66", mcmc.method="MH", B0=0,
+ baseline="D66", mcmc.method="IndMH", B0=0,
verbose=500, mcmc=100000, thin=10, tune=0.5,
data=Nethvote)
@@ -219,7 +219,7 @@ that can be used to analyze the posterior sample.
choicevar(distVVD, "sqdist", "4") +
choicevar(distCDA, "sqdist", "1") +
relig + class + income + educ + age + urban,
- baseline="2", mcmc.method="MH", B0=0,
+ baseline="2", mcmc.method="IndMH", B0=0,
verbose=500, mcmc=100000, thin=10, tune=0.5,
data=Nethvote)
diff --git a/man/MCMCpoisson.Rd b/man/MCMCpoisson.Rd
index d41c4ed..d51ace3 100644
--- a/man/MCMCpoisson.Rd
+++ b/man/MCMCpoisson.Rd
@@ -11,7 +11,7 @@
}
\usage{
-MCMCpoisson(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
+MCMCpoisson(formula, data = NULL, burnin = 1000, mcmc = 10000,
thin = 1, tune = 1.1, verbose = 0, seed = NA, beta.start = NA,
b0 = 0, B0 = 0, marginal.likelihood = c("none", "Laplace"), ...) }
diff --git a/man/MCMCprobit.Rd b/man/MCMCprobit.Rd
index 04552ca..52287b5 100644
--- a/man/MCMCprobit.Rd
+++ b/man/MCMCprobit.Rd
@@ -11,7 +11,7 @@
}
\usage{
-MCMCprobit(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
+MCMCprobit(formula, data = NULL, burnin = 1000, mcmc = 10000,
thin = 1, verbose = 0, seed = NA, beta.start = NA,
b0 = 0, B0 = 0, bayes.resid = FALSE,
marginal.likelihood=c("none", "Laplace"), ...) }
diff --git a/man/MCMCregress.Rd b/man/MCMCregress.Rd
index fb694f2..0bce653 100644
--- a/man/MCMCregress.Rd
+++ b/man/MCMCregress.Rd
@@ -13,7 +13,7 @@
}
\usage{
-MCMCregress(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
+MCMCregress(formula, data = NULL, burnin = 1000, mcmc = 10000,
thin = 1, verbose = 0, seed = NA, beta.start = NA,
b0 = 0, B0 = 0, c0 = 0.001, d0 = 0.001,
marginal.likelihood = c("none", "Laplace", "Chib95"), ...) }
diff --git a/man/MCMCtobit.Rd b/man/MCMCtobit.Rd
index 0db4277..9494768 100644
--- a/man/MCMCtobit.Rd
+++ b/man/MCMCtobit.Rd
@@ -13,7 +13,7 @@
analyzed with functions provided in the coda package.
}
\usage{
-MCMCtobit(formula, data = parent.frame(), below = 0, above = Inf,
+MCMCtobit(formula, data = NULL, below = 0, above = Inf,
burnin = 1000, mcmc = 10000, thin = 1, verbose = 0, seed = NA,
beta.start = NA, b0 = 0, B0 = 0, c0 = 0.001, d0 = 0.001, ...)
}
diff --git a/man/PErisk.Rd b/man/PErisk.Rd
index 7027b91..d032452 100644
--- a/man/PErisk.Rd
+++ b/man/PErisk.Rd
@@ -10,7 +10,7 @@ Political Economic Risk Data from 62 Countries in 1987.
A data frame with 62 observations on the following 9 variables. All
data points are from 1987. See Quinn (2004) for more details.
\describe{
- \item{country}{a factor with levels \code{Argentina} \code{Australia} \code{Austria} \code{Bangladesh} \code{Belgium} \code{Bolivia} \code{Botswana} \code{Brazil} \code{Burma} \code{Cameroon} \code{Canada} \code{Chile} \code{Colombia} \code{Congo-Kinshasa} \code{Costa Rica} \code{Cote d'Ivoire} \code{Denmark} \code{Dominican Republic} \code{Ecuador} \code{Finland} \code{Gambia, The} \code{Ghana} \code{Greece} \code{Hungary} \code{India} \code{Indonesia} \code{Iran} \code{Ireland} \co [...]
+ \item{country}{a factor with levels \code{Argentina} through \code{Zimbabwe}}
\item{courts}{an ordered factor with levels \code{0} <
\code{1}.\code{courts} is an indicator of whether the country in
question is judged to have an independent judiciary. From Henisz
diff --git a/src/MCMCSVDreg.cc b/src/MCMCSVDreg.cc
index dec6279..36ad39f 100644
--- a/src/MCMCSVDreg.cc
+++ b/src/MCMCSVDreg.cc
@@ -1,3 +1,4 @@
+//////////////////////////////////////////////////////////////////////////
// MCMCSVDreg.R samples from the posterior distribution of a Gaussian
// linear regression model in which the X matrix has been decomposed
// with an SVD. Useful for prediction when number of columns of X
@@ -20,10 +21,13 @@
// PUBLIC LICENSE Version 2, June 1991. See the package LICENSE
// file for more information.
//
-// Copyright (C) 2004 Andrew D. Martin and Kevin M. Quinn
-//
// KQ 9/9/2005
// KQ 8/1/2007 ported to Scythe 1.0.2
+//
+// Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+// Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+// and Jong Hee Park
+//////////////////////////////////////////////////////////////////////////
#ifndef MCMCSVDREG_CC
diff --git a/src/MCMCdynamicEI.cc b/src/MCMCdynamicEI.cc
index 5d604e1..58d2450 100644
--- a/src/MCMCdynamicEI.cc
+++ b/src/MCMCdynamicEI.cc
@@ -1,3 +1,4 @@
+//////////////////////////////////////////////////////////////////////////
// fits a model derived from Wakefield's baseline model for
// ecological inference in which logit(p_i) follows a random walk in time
// a priori. The model is fit using Wakefield's normal approximation
@@ -11,6 +12,13 @@
// KQ 7/20/2004 [minor changes regarding output and user interrupts]
// ADM 7/24/2004 [updated to new Scythe version]
// KQ 7/30/2007 [updated to Scythe 1.0.X]
+//
+//
+// Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+// Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+// and Jong Hee Park
+//////////////////////////////////////////////////////////////////////////
+
#ifndef MCMCDYNAMICEI_CC
diff --git a/src/MCMCfactanal.cc b/src/MCMCfactanal.cc
index 723e5ba..54ae8c3 100644
--- a/src/MCMCfactanal.cc
+++ b/src/MCMCfactanal.cc
@@ -1,3 +1,4 @@
+//////////////////////////////////////////////////////////////////////////
// MCMCfactanal.cc is C++ code to estimate a factor analysis model
//
// Andrew D. Martin
@@ -14,13 +15,16 @@
// PUBLIC LICENSE Version 2, June 1991. See the package LICENSE
// file for more information.
//
-// Copyright (C) 2004 Andrew D. Martin and Kevin M. Quinn
-//
// revised version of older MCMCfactanal 5/11/2004 KQ
// updated to new verion of scythe 7/25/2004 ADM
// updated to Scythe 1.0.X 7/10/2007 ADM
// finished update to Scythe 1.0.X 7/30/2007 KQ
//
+// Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+// Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+// and Jong Hee Park
+//////////////////////////////////////////////////////////////////////////
+
#ifndef MCMCFACTANAL_CC
#define MCMCFACTANAL_CC
diff --git a/src/MCMCfcds.h b/src/MCMCfcds.h
index e3df141..e77e9c8 100644
--- a/src/MCMCfcds.h
+++ b/src/MCMCfcds.h
@@ -1,3 +1,4 @@
+//////////////////////////////////////////////////////////////////////////
// MCMCfcds.h is the header file for MCMCfcds.cc. It contains declarations
// for a number of functions that produce draws from full conditional
// distributions.
@@ -16,10 +17,14 @@
// PUBLIC LICENSE Version 2, June 1991. See the package LICENSE
// file for more information.
//
-// Copyright (C) 2004 Andrew D. Martin and Kevin M. Quinn
-//
// KQ 6/10/2004
// DBP 7/01/2007 [ported to scythe 1.0.x (partial)]
+//
+// Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+// Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+// and Jong Hee Park
+//////////////////////////////////////////////////////////////////////////
+
#ifndef MCMCFCDS_H
@@ -402,14 +407,96 @@ NormNormfactanal_Lambda_draw(Matrix<>& Lambda,
}
+// update ability parameters (ideal points) for one dimensional
+// Hierarchical item response model.
+// 2008-11-18 now deals with PX alpha and holds on to thetahat
+template <typename RNGTYPE>
+void hirt_theta_update1 (Matrix<>& theta, Matrix<>& thetahat,
+ const Matrix<>& Z,
+ const Matrix<>& eta,
+ const Matrix<>& beta, const Matrix<>& Xj,
+ const double& sigma2,
+ const double& alpha,
+ rng<RNGTYPE>& stream)
+{
+
+ const unsigned int J = Z.rows();
+ const unsigned int K = Z.cols();
+ // Get level1 prior mean
+ const Matrix<double> Xbeta = (Xj * beta);
+
+ // a and b are backwards here, different parameterizations
+ // common for edu people and us.
+ const Matrix<> b = eta(_, 0); // location
+ const Matrix<> a = eta(_, 1); // relevance
+ // calculate the posterior variance outside the justice specific loop
+ const double sig2_inv = 1.0 / sigma2;
+ const Matrix<double> apa = crossprod(a);
+ const Matrix<double> theta_post_var = scythe::invpd(apa + sig2_inv);
+ const double theta_post_sd = std::sqrt(theta_post_var[0]);
+ // sample for each justice
+ for (unsigned int j = 0; j < J; ++j) {
+ thetahat(j) = 0.0;
+ for (unsigned int k = 0; k < K; ++k) {
+ thetahat(j) += a[k] * ( Z(j,k) + b[k]); // bill contribution
+ }
+ thetahat(j) += (Xbeta[j] / sigma2); // j prior level1 contribution
+
+ thetahat(j) *= theta_post_var[0];
+ const double t = thetahat(j) / alpha;
+ theta(j) = stream.rnorm( t , theta_post_sd );
+ }
+}
+// update item (case, roll call) parameters for item response model
+// note: works only for one-dimensional case
+// updated for PX (alpha) and hold on to etahat 2008-11-18
+template <typename RNGTYPE>
+void hirt_eta_update1 (Matrix<>& eta, Matrix<>& etahat, const Matrix<>& Z,
+ const Matrix<>& theta, const Matrix<>& AB0,
+ const Matrix<>& AB0ab0,
+ const double& alpha,
+ rng<RNGTYPE>& stream)
+{
+
+ // define constants
+ const unsigned int J = theta.rows();
+ const unsigned int K = Z.cols();
+ // perform update
+ //const Matrix<double> Ttheta_star = t(cbind(-1.0*ones<double>(J,1),theta)); // only needed for option 2
+ const Matrix<> tpt(2,2);
+ for (unsigned int i = 0; i < J; ++i) {
+ const double theta_i = theta(i);
+ tpt(0,1) -= theta_i;
+ tpt(1,1) += std::pow(theta_i, 2.0);
+ }
+ tpt(1,0) = tpt(0,1);
+ tpt(0,0) = J;
+ const Matrix<> eta_post_var = invpd(tpt + AB0);
+ const Matrix<> eta_post_C = cholesky(eta_post_var);
+ for (unsigned int k = 0; k < K; ++k) {
+ const Matrix<> TZ(2, 1);
+ for (unsigned int j = 0; j < J; ++j) {
+ TZ[0] -= Z(j,k);
+ TZ[1] += Z(j,k) * theta[j];
+ }
+ Matrix<> eta_post_mean = eta_post_var * (TZ + AB0ab0);
+ etahat(k,0) = eta_post_mean(0);
+ etahat(k,1) = eta_post_mean(1);
+ eta_post_mean /= alpha;
+ const Matrix<> new_eta = gaxpy(eta_post_C, stream.rnorm(2, 1, 0, 1),
+ (eta_post_mean) );
+ eta(k,0) = new_eta(0);
+ eta(k,1) = new_eta(1);
+ //Rprintf("\n\a: %3.1f,b:%3.1f ",eta(k,0),eta(k,1));
-
+ }
+}
diff --git a/src/MCMChierEI.cc b/src/MCMChierEI.cc
index f290920..855bc8f 100644
--- a/src/MCMChierEI.cc
+++ b/src/MCMChierEI.cc
@@ -1,3 +1,4 @@
+//////////////////////////////////////////////////////////////////////////
// fits Wakefield's hierarchical model for ecological inference using
// Wakefield's normal approximation to the binomial convolution likelihood
// and slice sampling and Gibbs sampling to sample from the posterior
@@ -9,6 +10,11 @@
// KQ 8/10/2004 bug fix and major overhaul of sampling scheme
// ADM 7/??/2007 ported to Scythe 1.0.X
// KQ 7/29/2007 some minor fixes (still not fully tested)
+//
+// Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+// Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+// and Jong Hee Park
+//////////////////////////////////////////////////////////////////////////
#ifndef MCMCHIEREI_CC
diff --git a/src/MCMCirt1d.cc b/src/MCMCirt1d.cc
index c7889ac..74e7487 100644
--- a/src/MCMCirt1d.cc
+++ b/src/MCMCirt1d.cc
@@ -1,3 +1,4 @@
+//////////////////////////////////////////////////////////////////////////
// MCMCirt1d.cc is C++ code to estimate a one-dimensional item response
// theory model.
//
@@ -6,6 +7,11 @@
// completely rewritten and optimized for the 1-d case 8/2/2004 KQ
// storage changed to save memory KQ 1/27/2006
// DBP 7/3/07 [ported to scythe 1.0.x]
+//
+// Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+// Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+// and Jong Hee Park
+//////////////////////////////////////////////////////////////////////////
#ifndef MCMCIRT1D_CC
#define MCMCIRT1D_CC
diff --git a/src/MCMCirtHier1d.cc b/src/MCMCirtHier1d.cc
new file mode 100644
index 0000000..bca391f
--- /dev/null
+++ b/src/MCMCirtHier1d.cc
@@ -0,0 +1,350 @@
+// MCMCirtHier1d.cc is C++ code to estimate a one-dimensional item response
+// theory model with subject-level predictors beta and common variance
+// sigma2.
+//
+// ADM and KQ 1/15/2003
+// ADM 7/28/2004 [updated to new Scythe version]
+// completely rewritten and optimized for the 1-d case 8/2/2004 KQ
+// storage changed to save memory KQ 1/27/2006
+// DBP 7/3/07 [ported to scythe 1.0.x]
+//
+// MJM added second level and marginal likelihood thereof
+// MJM implemented parameter expansion (alpha) for latent variance 2008-11-18
+
+#ifndef MCMCIRTHIER1D_CC
+#define MCMCIRTHIER1D_CC
+
+#include "MCMCrng.h"
+#include "matrix.h"
+#include "distributions.h"
+#include "stat.h"
+#include "la.h"
+#include "ide.h"
+#include "smath.h"
+#include "rng.h"
+#include "mersenne.h"
+#include "lecuyer.h"
+
+#include <R.h> // needed to use Rprintf()
+#include <R_ext/Utils.h> // needed to allow user interrupts
+
+#include "MCMCfcds.h"
+
+using namespace std;
+using namespace scythe;
+
+/* This fn basically defined in MCMCregress.cc
+and should be moved to Scythe proper */
+static double lndigamma(double theta, double a, double b) {
+ double logf = a * log(b) - lngammafn(a) + -(a+1) * log(theta) +
+ -b/theta;
+ return logf;
+ //pow(b, a) / gammafn(a) * pow(theta, -(a+1)) * exp(-b/theta);
+}
+
+
+// Parameter-Expanded Latent Data Update for
+// 1d IRT. mjm, 2008-11-18
+template <typename RNGTYPE>
+double irt_W_update(Matrix<>& Z, const Matrix<>& X, const Matrix<>& theta,
+ const Matrix<>& eta, const double& alpha,
+ const double& px_a0, const double& px_b0,
+ const Matrix<>& etahat, const Matrix<>& thetahat,
+ rng<RNGTYPE>& stream) {
+
+ // define constants
+ const unsigned int J = theta.rows();
+ const unsigned int K = eta.rows();
+
+ double RSS=0.0;
+ int df=0;
+ // perform update from truncated Normal / standard Normals
+ for (unsigned int i = 0; i < J; ++i) {
+ for (unsigned int j = 0; j < K; ++j){
+ const double Z_mean = alpha*( -eta(j,0) + theta(i) * eta(j,1) );
+ const double Zhat = ( -etahat(j,0) + thetahat(i) * etahat(j,1) );
+ if (X(i,j) == 1) {
+ Z(i,j) = stream.rtbnorm_combo(Z_mean, alpha, 0);
+ ++df;
+ } else if (X(i,j) == 0) {
+ Z(i,j) = stream.rtanorm_combo(Z_mean, alpha, 0);
+ ++df;
+ } else {
+ Z(i,j) = stream.rnorm(Z_mean, std::pow(alpha, 2.0));
+ }
+ Z(i,j) /= alpha;
+ const double e = Z(i,j) - Zhat;
+ RSS += std::pow(e , 2.0);
+ }
+ }
+ // Liu and Wu 1999, p1272:
+ // draw a0 ~ IG(a0,b0). Then draw a1 ~ IG(nu0+a0RSS/2, a0+n/2)
+ const double c_post = (px_a0 + df) * 0.5;
+ const double d_post = (px_b0 + RSS) * 0.5;
+ double alpha1 = stream.rigamma(c_post,d_post);
+
+ //Rprintf("\nRSS: %5f alpha0: %5f alpha1: %5f\n",RSS,alpha,alpha1);
+ return(std::sqrt(alpha1/alpha));
+}
+
+template <typename RNGTYPE>
+void hirt_level0_metrop (Matrix<>& theta, Matrix<>& thetahat,
+ const Matrix<>& Z,
+ const Matrix<>& eta,
+ const Matrix<>& beta, const Matrix<>& Xj,
+ const double& sigma2,
+ const double& alpha,
+ rng<RNGTYPE>& stream){
+ //construct proposal sum(Zjk ~N( eta1k*(-eta0k + thetaj), alpha? ) )
+ //proposal logdensity
+ //current logdensity
+ //exp( prop - cur)
+ //runif & accept if > ratio
+ //note the acceptance
+}
+
+
+/* MCMCirt1d implementation. */
+template <typename RNGTYPE>
+void MCMCirtHier1d_impl (rng<RNGTYPE>& stream,
+ const Matrix<int>& X, Matrix<>& theta,
+ Matrix<>& eta,
+ Matrix<>& thetahat, Matrix<>& etahat,
+ const Matrix<>& ab0,
+ const Matrix<>& AB0, const Matrix<>& Xj,
+ Matrix<>& beta, const Matrix<>& b0,
+ const Matrix<>& B0,
+ const double c0, const double d0,
+ unsigned int burnin,
+ unsigned int mcmc, unsigned int thin, unsigned int verbose,
+ bool storea, bool storei,
+ double* sampledata, unsigned int samplesize,
+ bool chib, double* logmarglike,
+ bool px, const double px_a0, const double px_b0,
+ bool metromix
+ )
+{
+ // constants
+ const unsigned int J = X.rows(); // # subjects (justices, legislators)
+ const unsigned int K = X.cols(); // number of items (cases, roll calls)
+ const unsigned int L = Xj.cols(); // covariates Xj for mutheta
+ const unsigned int tot_iter = burnin + mcmc;
+ const unsigned int nsamp = mcmc / thin;
+
+ // storage matrices (col major order)
+ Matrix<> theta_store;
+ Matrix<> eta_store;
+ Matrix<> beta_store(nsamp,L);
+ Matrix<> sigma2_store(nsamp,1);
+ if (storea)
+ theta_store = Matrix<>(nsamp, J);
+
+ if (storei)
+ eta_store = Matrix<>(nsamp, K*2);
+
+ // starting values
+ Matrix<> Z(J, K);
+
+ // pre-compute
+ const Matrix<> AB0ab0 = AB0 * ab0;
+ const Matrix<> XpX = crossprod(Xj);
+ Matrix<> XpY;
+ double alpha;
+ if(!px) {
+ alpha = 1.0;
+ } else {
+ alpha = stream.rigamma(px_a0,px_b0);
+ }
+ double sigma2 = NormIGregress_sigma2_draw ( Xj, theta, beta, c0, d0, stream);
+
+ unsigned int count = 0;
+ // MCMC sampling occurs in this for loop
+ for (unsigned int iter = 0; iter < tot_iter; ++iter){
+ // sample latent data (Z) OR parameter-expanded data (W)
+ if(!px) {
+ irt_Z_update1(Z, X, theta,eta,stream);
+ } else {
+ // alpha here is the post variance of W
+ alpha = irt_W_update(Z, X, theta,eta,alpha,px_a0,px_b0,
+ etahat,thetahat,stream);
+ }
+
+ // Following Tierney(1994) it would be nice to 'restart' the chain
+ // using a metropolis jump on the entire level0, something like
+ // once every thin*10 iterations or so.
+ if( metromix && (iter % (thin*10) == 0)) {
+ // hirt_level0_metrop(theta, thetahat, Z, eta, etahat,
+ // beta, Xj, sigma2, alpha, stream);
+ } else {
+ // sample ability (ideal points) (theta)
+ hirt_theta_update1(theta, thetahat, Z, eta, beta, Xj, sigma2,
+ alpha, stream);
+
+ // sample item (case, bill) parameters (eta)
+ hirt_eta_update1(eta, etahat, Z, theta, AB0, AB0ab0,
+ alpha, stream);
+ }
+
+ XpY = t(Xj) * theta;
+ beta = NormNormregress_beta_draw(XpX, XpY, b0, B0, sigma2, stream);
+
+ // update level2 sigma2
+ sigma2 = NormIGregress_sigma2_draw (Xj, theta, beta, c0, d0, stream);
+
+ // print results to screen
+ if (verbose > 0 && iter % verbose == 0) {
+ Rprintf("\n\nMCMCirt1d iteration %i of %i \n", (iter+1), tot_iter);
+ //Rprintf("theta = \n");
+ //for (int j=0; j<J; ++j)
+ //Rprintf("%10.5f\n", theta[j]);
+ }
+
+ // store results
+ if ((iter >= burnin) && ((iter % thin == 0))) {
+
+ // store ideal points
+ if (storea)
+ theta_store(count, _) = theta;
+
+ // store bill parameters
+ if (storei)
+ eta_store(count, _) = t(eta);
+
+ beta_store(count, _) = t(beta);
+ sigma2_store(count,0) = sigma2;
+ // store beta
+ count++;
+ }
+
+ R_CheckUserInterrupt(); // allow user interrupts
+
+ } // end Gibbs loop
+
+ // return output
+ Matrix<> output;
+ if(! storei && storea) { // only theta
+ output = theta_store;
+ } else if (storei && ! storea){ // only eta
+ output = eta_store;
+ } else { // everything
+ output = cbind(theta_store, eta_store);
+ }
+ output = cbind(output, beta_store); // always return beta,
+ output = cbind(output, sigma2_store);// and sigma2.
+
+ for (unsigned int i = 0; i < samplesize; ++i)
+ sampledata[i] = output[i];
+
+ // BEGIN MARGINAL LIKELIHOOD COMPUTATION
+ if (chib == 1) {
+ // marginal likelihood calculation stuff starts here
+ const double sigma2star = meanc(sigma2_store)(0);
+ const Matrix<> thetastar = t(meanc(theta_store));
+ Matrix<> betastar = t(meanc(beta_store));
+ double sigma2fcdsum = 0.0;
+ XpY = t(Xj) * thetastar;
+ // second set of Gibbs scans
+ for (unsigned int iter = 0; iter < tot_iter; ++iter) {
+ double sigma2 = NormIGregress_sigma2_draw (Xj, thetastar, beta, c0, d0,
+ stream);
+
+ beta = NormNormregress_beta_draw (XpX, XpY, b0, B0, sigma2,
+ stream);
+ const Matrix<> e = gaxpy(Xj, (-1*beta), thetastar);
+ const Matrix<> SSE = crossprod (e);
+ const double c_post = (c0 + X.rows ()) * 0.5;
+ const double d_post = (d0 + SSE(0)) * 0.5;
+ sigma2fcdsum += lndigamma(sigma2star, c_post, d_post);
+
+ // print output to stdout
+ if(verbose > 0 && iter % verbose == 0) {
+ Rprintf("\n\nMCMCregress (reduced) iteration %i of %i \n",
+ (iter+1), tot_iter);
+ }
+ R_CheckUserInterrupt(); // allow user interrupts
+ } // end MCMC loop
+
+ double sigma2fcdmean = sigma2fcdsum / static_cast<double>(tot_iter);
+ const double sig2_inv = 1.0 / sigma2star;
+ const Matrix<> sig_beta = invpd (B0 + XpX * sig2_inv);
+ const Matrix<> betahat = sig_beta * gaxpy(B0, b0, XpY*sig2_inv);
+ double logbetafcd = lndmvn(betastar, betahat, sig_beta);
+ if(L==1) { logbetafcd = lndnorm(betastar[0], betahat[0],sig_beta[0]);}
+ Rprintf("logPost: beta %10.5f sigma2 %10.5f\n",logbetafcd,sigma2fcdmean);
+ // calculate loglikelihood at (betastar, sigma2star)
+ double sigmastar = sqrt(sigma2star);
+ Matrix<> xi = Xj * betastar;
+ double loglike = 0.0;
+ // is there a reason to do this and not lndmvn?
+ for (unsigned int i = 0; i < L; ++i) {
+ loglike += lndnorm(thetastar(i), xi(i), sigmastar);
+ }
+ Rprintf("logLike: %10.5f\n",loglike);
+
+ // calculate log prior ordinate
+ const double logpriorsig = (lndigamma(sigma2star, c0/2.0, d0/2.0));
+ double logpriorbeta = lndmvn(betastar, b0, invpd(B0));
+ if (L==1) { logpriorbeta = lndnorm(betastar[0], b0[0], 1.0/B0(0)); }
+ const double logprior = logpriorsig+logpriorbeta;
+ Rprintf("logPrior: %10.5f\n",logprior);
+
+ // put pieces together and print the marginal likelihood
+ logmarglike[0] = loglike + logprior - logbetafcd - (sigma2fcdmean);
+ Rprintf("logM: %10.5f\n",logmarglike[0]);
+
+ }
+}
+
+extern "C" {
+
+ void
+ MCMCirtHier1d(double* sampledata, const int* samplerow, const int* samplecol,
+ const int* Xdata, const int* Xrow, const int* Xcol,
+ const int* burnin, const int* mcmc, const int* thin,
+ const int *uselecuyer, const int *seedarray,
+ const int *lecuyerstream,
+ const int* verbose, const double* thetastartdata,
+ const int* thetastartrow, const int* thetastartcol,
+ const double* astartdata, const int* astartrow, const int* astartcol,
+ const double* bstartdata, const int* bstartrow, const int* bstartcol,
+ const double* ab0data, const int* ab0row, const int* ab0col,
+ const double* AB0data, const int* AB0row, const int* AB0col,
+ const double* Xjdata, const int* Xjrow, const int* Xjcol,
+ 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* storei, const int* storea,
+ double* logmarglikeholder, const int* chib, const int* px,
+ const double* px_a0, const double* px_b0)
+ {
+ // put together matrices
+ const Matrix<int> X(*Xrow, *Xcol, Xdata);
+ Matrix<> theta(*thetastartrow, *thetastartcol, thetastartdata);
+ Matrix<>thetahat(*thetastartrow, *thetastartcol, thetastartdata);
+ Matrix<> a(*astartrow, *astartcol, astartdata);
+ Matrix<> b(*bstartrow, *bstartcol, bstartdata);
+ const Matrix<> ab0(*ab0row, *ab0col, ab0data);
+ const Matrix<> AB0(*AB0row, *AB0col, AB0data);
+ Matrix<> eta = cbind(a,b);
+ Matrix<> etahat = cbind(a,b);
+ const Matrix<> Xj(*Xjrow, *Xjcol, Xjdata);
+ Matrix<> beta(*betastartrow, *betastartcol, betastartdata);
+ const Matrix<> b0(*b0row,*b0col,b0data);
+ const Matrix<> B0(*B0row,*B0col,B0data);
+ const int samplesize = (*samplerow) * (*samplecol);
+ const bool metromix = 0;
+
+ MCMCPACK_PASSRNG2MODEL(MCMCirtHier1d_impl, X, theta,
+ eta,
+ thetahat, etahat, ab0, AB0,
+ Xj, beta, b0, B0, *c0, *d0,
+ *burnin, *mcmc, *thin,
+ *verbose, *storea, *storei,
+ sampledata, samplesize,
+ *chib, logmarglikeholder, *px, *px_a0, *px_b0,
+ metromix);
+ }
+}
+
+#endif
diff --git a/src/MCMCirtKdRob.cc b/src/MCMCirtKdRob.cc
index dc99bf5..6400c59 100644
--- a/src/MCMCirtKdRob.cc
+++ b/src/MCMCirtKdRob.cc
@@ -1,3 +1,4 @@
+//////////////////////////////////////////////////////////////////////////
// MCMCirtKdRob.cc is C++ code to estimate a robust K-dimensional
// item response model
//
@@ -15,11 +16,14 @@
// PUBLIC LICENSE Version 2, June 1991. See the package LICENSE
// file for more information.
//
-// Copyright (C) 2004 Andrew D. Martin and Kevin M. Quinn
-//
// 2/18/2005 KQ
// 8/1/2007 ported to Scythe 1.0.2 KQ
//
+// Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+// Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+// and Jong Hee Park
+//////////////////////////////////////////////////////////////////////////
+
#ifndef MCMCIRTKDROB_CC
diff --git a/src/MCMClogit.cc b/src/MCMClogit.cc
index dfceae4..16396cb 100644
--- a/src/MCMClogit.cc
+++ b/src/MCMClogit.cc
@@ -1,3 +1,4 @@
+//////////////////////////////////////////////////////////////////////////
// MCMClogit.cc is C++ code to estimate a logistic regression model with
// a multivariate normal prior
//
@@ -14,12 +15,15 @@
// 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/25/2004 KQ
//
-//
+// Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+// Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+// and Jong Hee Park
+//////////////////////////////////////////////////////////////////////////
+
+
#ifndef MCMCLOGIT_CC
#define MCMCLOGIT_CC
diff --git a/src/MCMClogituserprior.cc b/src/MCMClogituserprior.cc
index 6808b8e..43d1971 100644
--- a/src/MCMClogituserprior.cc
+++ b/src/MCMClogituserprior.cc
@@ -1,3 +1,4 @@
+//////////////////////////////////////////////////////////////////////////
// MCMClogituserprior.cc samples from the posterior distribution of a
// logistic regression model with user-written prior coded in R using a
// random walk Metropolis algorithm
@@ -16,11 +17,13 @@
// PUBLIC LICENSE Version 2, June 1991. See the package LICENSE
// file for more information.
//
-// Copyright (C) 2004 Andrew D. Martin and Kevin M. Quinn
-//
// KQ 8/17/2005 (based on current version of MCMCmetrop1R.cc)
//
-//
+// Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+// Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+// and Jong Hee Park
+//////////////////////////////////////////////////////////////////////////
+
#ifndef MCMCLOGITUSERPRIOR_CC
diff --git a/src/MCMCmetrop1R.cc b/src/MCMCmetrop1R.cc
index a7be8ac..e2585ef 100644
--- a/src/MCMCmetrop1R.cc
+++ b/src/MCMCmetrop1R.cc
@@ -1,3 +1,4 @@
+//////////////////////////////////////////////////////////////////////////
// MCMCmetrop1R.cc samples from a user-written posterior code in R using a
// random walk Metropolis algorithm
//
@@ -15,10 +16,13 @@
// PUBLIC LICENSE Version 2, June 1991. See the package LICENSE
// file for more information.
//
-// Copyright (C) 2004 Andrew D. Martin and Kevin M. Quinn
-//
// KQ 6/24/2004
// updated to work with new Scythe and RNGs ADM 7/24/2004
+//
+// Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+// Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+// and Jong Hee Park
+//////////////////////////////////////////////////////////////////////////
#ifndef MCMCMETROP1R_CC
#define MCMCMETROP1R_CC
@@ -50,16 +54,20 @@ double user_fun_eval(SEXP fun, SEXP theta, SEXP myframe) {
if(!isEnvironment(myframe)) error("myframe must be an environment");
PROTECT(R_fcall = lang2(fun, R_NilValue));
SETCADR(R_fcall, theta);
- SEXP funval = eval(R_fcall, myframe);
+ SEXP funval;
+ PROTECT(funval = eval(R_fcall, myframe));
+
if (!isReal(funval)) error("`fun' must return a double");
double fv = REAL(funval)[0];
- UNPROTECT(1);
+ if (fv == R_PosInf) error("`fun' returned +Inf");
+ if (R_IsNaN(fv) || R_IsNA(fv)) error("`fun' returned NaN or NA");
+ UNPROTECT(2);
return fv;
}
template <typename RNGTYPE>
-void MCMCmetrop1R_impl (rng<RNGTYPE>& stream, SEXP fun,
- SEXP theta, SEXP myframe,
+void MCMCmetrop1R_impl (rng<RNGTYPE>& stream, SEXP& fun,
+ SEXP& theta, SEXP& myframe,
unsigned int burnin,
unsigned int mcmc, unsigned int thin,
unsigned int verbose, bool logfun,
@@ -97,7 +105,7 @@ void MCMCmetrop1R_impl (rng<RNGTYPE>& stream, SEXP fun,
// put theta_can_M into a SEXP
SEXP theta_can;
- theta_can = PROTECT(allocVector(REALSXP, npar));
+ PROTECT(theta_can = allocVector(REALSXP, npar));
for (unsigned int i = 0; i < npar; ++i) {
REAL(theta_can)[i] = theta_can_M(i);
}
@@ -109,11 +117,15 @@ void MCMCmetrop1R_impl (rng<RNGTYPE>& stream, SEXP fun,
const double ratio = std::exp(userfun_can - userfun_cur);
if (stream() < ratio) {
- theta = theta_can;
+ for (unsigned int i = 0; i < npar; ++i) {
+ REAL(theta)[i] = theta_can_M(i);
+ }
+ // theta = theta_can;
theta_M = theta_can_M;
userfun_cur = userfun_can;
++accepts;
}
+ UNPROTECT(1);
// store values in matrices
if ((iter%thin) == 0 && iter >= burnin) {
@@ -124,27 +136,27 @@ void MCMCmetrop1R_impl (rng<RNGTYPE>& stream, SEXP fun,
if (verbose && iter % verbose == 0) {
Rprintf("MCMCmetrop1R iteration %i of %i \n", (iter+1), tot_iter);
+ Rprintf("function value = %10.5f\n", userfun_cur);
Rprintf("theta = \n");
for (unsigned int i = 0; i < npar; ++i)
Rprintf("%10.5f\n", REAL(theta)[i]);
- Rprintf("function value = %10.5f\n", userfun_cur);
Rprintf("Metropolis acceptance rate = %3.5f\n\n",
static_cast<double>(accepts) / static_cast<double>(iter+1));
}
- UNPROTECT(1);
R_CheckUserInterrupt(); // allow user interrupts
}
// put the sample into a SEXP and return it
- sample_SEXP = PROTECT(allocMatrix(REALSXP, nsamp, npar));
+ //sample_SEXP = PROTECT(allocMatrix(REALSXP, nsamp, npar));
for (unsigned int i = 0; i < nsamp; ++i) {
for (unsigned int j = 0; j < npar; ++j) {
REAL(sample_SEXP)[i + nsamp * j] = sample(i,j);
}
}
- UNPROTECT(1);
+ //UNPROTECT(1);
+
// print the the acceptance rate to the console in a way that
// everyone (even Windows users) can see
@@ -179,13 +191,16 @@ extern "C" {
Matrix <> propvar (propvar_nc, propvar_nr, propvar_data);
propvar = t(propvar);
+ const unsigned int npar = length(theta);
+ const unsigned int nsamp = INTEGER(mcmc_R)[0] / INTEGER(thin_R)[0];
SEXP sample_SEXP;
+ PROTECT(sample_SEXP = allocMatrix(REALSXP, nsamp, npar));
MCMCPACK_PASSRNG2MODEL(MCMCmetrop1R_impl, fun, theta, myframe,
INTEGER(burnin_R)[0], INTEGER(mcmc_R)[0],
INTEGER(thin_R)[0],
INTEGER(verbose)[0], INTEGER(logfun)[0],
- propvar, sample_SEXP)
-
+ propvar, sample_SEXP);
+ UNPROTECT(1);
// return the sample
return sample_SEXP;
diff --git a/src/MCMCmixfactanal.cc b/src/MCMCmixfactanal.cc
index 207cf80..c049fa3 100644
--- a/src/MCMCmixfactanal.cc
+++ b/src/MCMCmixfactanal.cc
@@ -1,3 +1,4 @@
+//////////////////////////////////////////////////////////////////////////
// MCMCmixfactanal.cc is C++ code to estimate a mixed data
// factor analysis model
//
@@ -15,13 +16,16 @@
// PUBLIC LICENSE Version 2, June 1991. See the package LICENSE
// file for more information.
//
-// Copyright (C) 2004 Andrew D. Martin and Kevin M. Quinn
-//
// 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
//
+// Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+// Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+// and Jong Hee Park
+//////////////////////////////////////////////////////////////////////////
+
#ifndef MCMCMIXFACTANAL_CC
diff --git a/src/MCMCmnl.h b/src/MCMCmnl.h
index d5a2ec5..a98f6ce 100644
--- a/src/MCMCmnl.h
+++ b/src/MCMCmnl.h
@@ -1,3 +1,4 @@
+//////////////////////////////////////////////////////////////////////////
// MCMCmnl.h contains multinomial logit functions called by both the
// metropolis hastings and slice sampling implemenetation of MCMCmnl,
// such as a function that returns the log posterior.
@@ -16,9 +17,13 @@
// PUBLIC LICENSE Version 2, June 1991. See the package LICENSE
// file for more information.
//
-// Copyright (C) 2007 Andrew D. Martin and Kevin M. Quinn
-//
// DBP 7/27/2007
+//
+// Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+// Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+// and Jong Hee Park
+//////////////////////////////////////////////////////////////////////////
+
#ifndef MCMCMNL_H
#define MCMCMNL_H
diff --git a/src/MCMCmnlMH.cc b/src/MCMCmnlMH.cc
index b4bf93e..8a5616d 100644
--- a/src/MCMCmnlMH.cc
+++ b/src/MCMCmnlMH.cc
@@ -1,3 +1,4 @@
+//////////////////////////////////////////////////////////////////////////
// MCMCmnlMH.cc samples from the posterior distribution of a multinomial
// logit model using a random walk Metropolis algorithm.
//
@@ -19,11 +20,15 @@
// 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)
// 7/27/2007 DBP ported to scythe 1.0
+//
+// Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+// Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+// and Jong Hee Park
+//////////////////////////////////////////////////////////////////////////
+
#ifndef MCMCMNLMH_CC
#define MCMCMNLMH_CC
diff --git a/src/MCMCmnlslice.cc b/src/MCMCmnlslice.cc
index df9385b..892d43c 100644
--- a/src/MCMCmnlslice.cc
+++ b/src/MCMCmnlslice.cc
@@ -1,3 +1,4 @@
+//////////////////////////////////////////////////////////////////////////
// MCMCmnlslice.cc DESCRIPTION HERE
//
// The initial version of this file was generated by the
@@ -18,12 +19,16 @@
// 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
//
// 7/28/07 DBP ported to scythe 1.0
+//
+// Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+// Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+// and Jong Hee Park
+//////////////////////////////////////////////////////////////////////////
+
#ifndef MCMCMNLSLICE_CC
#define MCMCMNLSLICE_CC
diff --git a/src/MCMCoprobit.cc b/src/MCMCoprobit.cc
index 6b10749..7cf4192 100644
--- a/src/MCMCoprobit.cc
+++ b/src/MCMCoprobit.cc
@@ -1,3 +1,4 @@
+//////////////////////////////////////////////////////////////////////////
// MCMCoprobit.cc is C++ code to estimate a ordinalprobit regression
// model with a multivariate normal prior
//
@@ -15,12 +16,16 @@
// 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
// updated to Scythe 1.0.X 7/10/2007 ADM
// Albert and Chib method added 9/20/2007 JHP
+//
+// Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+// Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+// and Jong Hee Park
+//////////////////////////////////////////////////////////////////////////
+
#ifndef MCMCOPROBIT_CC
#define MCMCOPROBIT_CC
diff --git a/src/MCMCordfactanal.cc b/src/MCMCordfactanal.cc
index d0b951a..2056e6c 100644
--- a/src/MCMCordfactanal.cc
+++ b/src/MCMCordfactanal.cc
@@ -1,3 +1,4 @@
+//////////////////////////////////////////////////////////////////////////
// MCMCordfactanal.cc is C++ code to estimate an ordinal data
// factor analysis model
//
@@ -15,13 +16,16 @@
// PUBLIC LICENSE Version 2, June 1991. See the package LICENSE
// file for more information.
//
-// Copyright (C) 2004 Andrew D. Martin and Kevin M. Quinn
-//
// 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
//
+// Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+// Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+// and Jong Hee Park
+//////////////////////////////////////////////////////////////////////////
+
#ifndef MCMCORDFACTANAL_CC
#define MCMCORDFACTANAL_CC
diff --git a/src/MCMCpoisson.cc b/src/MCMCpoisson.cc
index aebf7f5..8b3851b 100644
--- a/src/MCMCpoisson.cc
+++ b/src/MCMCpoisson.cc
@@ -1,3 +1,4 @@
+//////////////////////////////////////////////////////////////////////////
// MCMCpoisson.cc is C++ code to estimate a Poisson regression model with
// a multivariate normal prior
//
@@ -15,11 +16,15 @@
// 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
// updated to Scythe 1.0.X 7/7/2007 ADM
//
+// Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+// Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+// and Jong Hee Park
+//////////////////////////////////////////////////////////////////////////
+
+
#ifndef MCMCPOISSON_CC
#define MCMCPOISSON_CC
diff --git a/src/MCMCpoissonChangepoint.cc b/src/MCMCpoissonChangepoint.cc
index da68951..17bc6dd 100644
--- a/src/MCMCpoissonChangepoint.cc
+++ b/src/MCMCpoissonChangepoint.cc
@@ -1,3 +1,4 @@
+//////////////////////////////////////////////////////////////////////////
// MCMCpoissonChange.cc is C++ code to estimate a Poisson changepoint model
// with a gamma prior
//
@@ -10,7 +11,10 @@
// PUBLIC LICENSE Version 2, June 1991. See the package LICENSE
// file for more information.
//
-// Copyright (C) 2004 Andrew D. Martin and Kevin M. Quinn
+// Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+// Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+// and Jong Hee Park
+//////////////////////////////////////////////////////////////////////////
#ifndef MCMCPOISSONCHANGEPOINT_CC
#define MCMCPOISSONCHANGEPOINT_CC
diff --git a/src/MCMCprobit.cc b/src/MCMCprobit.cc
index fad1e0e..d4c954c 100644
--- a/src/MCMCprobit.cc
+++ b/src/MCMCprobit.cc
@@ -1,3 +1,4 @@
+//////////////////////////////////////////////////////////////////////////
// MCMCprobit.cc is C++ code to estimate a probit regression model with
// a multivariate normal prior
//
@@ -15,11 +16,14 @@
// 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
// updated to Scythe 1.0.X 7/7/2007 ADM
//
+// Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+// Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+// and Jong Hee Park
+//////////////////////////////////////////////////////////////////////////
+
#ifndef MCMCPROBIT_CC
#define MCMCPROBIT_CC
diff --git a/src/MCMCprobitres.cc b/src/MCMCprobitres.cc
index caa8079..e893c83 100644
--- a/src/MCMCprobitres.cc
+++ b/src/MCMCprobitres.cc
@@ -1,3 +1,4 @@
+//////////////////////////////////////////////////////////////////////////
// MCMCprobitres.cc is a program that simulates draws from the posterior
// density of a probit regression model and returns latent residuals.
//
@@ -15,10 +16,14 @@
// 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
// updated to Scythe 1.0.X 7/28/2007 KQ
+//
+// Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+// Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+// and Jong Hee Park
+//////////////////////////////////////////////////////////////////////////
+
#ifndef MCMCPROBITRES_CC
#define MCMCPROBITRES_CC
diff --git a/src/MCMCregress.cc b/src/MCMCregress.cc
index a338bf4..32eb701 100644
--- a/src/MCMCregress.cc
+++ b/src/MCMCregress.cc
@@ -1,3 +1,4 @@
+//////////////////////////////////////////////////////////////////////////
// MCMCregress.cc is a program that simualates draws from the posterior
// density of a linear regression model with Gaussian errors.
//
@@ -19,8 +20,6 @@
// 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 Fri Jul 23 15:07:21 2004
//
// ADM and KQ 10/10/2002 [ported to Scythe0.3]
@@ -29,6 +28,11 @@
// ADM 7/22/04 [modified to work with new Scythe and rngs]
// DBP 7/1/07 [ported to scythe 1.0.x]
//
+// Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+// Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+// and Jong Hee Park
+//////////////////////////////////////////////////////////////////////////
+
#ifndef MCMCREGRESS_CC
#define MCMCREGRESS_CC
diff --git a/src/MCMCrng.h b/src/MCMCrng.h
index d4207b8..8e4664a 100644
--- a/src/MCMCrng.h
+++ b/src/MCMCrng.h
@@ -1,3 +1,4 @@
+//////////////////////////////////////////////////////////////////////////
// MCMCrng.h is the header file for MCMCrng.cc. It contains
// functions used to handle random number generator streams.
//
@@ -14,10 +15,13 @@
// 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
//
// ADM 7/22/04
+//
+// Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+// Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+// and Jong Hee Park
+//////////////////////////////////////////////////////////////////////////
#ifndef MCMCRNG_H
#define MCMCRNG_H
diff --git a/src/MCMCtobit.cc b/src/MCMCtobit.cc
index c80654f..c2c2cab 100644
--- a/src/MCMCtobit.cc
+++ b/src/MCMCtobit.cc
@@ -1,3 +1,4 @@
+//////////////////////////////////////////////////////////////////////////
// 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.
@@ -20,12 +21,16 @@
// 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 Tue Sep 14 00:50:08 2004
// ADM and KQ 10/10/2002 [ported to Scythe0.3]
// BG 09/18/2004 [updated to new specification, added above censoring]
// ADM 7/7/2007 [updated to Scythe 1.0.X]
+//
+// Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+// Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+// and Jong Hee Park
+//////////////////////////////////////////////////////////////////////////
+
#ifndef MCMCTOBIT_CC
#define MCMCTOBIT_CC
--
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