[r-cran-mcmcpack] 05/90: Imported Upstream version 0.6-4
Andreas Tille
tille at debian.org
Fri Dec 16 09:07:33 UTC 2016
This is an automated email from the git hooks/post-receive script.
tille pushed a commit to branch master
in repository r-cran-mcmcpack.
commit fb01eda9df453aad7fa527d4932957471b3a58dd
Author: Andreas Tille <tille at debian.org>
Date: Fri Dec 16 08:07:02 2016 +0100
Imported Upstream version 0.6-4
---
DESCRIPTION | 8 ++---
HISTORY | 39 ++++++++++++++++++++++
R/MCMCdynamicEI.R | 4 +--
R/MCMCfactanal.R | 5 +--
R/MCMChierEI.R | 4 +--
R/MCMCirt1d.R | 15 ++++++---
R/MCMCirtKd.R | 6 ++--
R/MCMClogit.R | 2 +-
R/MCMCmetrop1R.R | 6 ++--
R/MCMCmixfactanal.R | 8 ++---
R/MCMCmnl.R | 2 +-
R/MCMCoprobit.R | 4 +--
R/MCMCordfactanal.R | 19 ++++++-----
R/MCMCpanel.R | 2 +-
R/MCMCpoisson.R | 2 +-
R/MCMCprobit.R | 2 +-
R/MCMCregress.R | 2 +-
R/MCMCtobit.R | 2 +-
R/distn.R | 30 +++++++++--------
R/hidden.R | 41 ++++++++++++++++-------
man/MCMCdynamicEI.Rd | 17 +++++-----
man/MCMCfactanal.Rd | 12 ++++---
man/MCMChierEI.Rd | 9 ++---
man/MCMCirt1d.Rd | 26 +++++++++++----
man/MCMCirtKd.Rd | 89 +++++++++++++++++++++++++++++---------------------
man/MCMClogit.Rd | 9 ++---
man/MCMCmetrop1R.Rd | 25 +++++++++-----
man/MCMCmixfactanal.Rd | 12 ++++---
man/MCMCmnl.Rd | 17 +++++-----
man/MCMCoprobit.Rd | 7 ++--
man/MCMCordfactanal.Rd | 10 +++---
man/MCMCpanel.Rd | 7 ++--
man/MCMCpoisson.Rd | 7 ++--
man/MCMCprobit.Rd | 7 ++--
man/MCMCregress.Rd | 12 +++----
man/MCMCtobit.Rd | 12 ++++---
man/invgamma.Rd | 15 ++++++---
man/iwishart.Rd | 10 ++++--
man/wishart.Rd | 8 +++--
src/MCMCdynamicEI.cc | 9 +++--
src/MCMCfactanal.cc | 7 ++--
src/MCMChierEI.cc | 7 ++--
src/MCMCirt1d.cc | 8 +++--
src/MCMClogit.cc | 6 ++--
src/MCMCmetrop1R.cc | 5 +--
src/MCMCmixfactanal.cc | 5 +--
src/MCMCmnlMH.cc | 8 +++--
src/MCMCmnlslice.cc | 5 +--
src/MCMCoprobit.cc | 5 ++-
src/MCMCordfactanal.cc | 15 ++++++---
src/MCMCpanel.cc | 14 ++++----
src/MCMCpoisson.cc | 7 ++--
src/MCMCprobit.cc | 5 +--
src/MCMCprobitres.cc | 5 +--
src/MCMCregress.cc | 5 +--
src/MCMCtobit.cc | 6 ++--
src/Makevars | 2 +-
src/lecuyer.h | 7 ++--
58 files changed, 405 insertions(+), 240 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index 46fc5b5..2446cdb 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,11 +1,11 @@
Package: MCMCpack
-Version: 0.6-1
-Date: 2005-2-18
+Version: 0.6-4
+Date: 2005-3-17
Title: Markov chain Monte Carlo (MCMC) Package
Author: Andrew D. Martin <admartin at wustl.edu>, and
Kevin M. Quinn <kevin_quinn at harvard.edu>
Maintainer: Andrew D. Martin <admartin at wustl.edu>
-Depends: R (>= 2.0.0), coda (>= 0.8-1), MASS
+Depends: R (>= 2.0.1), coda (>= 0.9-2), MASS
Description: This package contains functions for posterior
simulation for a number of statistical models. All simulation
is done in compiled C++ written in the Scythe Statistical
@@ -17,4 +17,4 @@ Description: This package contains functions for posterior
sampling algorithm, and tools for visualization.
License: GPL version 2 or newer
URL: http://mcmcpack.wustl.edu
-Packaged: Sun Feb 20 22:22:48 2005; hornik
+Packaged: Thu Mar 17 09:42:26 2005; adm
diff --git a/HISTORY b/HISTORY
index fa62d93..dca215b 100644
--- a/HISTORY
+++ b/HISTORY
@@ -2,6 +2,45 @@
// Changes and Bug Fixes
//
+0.6-3 to 0.6-4
+ * fixed bug with verbose in MCMCmetrop1R pointed out by a referee for
+ Rnews
+ * added a little clarification to the lecuyer.h file about the dual-
+ licensing scheme Pierre L'Ecuyer agreed to on 7 August 2004.
+
+0.6-2 to 0.6-3
+ * cleaned up the docs for MCMCirtKd
+ * MCMCirtKd now fits a model with a difficulty parameter rather than
+ a negative difficulty (easiness ???) parameter
+ * MCMCirtKd and MCMCordfactanal now do a better job of dropping
+ variables that have 0 variance.
+ * fixed Senate examples for MCMCirt1d and MCMCirtKd. EVENTUALLY THE
+ SENATE DATA SHOULD BE FIXED DIRECTLY SO THAT THERE ARE NO DUPLICATE
+ NAMES IN THE MEMBER VARIABLE.
+ * users can interrupt all estimation algorithms now with CTL-C
+ * changed the behavior of the verbose switch. Now if verbose is greater
+ than 0 output is printed every verboseth iteration.
+ * changed the settings when mcmc objects are created so that
+ "start=burnin+1" and "end=mcmc+burnin"
+ * fixed documentation of MCMCmetrop1R to make clear how data are passed
+ to the log-posterior function
+ * fixed factor.score.start.check to not return an error when a vector of
+ starting values is passed
+ * improved the way factor.score.start.check enforces constraints
+
+0.6-1 to 0.6-2
+ * fixed documentation for rinvgamma() and dinvgamma() so that it is clear
+ that these functions take shape and scale parameters as arguments.
+ * fixed documentation and R code for riwish(), diwish(), rwish(), and
+ dwish() so that it is more clear how these distributions are parameterized.
+ * fixed line endings in MCMCtobit.cc
+ * changed how MCMCordfactanal and MCMCmixfactanal report MH acceptance
+ rates. They now report separate rates for each manifest variable.
+ * fixed a bug with starting values in the IRT and factor models. If no
+ starting values were passed, those created in factor.score.start.check()
+ did not include the constraints (hard or soft), so the check failed. Now
+ the constructed starting values meet the constraints.
+
0.5-2 to 0.6-1
* added a tobit model for a linear model with censoring in MCMCtobit()
* added multinomial logit model in MCMCmnl()
diff --git a/R/MCMCdynamicEI.R b/R/MCMCdynamicEI.R
index d349872..ac43d71 100644
--- a/R/MCMCdynamicEI.R
+++ b/R/MCMCdynamicEI.R
@@ -5,7 +5,7 @@
"MCMCdynamicEI" <-
function(r0, r1, c0, c1, burnin=5000, mcmc=50000,
- thin=1, verbose=FALSE, seed=NA,
+ thin=1, verbose=0, seed=NA,
W=0, a0=0.825, b0=0.0105, a1=0.825,
b1=0.0105, ...){
@@ -114,7 +114,7 @@
sample <- matrix(C.sample$samdata, C.sample$samrow, C.sample$samcol,
byrow=TRUE)
- output <- mcmc(data=sample, start=1, end=mcmc, thin=thin)
+ output <- mcmc(data=sample, start=(burnin+1), end=burnin+mcmc, thin=thin)
p0names <- paste("p0table", 1:ntables, sep="")
p1names <- paste("p1table", 1:ntables, sep="")
varnames(output) <- c(p0names, p1names, "sigma^2_0", "sigma^2_1")
diff --git a/R/MCMCfactanal.R b/R/MCMCfactanal.R
index 77c568f..670fa78 100644
--- a/R/MCMCfactanal.R
+++ b/R/MCMCfactanal.R
@@ -27,7 +27,7 @@
"MCMCfactanal" <-
function(x, factors, lambda.constraints=list(),
data=parent.frame(), burnin = 1000, mcmc = 20000,
- thin=1, verbose = FALSE, seed = NA,
+ thin=1, verbose = 0, seed = NA,
lambda.start = NA, psi.start = NA,
l0=0, L0=0, a0=0.001, b0=0.001,
store.scores = FALSE, std.var=TRUE, ... ) {
@@ -143,7 +143,8 @@
output.var <- diag(var(output.df))
output.df <- output.df[,output.var != 0]
- output <- mcmc(as.matrix(output.df), start=1, end=mcmc, thin=thin)
+ output <- mcmc(as.matrix(output.df), start=burnin+1, end=mcmc+burnin,
+ thin=thin)
## add constraint info so this isn't lost
attr(output, "constraints") <- lambda.constraints
diff --git a/R/MCMChierEI.R b/R/MCMChierEI.R
index d07fc3f..31d50e1 100644
--- a/R/MCMChierEI.R
+++ b/R/MCMChierEI.R
@@ -5,7 +5,7 @@
"MCMChierEI" <-
function(r0, r1, c0, c1, burnin=5000, mcmc=50000, thin=1,
- verbose=FALSE, seed=NA,
+ verbose=0, seed=NA,
m0=0, M0=2.287656,
m1=0, M1=2.287656,
a0=0.825, b0=0.0105,
@@ -121,7 +121,7 @@
sample <- matrix(C.sample$samdata, C.sample$samrow, C.sample$samcol,
byrow=TRUE)
- output <- mcmc(data=sample, start=1, end=mcmc, thin=thin)
+ output <- mcmc(data=sample, start=burnin+1, end=burnin+mcmc, thin=thin)
p0names <- paste("p0table", 1:ntables, sep="")
p1names <- paste("p1table", 1:ntables, sep="")
varnames(output) <- c(p0names, p1names, "mu0", "mu1", "sigma^2.0",
diff --git a/R/MCMCirt1d.R b/R/MCMCirt1d.R
index a9020e6..21b0d1e 100644
--- a/R/MCMCirt1d.R
+++ b/R/MCMCirt1d.R
@@ -6,15 +6,22 @@
"MCMCirt1d" <-
function(datamatrix, theta.constraints=list(), burnin = 1000,
- mcmc = 20000, thin=1, verbose = FALSE, seed = NA,
+ mcmc = 20000, thin=1, verbose = 0, seed = NA,
theta.start = NA, alpha.start = NA, beta.start = NA, t0 = 0,
- T0 = 1, ab0=0, AB0=.25, store.item = FALSE, ... ) {
+ T0 = 1, ab0=0, AB0=.25, store.item = FALSE,
+ drop.constant.items=TRUE, ... ) {
## checks
check.offset(list(...))
check.mcmc.parameters(burnin, mcmc, thin)
- ## check vote matrix and convert to work with C++ code
+ ## 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
@@ -164,7 +171,7 @@
sample <- matrix(posterior$sampledata, posterior$samplerow,
posterior$samplecol,
byrow=TRUE)
- output <- mcmc(data=sample, start=1, end=mcmc, thin=thin)
+ output <- mcmc(data=sample, start=burnin+1, end=burnin+mcmc, thin=thin)
if(store.item == FALSE) {
names <- theta.names
diff --git a/R/MCMCirtKd.R b/R/MCMCirtKd.R
index 4c9115d..855dc97 100644
--- a/R/MCMCirtKd.R
+++ b/R/MCMCirtKd.R
@@ -16,10 +16,10 @@
"MCMCirtKd" <-
function(datamatrix, dimensions, item.constraints=list(),
burnin = 1000, mcmc = 10000,
- thin=1, verbose = FALSE, seed = NA,
+ thin=1, verbose = 0, seed = NA,
alphabeta.start = NA, b0=0, B0=0,
store.item=FALSE, store.ability=TRUE,
- drop.constantvars=TRUE, ... ) {
+ drop.constant.items=TRUE, ... ) {
datamatrix <- as.matrix(datamatrix)
@@ -30,7 +30,7 @@
lambda.start=alphabeta.start,
l0=b0, L0=B0, store.lambda=store.item,
store.scores=store.ability,
- drop.constantvars=drop.constantvars,
+ drop.constantvars=drop.constant.items,
model="MCMCirtKd")
return(post)
}
diff --git a/R/MCMClogit.R b/R/MCMClogit.R
index 7ccc1f4..8aabcd9 100644
--- a/R/MCMClogit.R
+++ b/R/MCMClogit.R
@@ -9,7 +9,7 @@
"MCMClogit" <-
function(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
- thin=1, tune=1.1, verbose = FALSE, seed = NA, beta.start = NA,
+ thin=1, tune=1.1, verbose = 0, seed = NA, beta.start = NA,
b0 = 0, B0 = 0, ...) {
## checks
diff --git a/R/MCMCmetrop1R.R b/R/MCMCmetrop1R.R
index 356585d..4bc4da8 100644
--- a/R/MCMCmetrop1R.R
+++ b/R/MCMCmetrop1R.R
@@ -6,7 +6,7 @@
"MCMCmetrop1R" <- function(fun, theta.init,
burnin=500, mcmc=20000, thin=1,
- tune=1, verbose=TRUE, seed=NA, logfun=TRUE,
+ tune=1, verbose=0, seed=NA, logfun=TRUE,
optim.trace=0, optim.REPORT=10,
optim.maxit=500, ...){
@@ -53,7 +53,7 @@
sample <- .Call("MCMCmetrop1R_cc", fun, as.double(theta.init),
my.env, as.integer(burnin), as.integer(mcmc),
as.integer(thin),
- as.logical(verbose),
+ as.integer(verbose),
lecuyer=as.integer(lecuyer),
seedarray=as.integer(seed.array),
lecuyerstream=as.integer(lecuyer.stream),
@@ -62,7 +62,7 @@
PACKAGE="MCMCpack")
## turn sample into an mcmc object
- sample <- mcmc(data=sample, start=1, end = mcmc, thin=thin)
+ sample <- mcmc(data=sample, start=burnin+1, end=burnin+mcmc, thin=thin)
return(sample)
}
diff --git a/R/MCMCmixfactanal.R b/R/MCMCmixfactanal.R
index efa88b7..2cb04dc 100644
--- a/R/MCMCmixfactanal.R
+++ b/R/MCMCmixfactanal.R
@@ -28,7 +28,7 @@
"MCMCmixfactanal" <-
function(x, factors, lambda.constraints=list(),
data=parent.frame(), burnin = 1000, mcmc = 20000,
- thin=1, tune=NA, verbose = FALSE, seed = NA,
+ thin=1, tune=NA, verbose = 0, seed = NA,
lambda.start = NA, psi.start=NA,
l0=0, L0=0, a0=0.001, b0=0.001,
store.lambda=TRUE, store.scores=FALSE,
@@ -351,14 +351,14 @@
output.df <- as.data.frame(as.matrix(output))
output.var <- diag(var(output.df))
output.df <- output.df[,output.var != 0]
- output <- mcmc(as.matrix(output.df), start=1, end=mcmc, thin=thin)
+ output <- mcmc(as.matrix(output.df), start=burnin+1, end=burnin+mcmc,
+ thin=thin)
# add constraint info so this isn't lost
attr(output, "constraints") <- lambda.constraints
attr(output, "n.manifest") <- K
attr(output, "n.factors") <- factors
- attr(output, "accept.rate") <- posterior$accepts /
- ((posterior$burnin+posterior$mcmc)*n.ord.ge3)
+ attr(output, "accept.rates") <- t(accepts) / (posterior$burnin+posterior$mcmc)
attr(output,"title") <-
"MCMCpack Mixed Data Factor Analysis Posterior Density Sample"
diff --git a/R/MCMCmnl.R b/R/MCMCmnl.R
index 879da08..3827222 100644
--- a/R/MCMCmnl.R
+++ b/R/MCMCmnl.R
@@ -189,7 +189,7 @@
"MCMCmnl" <-
function(formula, baseline=NULL, data = parent.frame(),
burnin = 1000, mcmc = 10000, thin=1,
- mcmc.method = "MH", tune = 1.1, verbose = FALSE, seed = NA,
+ mcmc.method = "MH", tune = 1.1, verbose = 0, seed = NA,
beta.start = NA, b0 = 0, B0 = 0, ...) {
## checks
diff --git a/R/MCMCoprobit.R b/R/MCMCoprobit.R
index 635913e..f072c47 100644
--- a/R/MCMCoprobit.R
+++ b/R/MCMCoprobit.R
@@ -8,7 +8,7 @@
"MCMCoprobit" <-
function(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
- thin = 1, tune = NA, verbose = FALSE, seed = NA, beta.start = NA,
+ thin = 1, tune = NA, verbose = 0, seed = NA, beta.start = NA,
b0 = 0, B0 = 0, ...) {
@@ -128,7 +128,7 @@
sample <- matrix(posterior$sampledata, posterior$samplerow,
posterior$samplecol, byrow=TRUE)
sample <- sample[,c(1:K, (K+3):(K+ncat))]
- output <- mcmc(data=sample, start=1, end=mcmc, thin=thin)
+ output <- mcmc(data=sample, start=burnin+1, end=burnin+mcmc, thin=thin)
xnames <- c(X.names, paste("gamma", 2:(ncat-1), sep=""))
varnames(output) <- xnames
attr(output, "title") <- "MCMCoprobit Posterior Density Sample"
diff --git a/R/MCMCordfactanal.R b/R/MCMCordfactanal.R
index ed75b0c..11a7daf 100644
--- a/R/MCMCordfactanal.R
+++ b/R/MCMCordfactanal.R
@@ -26,7 +26,7 @@
"MCMCordfactanal" <-
function(x, factors, lambda.constraints=list(),
data=parent.frame(), burnin = 1000, mcmc = 20000,
- thin=1, tune=NA, verbose = FALSE, seed = NA,
+ thin=1, tune=NA, verbose = 0, seed = NA,
lambda.start = NA, l0=0, L0=0,
store.lambda=TRUE, store.scores=FALSE,
drop.constantvars=TRUE, ... ) {
@@ -55,9 +55,9 @@
if (is.matrix(x)){
if (drop.constantvars==TRUE){
x.col.var <- apply(x, 2, var, na.rm=TRUE)
- x <- x[,x.col.var!=0]
- x.row.var <- apply(x, 1, var, na.rm=TRUE)
- x <- x[x.row.var!=0,]
+ keep.inds <- x.col.var>0
+ keep.inds[is.na(keep.inds)] <- FALSE
+ x <- x[,keep.inds]
}
X <- as.data.frame(x)
xvars <- dimnames(X)[[2]]
@@ -129,7 +129,10 @@
holder <- form.factload.norm.prior(l0, L0, K, factors+1, X.names)
Lambda.prior.mean <- holder[[1]]
Lambda.prior.prec <- holder[[2]]
-
+ if (case.switch==2){# if MCMCirtKD make it a prior on the diff. param.
+ Lambda.prior.mean[,1] <- Lambda.prior.mean[,1] * -1
+ }
+
# seeds
seeds <- form.seeds(seed)
lecuyer <- seeds[[1]]
@@ -349,14 +352,14 @@
output.df <- as.data.frame(as.matrix(output))
output.var <- diag(var(output.df))
output.df <- output.df[,output.var != 0]
- output <- mcmc(as.matrix(output.df), start=1, end=mcmc, thin=thin)
+ output <- mcmc(as.matrix(output.df), start=burnin+1, end=burnin+mcmc,
+ thin=thin)
# add constraint info so this isn't lost
attr(output, "constraints") <- lambda.constraints
attr(output, "n.manifest") <- K
attr(output, "n.factors") <- factors
- attr(output, "accept.rate") <- posterior$accepts /
- ((posterior$burnin+posterior$mcmc)*K)
+ attr(output, "accept.rates") <- t(accepts) / (posterior$burnin+posterior$mcmc)
if(case.switch==1) {
attr(output,"title") <-
"MCMCpack Ordinal Data Factor Analysis Posterior Density Sample"
diff --git a/R/MCMCpanel.R b/R/MCMCpanel.R
index b8dc434..5b32ada 100644
--- a/R/MCMCpanel.R
+++ b/R/MCMCpanel.R
@@ -15,7 +15,7 @@
"MCMCpanel" <-
function(obs, Y, X, W, burnin = 1000, mcmc = 10000, thin = 5,
- verbose = FALSE, seed = NA, sigma2.start = NA,
+ verbose = 0, seed = NA, sigma2.start = NA,
D.start = NA, b0 = 0, B0 = 1, eta0, R0, nu0 = 0.001,
delta0 = 0.001, ...) {
diff --git a/R/MCMCpoisson.R b/R/MCMCpoisson.R
index 69e86c4..3eecf13 100644
--- a/R/MCMCpoisson.R
+++ b/R/MCMCpoisson.R
@@ -8,7 +8,7 @@
"MCMCpoisson" <-
function(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
- thin=1, tune=1.1, verbose = FALSE, seed = NA, beta.start = NA,
+ thin=1, tune=1.1, verbose = 0, seed = NA, beta.start = NA,
b0 = 0, B0 = 0, ...) {
## checks
diff --git a/R/MCMCprobit.R b/R/MCMCprobit.R
index ff15b94..2de612e 100644
--- a/R/MCMCprobit.R
+++ b/R/MCMCprobit.R
@@ -8,7 +8,7 @@
"MCMCprobit" <-
function(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
- thin = 1, verbose = FALSE, seed = NA, beta.start = NA,
+ thin = 1, verbose = 0, seed = NA, beta.start = NA,
b0 = 0, B0 = 0, bayes.resid=FALSE, ...) {
## checks
diff --git a/R/MCMCregress.R b/R/MCMCregress.R
index b32fb4d..d22cb8a 100644
--- a/R/MCMCregress.R
+++ b/R/MCMCregress.R
@@ -8,7 +8,7 @@
"MCMCregress" <-
function(formula, data=parent.frame(), burnin = 1000, mcmc = 10000,
- thin=1, verbose = FALSE, seed = NA, beta.start = NA,
+ thin=1, verbose = 0, seed = NA, beta.start = NA,
b0 = 0, B0 = 0, c0 = 0.001, d0 = 0.001, ...) {
# checks
diff --git a/R/MCMCtobit.R b/R/MCMCtobit.R
index 21476bb..633f3e4 100644
--- a/R/MCMCtobit.R
+++ b/R/MCMCtobit.R
@@ -1,7 +1,7 @@
"MCMCtobit" <-
function(formula, data=parent.frame(), below = 0.0, above = Inf,
burnin = 1000, mcmc = 10000,
- thin=1, verbose = FALSE, seed = NA, beta.start = NA,
+ thin=1, verbose = 0, seed = NA, beta.start = NA,
b0 = 0, B0 = 0, c0 = 0.001, d0 = 0.001, ...) {
# checks
diff --git a/R/distn.R b/R/distn.R
index 05ac9db..c69836d 100644
--- a/R/distn.R
+++ b/R/distn.R
@@ -120,7 +120,7 @@
"riwish" <-
function(v, S) {
- return(solve(rwish(v,S)))
+ return(solve(rwish(v,solve(S))))
}
# diwish evaluates the inverse Wishart pdf at positive definite
@@ -186,31 +186,33 @@
## Inverse Gamma
##
-# evaluate the inverse gamma density
-#
-# Kevin Rompala 5/6/2003
+## evaluate the inverse gamma density
+##
+## Kevin Rompala 5/6/2003
+## fixed KQ 3/8/2005
"dinvgamma" <-
- function(x, shape, rate = 1) {
+ function(x, shape, scale = 1) {
# error checking
- if(shape <= 0 | rate <=0) {
- stop("Shape or rate parameter negative in dinvgamma().\n")
+ if(shape <= 0 | scale <=0) {
+ stop("Shape or scale parameter negative in dinvgamma().\n")
}
alpha <- shape
- beta <- rate
+ beta <- scale
return(beta^alpha / gamma(alpha) * x^(-1*(alpha + 1)) * exp(-beta/x))
}
-# generate draws from the inverse gamma density (using
-# the gamma simulator)
-#
-# Kevin Rompala 5/6/2003
+## generate draws from the inverse gamma density (using
+## the gamma simulator)
+##
+## Kevin Rompala 5/6/2003
+## fixed KQ 3/8/2005
"rinvgamma" <-
- function(n, shape, rate = 1) {
- return(1 / rgamma(n, shape, rate))
+ function(n, shape, scale = 1) {
+ return(1 / rgamma(n, shape, scale))
}
##
diff --git a/R/hidden.R b/R/hidden.R
index 6a963e1..53f4245 100644
--- a/R/hidden.R
+++ b/R/hidden.R
@@ -254,11 +254,28 @@
N <- nrow(X)
## set value of theta.start
- if (is.na(theta.start)) {
+ if (max(is.na(theta.start))==1) {
theta.start <- factor.score.eigen.start(agree.mat(X), 1)
for (i in 1:factors){
theta.start[,i] <- prior.mean[i] + theta.start[,i] *
sqrt(1/prior.prec[i,i])
+
+ # make sure these are consistent with hard and soft constraints
+ for (j in 1:nrow(theta.start)){
+ if (eq.constraints[j,i] != -999){
+ if (theta.start[j,i] * eq.constraints[j,i] < 0){
+ theta.start[,i] <- -1*theta.start[,i]
+ }
+ }
+ if (theta.start[j,i] * ineq.constraints[j,i] < 0){
+ theta.start[,i] <- -1*theta.start[,i]
+ }
+ }
+ theta.start[eq.constraints[,i]!=-999,i] <-
+ eq.constraints[eq.constraints[,i]!=-999,i]
+ theta.start[ineq.constraints[,i]!=0,i] <-
+ abs(theta.start[ineq.constraints[,i]!=0,i]) *
+ ineq.constraints[ineq.constraints[,i]!=0,i]
}
}
else if(is.numeric(theta.start) && is.null(dim(theta.start))) {
@@ -471,17 +488,17 @@
# pull together the posterior density sample
"form.mcmc.object" <-
function(posterior.object, names, title) {
- holder <- matrix(posterior.object$sampledata,
- posterior.object$samplerow,
- posterior.object$samplecol,
- byrow=TRUE)
-
- output <- mcmc(data=holder, start=1,
- end=posterior.object$mcmc,
- thin=posterior.object$thin)
- varnames(output) <- as.list(names)
- attr(output,"title") <- title
- return(output)
+ holder <- matrix(posterior.object$sampledata,
+ posterior.object$samplerow,
+ posterior.object$samplecol,
+ byrow=TRUE)
+
+ output <- mcmc(data=holder, start=(posterior.object$burnin+1),
+ end=(posterior.object$burnin+posterior.object$mcmc),
+ thin=posterior.object$thin)
+ varnames(output) <- as.list(names)
+ attr(output,"title") <- title
+ return(output)
}
# form multivariate Normal prior
diff --git a/man/MCMCdynamicEI.Rd b/man/MCMCdynamicEI.Rd
index 1f000d6..1c35988 100644
--- a/man/MCMCdynamicEI.Rd
+++ b/man/MCMCdynamicEI.Rd
@@ -9,7 +9,7 @@
\usage{
MCMCdynamicEI(r0, r1, c0, c1, burnin=5000, mcmc=50000, thin=1,
- verbose=FALSE, seed=NA, W=0, a0=0.825,
+ verbose=0, seed=NA, W=0, a0=0.825,
b0=0.0105, a1=0.825, b1=0.0105, ...)
}
@@ -34,9 +34,10 @@ MCMCdynamicEI(r0, r1, c0, c1, burnin=5000, mcmc=50000, thin=1,
mcmc iterations must be divisible by this value.}
\item{verbose}{A switch which determines whether or not the progress of
- the sampler is printed to the screen. Information is printed
- if TRUE.}
-
+ 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
@@ -171,11 +172,11 @@ c1 <- (r0+r1) - c0
dtomogplot(r0, r1, c0, c1, delay=0.1)
## fit dynamic model
-post1 <- MCMCdynamicEI(r0,r1,c0,c1, mcmc=40000, thin=5, verbose=TRUE,
+post1 <- MCMCdynamicEI(r0,r1,c0,c1, mcmc=40000, thin=5, verbose=100,
seed=list(NA, 1))
## fit exchangeable hierarchical model
-post2 <- MCMChierEI(r0,r1,c0,c1, mcmc=40000, thin=5, verbose=TRUE,
+post2 <- MCMChierEI(r0,r1,c0,c1, mcmc=40000, thin=5, verbose=100,
seed=list(NA, 2))
p0meanDyn <- colMeans(post1)[1:n]
@@ -203,11 +204,11 @@ c1 <- (r0+r1) - c0
dtomogplot(r0, r1, c0, c1, delay=0.1)
## fit dynamic model
-post1 <- MCMCdynamicEI(r0,r1,c0,c1, mcmc=40000, thin=5, verbose=TRUE,
+post1 <- MCMCdynamicEI(r0,r1,c0,c1, mcmc=40000, thin=5, verbose=100,
seed=list(NA, 1))
## fit exchangeable hierarchical model
-post2 <- MCMChierEI(r0,r1,c0,c1, mcmc=40000, thin=5, verbose=TRUE,
+post2 <- MCMChierEI(r0,r1,c0,c1, mcmc=40000, thin=5, verbose=100,
seed=list(NA, 2))
p0meanDyn <- colMeans(post1)[1:n]
diff --git a/man/MCMCfactanal.Rd b/man/MCMCfactanal.Rd
index 158972c..3a00c4b 100644
--- a/man/MCMCfactanal.Rd
+++ b/man/MCMCfactanal.Rd
@@ -14,7 +14,7 @@
\usage{
MCMCfactanal(x, factors, lambda.constraints=list(),
data=parent.frame(), burnin = 1000, mcmc = 20000,
- thin=1, verbose = FALSE, seed = NA,
+ thin=1, verbose = 0, seed = NA,
lambda.start = NA, psi.start = NA,
l0=0, L0=0, a0=0.001, b0=0.001,
store.scores = FALSE, std.var=TRUE, ... )
@@ -47,8 +47,10 @@ MCMCfactanal(x, factors, lambda.constraints=list(),
iterations must be divisible by this value.}
\item{verbose}{A switch which determines whether or not the progress of
- the sampler is printed to the screen. If TRUE, the iteration number and
- the factor loadings and uniquenesses are printed to the screen.}
+ the sampler is printed to the screen. If \code{verbose} is greater
+ than 0 the iteration number and
+ the factor loadings and uniquenesses are printed to the screen every
+ \code{verbose}th iteration.}
\item{seed}{The seed for the random number generator. If NA, the Mersenne
Twister generator is used with default seed 12345; if an integer is
@@ -172,7 +174,7 @@ MCMCfactanal(x, factors, lambda.constraints=list(),
lambda.constraints=list(Examination=list(1,"+"),
Examination=list(2,"-"), Education=c(2,0),
Infant.Mortality=c(1,0)),
- verbose=FALSE, store.scores=FALSE, a0=1, b0=0.15,
+ verbose=0, store.scores=FALSE, a0=1, b0=0.15,
data=swiss, burnin=5000, mcmc=50000, thin=20)
plot(posterior)
summary(posterior)
@@ -187,7 +189,7 @@ MCMCfactanal(x, factors, lambda.constraints=list(),
lambda.constraints=list(Examination=list(1,"+"),
Examination=list(2,"-"), Education=c(2,0),
Infant.Mortality=c(1,0)),
- verbose=FALSE, store.scores=FALSE, a0=1, b0=0.15,
+ verbose=0, store.scores=FALSE, a0=1, b0=0.15,
burnin=5000, mcmc=50000, thin=20)
}
}
diff --git a/man/MCMChierEI.Rd b/man/MCMChierEI.Rd
index d1e010f..47918b2 100644
--- a/man/MCMChierEI.Rd
+++ b/man/MCMChierEI.Rd
@@ -9,7 +9,7 @@
\usage{
MCMChierEI(r0, r1, c0, c1, burnin=5000, mcmc=50000, thin=1,
- verbose=FALSE, seed=NA,
+ verbose=0, seed=NA,
m0=0, M0=2.287656, m1=0, M1=2.287656, a0=0.825, b0=0.0105,
a1=0.825, b1=0.0105, ...)
}
@@ -35,8 +35,9 @@ MCMChierEI(r0, r1, c0, c1, burnin=5000, mcmc=50000, thin=1,
mcmc iterations must be divisible by this value.}
\item{verbose}{A switch which determines whether or not the progress of
- the sampler is printed to the screen. Information is printed
- if TRUE.}
+ 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
@@ -171,7 +172,7 @@ c1 <- (r0+r1) - c0
tomogplot(r0, r1, c0, c1)
## fit exchangeable hierarchical model
-post <- MCMChierEI(r0,r1,c0,c1, mcmc=40000, thin=5, verbose=TRUE,
+post <- MCMChierEI(r0,r1,c0,c1, mcmc=40000, thin=5, verbose=100,
seed=list(NA, 1))
p0meanHier <- colMeans(post)[1:n]
diff --git a/man/MCMCirt1d.Rd b/man/MCMCirt1d.Rd
index 74ceede..7f0b222 100644
--- a/man/MCMCirt1d.Rd
+++ b/man/MCMCirt1d.Rd
@@ -19,9 +19,9 @@
\usage{
MCMCirt1d(datamatrix, theta.constraints=list(), burnin = 1000,
- mcmc = 20000, thin=1, verbose = FALSE, seed = NA, theta.start = NA,
+ mcmc = 20000, thin=1, verbose = 0, seed = NA, theta.start = NA,
alpha.start = NA, beta.start = NA, t0 = 0, T0 = 1, ab0=0, AB0=.25,
- store.item = FALSE, ... ) }
+ store.item = FALSE, drop.constant.items=TRUE, ... ) }
\arguments{
\item{datamatrix}{The matrix of data. Must be 0, 1, or missing values.
@@ -48,8 +48,9 @@ MCMCirt1d(datamatrix, theta.constraints=list(), burnin = 1000,
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 TRUE, the iteration number
- is printed to the screen every 100 iterations.}
+ 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
@@ -111,6 +112,10 @@ MCMCirt1d(datamatrix, theta.constraints=list(), burnin = 1000,
should only be used if the chain is thinned heavily, or for
applications with a small number of items}. By default, the
item parameters are not stored.}
+
+ \item{drop.constant.items}{A switch that determines whether or not
+ items that have no variation
+ should be deleted before fitting the model. Default = TRUE.}
\item{...}{further arguments to be passed}
}
@@ -187,7 +192,7 @@ MCMCirt1d(datamatrix, theta.constraints=list(), burnin = 1000,
posterior1 <- MCMCirt1d(t(SupremeCourt),
theta.constraints=list(Scalia="+", Ginsburg="-"),
B0.alpha=.2, B0.beta=.2,
- burnin=500, mcmc=100000, thin=20, verbose=TRUE,
+ burnin=500, mcmc=100000, thin=20, verbose=500,
store.item=TRUE)
geweke.diag(posterior1)
plot(posterior1)
@@ -195,11 +200,18 @@ MCMCirt1d(datamatrix, theta.constraints=list(), burnin = 1000,
## US Senate Example with equality constraints
data(Senate)
+ namestring <- as.character(Senate$member)
+ namestring[78] <- "CHAFEE1"
+ namestring[79] <- "CHAFEE2"
+ namestring[59] <- "SMITH.NH"
+ namestring[74] <- "SMITH.OR"
+ rownames(Senate) <- namestring
+
Sen.rollcalls <- Senate[,6:677]
- rownames(Sen.rollcalls) <- as.character(Senate$member)
+
posterior2 <- MCMCirt1d(Sen.rollcalls,
theta.constraints=list(KENNEDY=-2, HELMS=2),
- burnin=2000, mcmc=100000, thin=20, verbose=TRUE)
+ burnin=2000, mcmc=100000, thin=20, verbose=500)
geweke.diag(posterior2)
plot(posterior2)
summary(posterior2)
diff --git a/man/MCMCirtKd.Rd b/man/MCMCirtKd.Rd
index bc2296c..d218935 100644
--- a/man/MCMCirtKd.Rd
+++ b/man/MCMCirtKd.Rd
@@ -5,19 +5,19 @@
\description{
This function generates a posterior density sample from a
K-dimensional item response theory (IRT) model, with standard
- Normal priors on the subject abilities (ideal points), and
- Normal priors on the item parameters. The user
+ normal priors on the subject abilities (ideal points), and
+ normal priors on the item parameters. The user
supplies data and priors, and a sample from the posterior
- density is returned as an mcmc object, which can be
+ distribution is returned as an mcmc object, which can be
subsequently analyzed with functions provided in the coda
package.
}
\usage{
MCMCirtKd(datamatrix, dimensions, item.constraints=list(),
- burnin = 1000, mcmc = 10000, thin=1, verbose = FALSE, seed = NA,
+ burnin = 1000, mcmc = 10000, thin=1, verbose = 0, seed = NA,
alphabeta.start = NA, b0 = 0, B0=0, store.item = FALSE,
- store.ability=TRUE, drop.constantvars=TRUE, ... ) }
+ store.ability=TRUE, drop.constant.items=TRUE, ... ) }
\arguments{
\item{datamatrix}{The matrix of data. Must be 0, 1, or missing values.
@@ -32,16 +32,16 @@ MCMCirtKd(datamatrix, dimensions, item.constraints=list(),
rowname to be equal to c, \code{rowname=list(d,"+")} which will
constrain the dth item parameter for the item named rowname to be
positive, and\code{rowname=list(d, "-")} which will constrain the dth
- item parameter for the item named varname to be negative. If x is a
+ item parameter for the item named rowname to be negative. If x is a
matrix without row names defaults names of ``V1", ``V2", ... , etc
- will be used. In a d dimensional model, the first item parameter for
+ will be used. In a K dimensional model, the first item parameter for
item \eqn{i}{i} is the difficulty parameter (\eqn{\alpha_i}{alpha_i}),
the second item parameter is the discrimation parameter on dimension
- 1 (\eqn{\beta_{i,1}{beta_{i,1}}}), the third item parameter is the
+ 1 (\eqn{\beta_{i,1}}{beta_{i,1}}), the third item parameter is the
discrimation parameter on dimension 2
- (\eqn{\beta_{i,2}{beta_{i,2}}}), ..., and the (d+1)th item parameter
- is the discrimation parameter on dimension d
- (\eqn{\beta_{i,1}{beta_{i,1}}}).
+ (\eqn{\beta_{i,2}}{beta_{i,2}}), ..., and the (K+1)th item parameter
+ is the discrimation parameter on dimension K
+ (\eqn{\beta_{i,1}}{beta_{i,1}}).
The item difficulty parameters (\eqn{\alpha}{alpha}) should
generally not be constrained.
}
@@ -54,8 +54,9 @@ MCMCirtKd(datamatrix, dimensions, item.constraints=list(),
iterations must be divisible by this value.}
\item{verbose}{A switch which determines whether or not the progress of
- the sampler is printed to the screen. If TRUE, the iteration number and
- the subject abilities (ideal points) are printed to the screen.}
+ 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
@@ -88,38 +89,42 @@ MCMCirtKd(datamatrix, dimensions, item.constraints=list(),
is used as the prior mean for all items.}
\item{B0}{The prior precisions (inverse variances) of the
- independent Normal prior on the item parameters.
+ independent normal prior on the item parameters.
Can be either a scalar or a matrix of dimension
\eqn{(K+1) \times items}{(K+1) x items}.}
\item{store.item}{A switch that determines whether or not to
store the item parameters for posterior analysis.
- \emph{NOTE: This takes an enormous amount of memory, so
- should only be used if the chain is thinned heavily, or for
- applications with a small number of items}. By default, the
+ \emph{NOTE: In applications with many items
+ this takes an enormous amount of memory. If you have many items
+ and want to want to store the item parameters you may want to thin
+ the chain heavily}. By default, the
item parameters are not stored.}
\item{store.ability}{A switch that determines whether or not to store
- the subject abilities for posterior analysis. By default, the
- item parameters are all stored.}
-
- \item{drop.constantvars}{A switch that determines whether or not
- items and subjects that have no variation
+ the subject abilities for posterior analysis. \emph{NOTE: In
+ applications with many subjects this takes an enormous amount of
+ memory. If you have many subjects and want to want to store the ability
+ parameters you may want to thin the chain heavily}. By default, the
+ ability parameters are all stored.}
+
+ \item{drop.constant.items}{A switch that determines whether or not
+ items that have no variation
should be deleted before fitting the model. Default = TRUE.}
\item{...}{further arguments to be passed}
}
\value{
- An mcmc object that contains the posterior density sample. This
+ An mcmc object that contains the posterior sample. This
object can be summarized by functions provided by the coda package.
}
\details{
- \code{MCMCirtKd} simulates from the posterior density using
- standard Gibbs sampling using data augmentation (a Normal draw
- for the subject abilities, a multivariate Normal
- draw for the item parameters, and a truncated Normal draw for
+ \code{MCMCirtKd} simulates from the posterior distribution using
+ standard Gibbs sampling using data augmentation (a normal draw
+ for the subject abilities, a multivariate normal
+ draw for the item parameters, and a truncated normal draw for
the latent utilities). The simulation proper is done in
compiled C++ code to maximize efficiency. Please consult the
coda documentation for a comprehensive list of functions that
@@ -129,15 +134,14 @@ MCMCirtKd(datamatrix, dimensions, item.constraints=list(),
smaller than the typical default values in MCMCpack. This is
because fitting this model is extremely computationally
expensive. It does not mean that this small of a number of
- scans will yield good estimates. If the verbose option is
- chosen, output will be printed to the screen every fifty
- iterations. The priors of this model need to be proper for
+ scans will yield good estimates.
+ The priors of this model need to be proper for
identification purposes. The user is asked to provide prior
means and precisions \emph{(not variances)} for the item
parameters and the subject parameters.
The model takes the following form. We assume that each
- subject has an subject ability (ideal point) denoted
+ subject has an ability (ideal point) denoted
\eqn{\theta_j}{theta_j} \eqn{(K \times 1)}{(K x 1)},
and that each item has a difficulty
parameter \eqn{\alpha_i}{alpha_i} and discrimination parameter
@@ -145,17 +149,18 @@ MCMCirtKd(datamatrix, dimensions, item.constraints=list(),
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 +
+ dictated by an unobserved utility: \deqn{z_{i,j} = - \alpha_i +
\beta_i' \theta_j + \varepsilon_{i,j}}{z_ij = alpha_i +
- beta_i'*theta_j + epsilon_ij} Where the errors are assumed to be
- distributed standard Normal. The parameters of interest are
+ beta_i'*theta_j + epsilon_ij} Where the
+ \eqn{\varepsilon_{i,j}}{epsilon_ij}s are assumed to be
+ distributed standard normal. The parameters of interest are
the subject abilities (ideal points) and the item parameters.
We assume the following priors. For the subject abilities (ideal points)
- we assume independent standard Normal priors:
+ we assume independent standard normal priors:
\deqn{\theta_{j,k} \sim \mathcal{N}(0,1)}{theta_j,k ~ N(0, 1)}
These cannot be changed by the user.
- For each item parameter, we assume independent Normal priors:
+ For each item parameter, we assume independent normal priors:
\deqn{\left[\alpha_i, \beta_i \right]' \sim \mathcal{N}_{(K+1)}
(b_{0,i},B_{0,i})}{[alpha_i beta_i]' ~ N_(K+1) (b_0,i, B_0,i)}
Where \eqn{B_{0,i}}{B_0,i} is a diagonal matrix.
@@ -215,10 +220,18 @@ MCMCirtKd(datamatrix, dimensions, item.constraints=list(),
data(Senate)
- rownames(Senate) <- Senate$member
+ namestring <- as.character(Senate$member)
+ namestring[78] <- "CHAFEE1"
+ namestring[79] <- "CHAFEE2"
+ namestring[59] <- "SMITH.NH"
+ namestring[74] <- "SMITH.OR"
+ rownames(Senate) <- namestring
+
+ Sen.rollcalls <- Senate[,6:677]
+
# note that we need to transpose the data to get
# the bills on the rows
- posterior2 <- MCMCirtKd(Senate[,6:677], dimensions=2,
+ posterior2 <- MCMCirtKd(Sen.rollcalls, dimensions=2,
burnin=5000, mcmc=50000, thin=10,
item.constraints=list(rc2=list(2,"-"), rc2=c(3,0),
rc3=list(3,"-")),
diff --git a/man/MCMClogit.Rd b/man/MCMClogit.Rd
index d1ae4bd..dda7d63 100644
--- a/man/MCMClogit.Rd
+++ b/man/MCMClogit.Rd
@@ -12,7 +12,7 @@
\usage{
MCMClogit(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
- thin=1, tune=1.1, verbose = FALSE, seed = NA, beta.start = NA,
+ thin=1, tune=1.1, verbose = 0, seed = NA, beta.start = NA,
b0 = 0, B0 = 0, ...) }
\arguments{
@@ -34,9 +34,10 @@ MCMClogit(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
before using the posterior density sample for inference.}
\item{verbose}{A switch which determines whether or not the progress of
- the sampler is printed to the screen. If TRUE, the iteration number,
- the current beta vector, and the Metropolis acceptance rate are
- printed to the screen every 500 iterations.}
+ the sampler is printed to the screen. If \code{verbose} is greater
+ than 0 the iteration number,
+ the current beta vector, and the Metropolis acceptance rate are
+ printed to the screen every \code{verbose}th iteration.}
\item{seed}{The seed for the random number generator. If NA, the Mersenne
Twister generator is used with default seed 12345; if an integer is
diff --git a/man/MCMCmetrop1R.Rd b/man/MCMCmetrop1R.Rd
index 378c2a0..5b5f3c4 100644
--- a/man/MCMCmetrop1R.Rd
+++ b/man/MCMCmetrop1R.Rd
@@ -8,7 +8,7 @@
}
\usage{
MCMCmetrop1R(fun, theta.init, burnin = 500, mcmc = 20000, thin = 1,
- tune = 1, verbose = TRUE, seed=NA, logfun = TRUE,
+ tune = 1, verbose = 0, seed=NA, logfun = TRUE,
optim.trace = 0, optim.REPORT = 10, optim.maxit = 500, ...)
}
\arguments{
@@ -16,8 +16,8 @@ MCMCmetrop1R(fun, theta.init, burnin = 500, mcmc = 20000, thin = 1,
be a function defined in R and it should take a single
argument. Additional arguments can be passed implicitly by either
putting them in the global environment or by passing them as
- additional arguments to \code{MCMCmetrop1R()}. The examples below
- demonstrate both of these approaches. }
+ additional arguments to \code{MCMCmetrop1R()}. See the Details section
+ and the examples below. }
\item{theta.init}{Starting values for the sampling. Must be of the
appropriate dimension.}
@@ -33,9 +33,10 @@ MCMCmetrop1R(fun, theta.init, burnin = 500, mcmc = 20000, thin = 1,
\eqn{\theta}{theta}.}
\item{verbose}{A switch which determines whether or not the progress of
- the sampler is printed to the screen. If TRUE, the iteration number, the
+ the sampler is printed to the screen. If \code{verbose} is greater
+ than 0 the iteration number, the
\eqn{\theta}{theta} vector, the function value, and the Metropolis
- acceptance rate are sent to the screen every 500 iterations.}
+ acceptance rate are sent to the screen every \code{verbose}th iteration.}
\item{seed}{The seed for the random number generator. If NA, the Mersenne
Twister generator is used with default seed 12345; if an integer is
@@ -75,6 +76,14 @@ diagonal positive definite matrix formed from the \code{tune} and
mode. This last calculation is done via an initial call to
\code{optim}.
+Passing data into a (log)density function does not conform to standard R
+practice. The (log)density function should have only one argument -- the
+parameter vector. The data used by the (log)density need to be either in the
+global environment or defined in the call to \code{MCMCmetrop1R()}, which will
+create an environment that contains the data and then call the (log)density function in that environment. The names of the data objects
+are not checked, so be careful that these match the names used with the
+(log)density function. \emph{NOTE:} This interface will likely be changed in
+a future implementation of MCMCpack.
}
\value{
An mcmc object that contains the posterior density sample. This
@@ -102,7 +111,7 @@ mode. This last calculation is done via an initial call to
X=Xdata, y=yvector,
thin=1, mcmc=40000, burnin=500,
tune=c(1.5, 1.5, 1.5),
- verbose=TRUE, logfun=TRUE, optim.maxit=100)
+ verbose=500, logfun=TRUE, optim.maxit=100)
raftery.diag(post.samp)
plot(post.samp)
@@ -128,7 +137,7 @@ mode. This last calculation is done via an initial call to
post.samp <- MCMCmetrop1R(logitfun, theta.init=c(0,0,0),
thin=1, mcmc=40000, burnin=500,
tune=c(1.5, 1.5, 1.5),
- verbose=TRUE, logfun=TRUE, optim.maxit=100)
+ verbose=500, logfun=TRUE, optim.maxit=100)
raftery.diag(post.samp)
plot(post.samp)
@@ -159,7 +168,7 @@ mode. This last calculation is done via an initial call to
post.samp <- MCMCmetrop1R(negbinfun, theta.init=c(0,0,0,0), y=yy, X=XX,
thin=1, mcmc=35000, burnin=1000,
- tune=1.5, verbose=TRUE, logfun=TRUE,
+ tune=1.5, verbose=500, logfun=TRUE,
optim.maxit=500, seed=list(NA,1))
raftery.diag(post.samp)
plot(post.samp)
diff --git a/man/MCMCmixfactanal.Rd b/man/MCMCmixfactanal.Rd
index ae7b3eb..7434d38 100644
--- a/man/MCMCmixfactanal.Rd
+++ b/man/MCMCmixfactanal.Rd
@@ -15,7 +15,7 @@
\usage{
MCMCmixfactanal(x, factors, lambda.constraints=list(),
data=parent.frame(), burnin = 1000, mcmc = 20000,
- thin=1, tune=NA, verbose = FALSE, seed = NA,
+ thin=1, tune=NA, verbose = 0, seed = NA,
lambda.start = NA, psi.start=NA,
l0=0, L0=0, a0=0.001, b0=0.001,
store.lambda=TRUE, store.scores=FALSE,
@@ -64,8 +64,10 @@ MCMCmixfactanal(x, factors, lambda.constraints=list(),
strictly positive.}
\item{verbose}{A switch which determines whether or not the progress of
- the sampler is printed to the screen. If TRUE, the iteration number and
- the Metropolis-Hastings acceptance rate are printed to the screen.}
+ the sampler is printed to the screen. If \code{verbose} is great
+ than 0 the iteration number and
+ the Metropolis-Hastings acceptance rate are printed to the screen
+ every \code{verbose}th iteration.}
\item{seed}{The seed for the random number generator. If NA, the Mersenne
Twister generator is used with default seed 12345; if an integer is
@@ -242,7 +244,7 @@ post <- MCMCmixfactanal(~courts+barb2+prsexp2+prscorr2+gdpw2,
factors=1, data=PErisk,
lambda.constraints = list(courts=list(2,"-")),
burnin=5000, mcmc=1000000, thin=50,
- verbose=TRUE, L0=.25, store.lambda=TRUE,
+ verbose=500, L0=.25, store.lambda=TRUE,
store.scores=TRUE, tune=1.2)
plot(post)
summary(post)
@@ -284,7 +286,7 @@ post <- MCMCmixfactanal(~log.Price+log.MPG.city+
lambda.constraints=list(log.Horsepower=list(2,"+"),
log.Horsepower=c(3,0), weight=list(3,"+")),
factors=2,
- burnin=5000, mcmc=500000, thin=100, verbose=TRUE,
+ burnin=5000, mcmc=500000, thin=100, verbose=500,
L0=.25, tune=3.0)
plot(post)
summary(post)
diff --git a/man/MCMCmnl.Rd b/man/MCMCmnl.Rd
index 5f03fea..f9be68d 100644
--- a/man/MCMCmnl.Rd
+++ b/man/MCMCmnl.Rd
@@ -14,7 +14,7 @@
\usage{
MCMCmnl(formula, baseline = NULL, data = parent.frame(),
burnin = 1000, mcmc = 10000, thin = 1,
- mcmc.method = "MH", tune = 1.1, verbose = FALSE,
+ mcmc.method = "MH", tune = 1.1, verbose = 0,
seed = NA, beta.start = NA, b0 = 0, B0 = 0, ...) }
\arguments{
@@ -78,9 +78,10 @@ MCMCmnl(formula, baseline = NULL, data = parent.frame(),
\item{verbose}{A switch which determines whether or not the progress of
- the sampler is printed to the screen. If TRUE, the iteration number,
+ the sampler is printed to the screen. If \code{verbose} is greater
+ than 0 the iteration number,
the current beta vector, and the Metropolis acceptance rate are
- printed to the screen every 500 iterations. }
+ printed to the screen every \code{verbose}th iteration. }
\item{seed}{The seed for the random number generator. If NA, the Mersenne
Twister generator is used with default seed 12345; if an integer is
@@ -170,7 +171,7 @@ MCMCmnl(formula, baseline = NULL, data = parent.frame(),
choicevar(distVVD, "sqdist", "VVD") +
choicevar(distCDA, "sqdist", "CDA"),
baseline="D66", mcmc.method="MH", B0=0,
- verbose=TRUE, mcmc=100000, thin=10, tune=1.0,
+ verbose=500, mcmc=100000, thin=10, tune=1.0,
data=Nethvote)
plot(post1)
@@ -182,7 +183,7 @@ MCMCmnl(formula, baseline = NULL, data = parent.frame(),
post2<- MCMCmnl(vote ~
relig + class + income + educ + age + urban,
baseline="D66", mcmc.method="MH", B0=0,
- verbose=TRUE, mcmc=100000, thin=10, tune=0.5,
+ verbose=500, mcmc=100000, thin=10, tune=0.5,
data=Nethvote)
plot(post2)
@@ -198,7 +199,7 @@ MCMCmnl(formula, baseline = NULL, data = parent.frame(),
choicevar(distCDA, "sqdist", "CDA") +
relig + class + income + educ + age + urban,
baseline="D66", mcmc.method="MH", B0=0,
- verbose=TRUE, mcmc=100000, thin=10, tune=0.5,
+ verbose=500, mcmc=100000, thin=10, tune=0.5,
data=Nethvote)
plot(post3)
@@ -214,7 +215,7 @@ MCMCmnl(formula, baseline = NULL, data = parent.frame(),
choicevar(distCDA, "sqdist", "1") +
relig + class + income + educ + age + urban,
baseline="2", mcmc.method="MH", B0=0,
- verbose=TRUE, mcmc=100000, thin=10, tune=0.5,
+ verbose=500, mcmc=100000, thin=10, tune=0.5,
data=Nethvote)
@@ -254,7 +255,7 @@ MCMCmnl(formula, baseline = NULL, data = parent.frame(),
choicevar(xb, "x", "b") +
choicevar(xc, "x", "c") +
choicevar(xd, "x", "d") + z,
- baseline="c", verbose=TRUE,
+ baseline="c", verbose=500,
mcmc=100000, thin=10, tune=.85)
plot(post5)
diff --git a/man/MCMCoprobit.Rd b/man/MCMCoprobit.Rd
index 4d13702..b1a4a04 100644
--- a/man/MCMCoprobit.Rd
+++ b/man/MCMCoprobit.Rd
@@ -12,7 +12,7 @@
\usage{
MCMCoprobit(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
- thin=1, tune = NA, verbose = FALSE, seed = NA, beta.start = NA,
+ thin=1, tune = NA, verbose = 0, seed = NA, beta.start = NA,
b0 = 0, B0 = 0, ...) }
\arguments{
@@ -32,9 +32,10 @@ MCMCoprobit(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
number of categories in the response variable.}
\item{verbose}{A switch which determines whether or not the progress of
- the sampler is printed to the screen. If TRUE, the iteration
+ the sampler is printed to the screen. If \code{verbose} is greater
+ than 0 the iteration
number, the beta vector, and the Metropolis-Hastings acceptance rate
- are printed to the screen every 500 iterations.}
+ are printed to the screen every \code{verbose}th iteration.}
\item{seed}{The seed for the random number generator. If NA, the Mersenne
Twister generator is used with default seed 12345; if an integer is
diff --git a/man/MCMCordfactanal.Rd b/man/MCMCordfactanal.Rd
index 67e9cb8..a4f3b9f 100644
--- a/man/MCMCordfactanal.Rd
+++ b/man/MCMCordfactanal.Rd
@@ -14,7 +14,7 @@
\usage{
MCMCordfactanal(x, factors, lambda.constraints=list(),
data=parent.frame(), burnin = 1000, mcmc = 20000,
- thin=1, tune=NA, verbose = FALSE, seed = NA,
+ thin=1, tune=NA, verbose = 0, seed = NA,
lambda.start = NA, l0=0, L0=0,
store.lambda=TRUE, store.scores=FALSE,
drop.constantvars=TRUE, ... )
@@ -56,8 +56,10 @@ MCMCordfactanal(x, factors, lambda.constraints=list(),
strictly positive.}
\item{verbose}{A switch which determines whether or not the progress of
- the sampler is printed to the screen. If TRUE, the iteration number and
- the Metropolis-Hastings acceptance rate are printed to the screen.}
+ the sampler is printed to the screen. If \code{verbose} is greater
+ than 0 the iteration number and
+ the Metropolis-Hastings acceptance rate are printed to the screen
+ every \code{verbose}th iteration.}
\item{seed}{The seed for the random number generator. If NA, the Mersenne
Twister generator is used with default seed 12345; if an integer is
@@ -214,7 +216,7 @@ MCMCordfactanal(x, factors, lambda.constraints=list(),
posterior <- MCMCordfactanal(~Composition+Drawing+Colour+Expression,
data=new.painters, factors=1,
lambda.constraints=list(Drawing=list(2,"+")),
- burnin=5000, mcmc=500000, thin=200, verbose=TRUE,
+ burnin=5000, mcmc=500000, thin=200, verbose=500,
L0=0.5, store.lambda=TRUE,
store.scores=TRUE, tune=1.2)
plot(posterior)
diff --git a/man/MCMCpanel.Rd b/man/MCMCpanel.Rd
index 47781fa..11ad638 100644
--- a/man/MCMCpanel.Rd
+++ b/man/MCMCpanel.Rd
@@ -15,7 +15,7 @@
\usage{
MCMCpanel(obs, Y, X, W, burnin = 1000, mcmc = 10000, thin = 5,
- verbose = FALSE, seed = NA, sigma2.start = NA,
+ verbose = 0, seed = NA, sigma2.start = NA,
D.start = NA, b0 = 0, B0 = 1, eta0, R0, nu0 = 0.001,
delta0 = 0.001, ...)
}
@@ -52,8 +52,9 @@ MCMCpanel(obs, Y, X, W, burnin = 1000, mcmc = 10000, thin = 5,
\item{verbose}{A switch which determines whether or not the progress of
- the sampler is printed to the screen. If TRUE, the iteration number
- and parameters are printed to the screen.}
+ the sampler is printed to the screen. If \code{verbose} is greater
+ than 0 the iteration number
+ and parameters are printed to the screen every \code{verbose}th iteration.}
\item{sigma2.start}{The starting value for the conditional error
variance. Default value of NA uses the least squares estimates.}
diff --git a/man/MCMCpoisson.Rd b/man/MCMCpoisson.Rd
index dc0a2fe..2b7574d 100644
--- a/man/MCMCpoisson.Rd
+++ b/man/MCMCpoisson.Rd
@@ -12,7 +12,7 @@
\usage{
MCMCpoisson(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
- thin = 1, tune = 1.1, verbose = FALSE, seed = NA, beta.start = NA,
+ thin = 1, tune = 1.1, verbose = 0, seed = NA, beta.start = NA,
b0 = 0, B0 = 0, ...) }
\arguments{
@@ -35,9 +35,10 @@ MCMCpoisson(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
\item{verbose}{A switch which determines whether or not the progress of
- the sampler is printed to the screen. If TRUE, the iteration number,
+ the sampler is printed to the screen. If \code{verbose} is greater
+ than 0 the iteration number,
the current beta vector, and the Metropolis acceptance rate are
- printed to the screen every 500 iterations.}
+ printed to the screen every \code{verbose}th iteration.}
\item{seed}{The seed for the random number generator. If NA, the Mersenne
Twister generator is used with default seed 12345; if an integer is
diff --git a/man/MCMCprobit.Rd b/man/MCMCprobit.Rd
index 2ed7219..4c64805 100644
--- a/man/MCMCprobit.Rd
+++ b/man/MCMCprobit.Rd
@@ -12,7 +12,7 @@
\usage{
MCMCprobit(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
- thin = 1, verbose = FALSE, seed = NA, beta.start = NA,
+ thin = 1, verbose = 0, seed = NA, beta.start = NA,
b0 = 0, B0 = 0, bayes.resid = FALSE, ...) }
\arguments{
@@ -28,8 +28,9 @@ MCMCprobit(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
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 TRUE, the iteration number and
- the betas are printed to the screen.}
+ the sampler is printed to the screen. If \code{verbose} is greater
+ than 0 the iteration number and
+ the betas are printed to the screen every \code{verbose}th iteration.}
\item{seed}{The seed for the random number generator. If NA, the Mersenne
Twister generator is used with default seed 12345; if an integer is
diff --git a/man/MCMCregress.Rd b/man/MCMCregress.Rd
index 555577e..ba9307a 100644
--- a/man/MCMCregress.Rd
+++ b/man/MCMCregress.Rd
@@ -14,7 +14,7 @@
\usage{
MCMCregress(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
- thin = 1, verbose = FALSE, seed = NA, beta.start = NA,
+ thin = 1, verbose = 0, seed = NA, beta.start = NA,
b0 = 0, B0 = 0, c0 = 0.001, d0 = 0.001, ...) }
\arguments{
@@ -30,10 +30,10 @@ MCMCregress(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
MCMC iterations must be divisible by this value.}
\item{verbose}{A switch which determines whether or not the progress of
- the sampler is printed to the screen. If TRUE, the iteration number, the
- \eqn{\beta}{beta} vector, and the conditional error variance is printed to
- the screen
- every 500 iterations.}
+ the sampler is printed to the screen. If \code{verbose} is greater
+ than 0 the iteration number, the
+ \eqn{\beta}{beta} vector, and the error variance are printed to
+ the screen every \code{verbose}th iteration.}
\item{seed}{The seed for the random number generator. If NA, the Mersenne
Twister generator is used with default seed 12345; if an integer is
@@ -123,7 +123,7 @@ MCMCregress(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
\examples{
\dontrun{
line <- list(X = c(-2,-1,0,1,2), Y = c(1,3,3,3,5))
-posterior <- MCMCregress(Y~X, data=line, verbose=TRUE)
+posterior <- MCMCregress(Y~X, data=line, verbose=1000)
plot(posterior)
raftery.diag(posterior)
summary(posterior)
diff --git a/man/MCMCtobit.Rd b/man/MCMCtobit.Rd
index 100d9dd..80872fa 100644
--- a/man/MCMCtobit.Rd
+++ b/man/MCMCtobit.Rd
@@ -14,7 +14,7 @@
}
\usage{
MCMCtobit(formula, data = parent.frame(), below = 0, above = Inf,
- burnin = 1000, mcmc = 10000, thin = 1, verbose = FALSE, seed = NA,
+ burnin = 1000, mcmc = 10000, thin = 1, verbose = 0, seed = NA,
beta.start = NA, b0 = 0, B0 = 0, c0 = 0.001, d0 = 0.001, ...)
}
\arguments{
@@ -38,9 +38,10 @@ MCMCtobit(formula, data = parent.frame(), below = 0, above = Inf,
MCMC iterations must be divisible by this value.}
\item{verbose}{A switch which determines whether or not the progress of
- the sampler is printed to the screen. If TRUE, the iteration number, the
- \eqn{\beta}{beta} vector, and the conditional error variance is printed to
- the screen every 1000 iterations.}
+ the sampler is printed to the screen. If \code{verbose} is greater
+ than 0 the iteration number, the
+ \eqn{\beta}{beta} vector, and the error variance is printed to
+ the screen every \code{verbose}th iteration.}
\item{seed}{The seed for the random number generator. If NA, the Mersenne
Twister generator is used with default seed 12345; if an integer is
@@ -146,7 +147,8 @@ MCMCtobit(formula, data = parent.frame(), below = 0, above = Inf,
library(survival)
example(tobin)
summary(tfit)
-tfit.mcmc <- MCMCtobit(durable ~ age + quant, data=tobin, mcmc=30000, verbose=TRUE)
+tfit.mcmc <- MCMCtobit(durable ~ age + quant, data=tobin, mcmc=30000,
+ verbose=1000)
plot(tfit.mcmc)
raftery.diag(tfit.mcmc)
summary(tfit.mcmc)
diff --git a/man/invgamma.Rd b/man/invgamma.Rd
index 17667d5..466d53c 100644
--- a/man/invgamma.Rd
+++ b/man/invgamma.Rd
@@ -4,19 +4,19 @@
\alias{InvGamma}
\title{The Inverse Gamma Distribution}
\description{
- Density function and random generation from the inverse Gamma distribution.
+ Density function and random generation from the inverse gamma distribution.
}
\usage{
-rinvgamma(n, shape, rate = 1)
-dinvgamma(x, shape, rate = 1)
+rinvgamma(n, shape, scale = 1)
+dinvgamma(x, shape, scale = 1)
}
\arguments{
\item{x}{Scalar location to evaluate density.}
\item{n}{Number of draws from the distribution.}
\item{shape}{Scalar shape parameter.}
- \item{rate}{Scalar rate parameter (default value one).}
+ \item{scale}{Scalar scale parameter (default value one).}
}
\value{
@@ -25,6 +25,13 @@ dinvgamma(x, shape, rate = 1)
consistent with the Gamma Distribution in the stats package.
}
+\details{
+An inverse gamma random variable with shape \eqn{a}{a} and scale
+\eqn{b}{b} has mean \eqn{\frac{b}{a-1}}{b/(a-1)} (assuming \eqn{a>1}{a>1}) and
+variance \eqn{\frac{b^2}{(a-1)^2(a-2)}}{(b^2)/((a-1)^2 (a-2))} (assuming
+\eqn{a>2}{a>2}).
+}
+
\examples{
density <- dinvgamma(4.2, 1.1)
draws <- rinvgamma(10, 3.2)
diff --git a/man/iwishart.Rd b/man/iwishart.Rd
index f31b163..9992988 100644
--- a/man/iwishart.Rd
+++ b/man/iwishart.Rd
@@ -14,14 +14,20 @@
\arguments{
\item{W}{Positive definite matrix W \eqn{(p \times p)}{(p x p)}.}
- \item{v}{Inverse Wishart degrees of freedom (scalar).}
- \item{S}{Inverse Wishart scale matrix \eqn{(p \times p)}{(p x p)}.}}
+ \item{v}{Degrees of freedom (scalar).}
+ \item{S}{Scale matrix \eqn{(p \times p)}{(p x p)}.}}
\value{
\code{diwish} evaluates the density at positive definite matrix W.
\code{riwish} generates one random draw from the distribution.
}
+\details{
+ The mean of an inverse Wishart random variable with \code{v} degrees
+ of freedom and scale matrix \code{S} is \eqn{(v-p-1)^{-1}S}{1/(v-p-1)
+ S}.
+}
+
\examples{
density <- diwish(matrix(c(2,-.3,-.3,4),2,2), 3, matrix(c(1,.3,.3,1),2,2))
draw <- riwish(3, matrix(c(1,.3,.3,1),2,2))
diff --git a/man/wishart.Rd b/man/wishart.Rd
index f931d90..13e0ac6 100644
--- a/man/wishart.Rd
+++ b/man/wishart.Rd
@@ -14,8 +14,8 @@
\arguments{
\item{W}{Positive definite matrix W \eqn{(p \times p)}{(p x p)}.}
- \item{v}{Wishart degrees of freedom (scalar).}
- \item{S}{Wishart scale matrix \eqn{(p \times p)}{(p x p)}.}
+ \item{v}{Degrees of freedom (scalar).}
+ \item{S}{Inverse scale matrix \eqn{(p \times p)}{(p x p)}.}
}
\value{
@@ -23,6 +23,10 @@
\code{rwish} generates one random draw from the distribution.
}
+\details{
+The mean of a Wishart random variable with \code{v} degrees of freedom
+and inverse scale matrix \code{S} is \eqn{vS}{vS}.
+}
\examples{
density <- dwish(matrix(c(2,-.3,-.3,4),2,2), 3, matrix(c(1,.3,.3,1),2,2))
diff --git a/src/MCMCdynamicEI.cc b/src/MCMCdynamicEI.cc
index 46b1793..25e7bd1 100644
--- a/src/MCMCdynamicEI.cc
+++ b/src/MCMCdynamicEI.cc
@@ -323,6 +323,9 @@ extern "C"{
delta2 = (delta1 + SSE[0])*0.5;
sigma_theta1 = stream->rigamma(nu2, delta2);
}
+
+ // allow user interrupts
+ R_CheckUserInterrupt();
// @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
@@ -394,13 +397,13 @@ extern "C"{
// print output to screen
- if (verbose==1 && (iter%1000)==0){
+ if (verbose>0 && (iter%verbose)==0){
Rprintf("\nMCMCdynamicEI iteration %i of %i \n", (iter+1),
tot_iter);
}
-
+
// allow user interrupts
- void R_CheckUserInterrupt(void);
+ R_CheckUserInterrupt();
}
delete stream; // clean up random number stream
diff --git a/src/MCMCfactanal.cc b/src/MCMCfactanal.cc
index 5e4c200..d6a63ee 100644
--- a/src/MCMCfactanal.cc
+++ b/src/MCMCfactanal.cc
@@ -124,7 +124,7 @@ extern "C" {
// print results to screen
- if (iter % 500 == 0 && verbose[0] == 1){
+ if (iter % verbose[0] == 0 && verbose[0] > 0){
Rprintf("\n\nMCMCfactanal iteration %i of %i \n", (iter+1), tot_iter);
Rprintf("Lambda = \n");
for (int i=0; i<K; ++i){
@@ -157,9 +157,10 @@ extern "C" {
}
count++;
}
-
+
// allow user interrupts
- void R_CheckUserInterrupt(void);
+ R_CheckUserInterrupt();
+
} // end Gibbs loop
delete stream; // clean up random number stream
diff --git a/src/MCMChierEI.cc b/src/MCMChierEI.cc
index 2022d7b..f77d381 100644
--- a/src/MCMChierEI.cc
+++ b/src/MCMChierEI.cc
@@ -323,6 +323,9 @@ extern "C"{
nu2 = (nu1 + ntables)*0.5;
delta2 = (delta1 + SSE[0])*0.5;
sigma1 = stream->rigamma(nu2,delta2);
+
+ // allow user interrupts
+ R_CheckUserInterrupt();
}
// @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
@@ -413,13 +416,13 @@ extern "C"{
}
// print output to screen
- if (verbose==1 && (iter%1000)==0){
+ if (verbose>0 && (iter%verbose)==0){
Rprintf("\nMCMChierEI iteration %i of %i \n", (iter+1),
tot_iter);
}
// allow user interrupts
- void R_CheckUserInterrupt(void);
+ R_CheckUserInterrupt();
}
delete stream; // clean up random number stream
diff --git a/src/MCMCirt1d.cc b/src/MCMCirt1d.cc
index 5d139e2..c1848bc 100644
--- a/src/MCMCirt1d.cc
+++ b/src/MCMCirt1d.cc
@@ -17,6 +17,9 @@
#include <R.h> // needed to use Rprintf()
#include <R_ext/Utils.h> // needed to allow user interrupts
+using namespace SCYTHE;
+using namespace std;
+
extern "C"{
// function called by R to fit model
@@ -112,7 +115,7 @@ extern "C"{
theta_ineq, stream);
// print results to screen
- if (*verbose == 1 && iter % 100 == 0){
+ if (*verbose > 0 && iter % *verbose == 0){
Rprintf("\n\nMCMCirt1d iteration %i of %i \n",
(iter+1), tot_iter);
//Rprintf("theta = \n");
@@ -132,8 +135,9 @@ extern "C"{
eta_store(count, l) = eta[l];
count++;
}
+
+ R_CheckUserInterrupt(); // allow user interrupts
- void R_CheckUserInterrupt(void); // allow user interrupts
} // end Gibbs loop
delete stream; // clean up random number stream
diff --git a/src/MCMClogit.cc b/src/MCMClogit.cc
index 65b1b7c..6104744 100644
--- a/src/MCMClogit.cc
+++ b/src/MCMClogit.cc
@@ -117,13 +117,14 @@ extern "C"{
// store values in matrices
if (iter >= *burnin && ((iter % *thin)==0)){
+
for (int j = 0; j < k; j++)
storemat(count, j) = beta[j];
++count;
}
// print output to stdout
- if(*verbose == 1 && iter % 500 == 0){
+ if(*verbose > 0 && iter % *verbose == 0){
Rprintf("\n\nMCMClogit iteration %i of %i \n", (iter+1), tot_iter);
Rprintf("beta = \n");
for (int j=0; j<k; ++j)
@@ -133,7 +134,8 @@ extern "C"{
static_cast<double>(iter+1));
}
- void R_CheckUserInterrupt(void); // allow user interrupts
+ R_CheckUserInterrupt(); // allow user interrupts
+
}// end MCMC loop
delete stream; // clean up random number stream
diff --git a/src/MCMCmetrop1R.cc b/src/MCMCmetrop1R.cc
index c781019..07ff7a3 100644
--- a/src/MCMCmetrop1R.cc
+++ b/src/MCMCmetrop1R.cc
@@ -137,7 +137,7 @@ extern "C" {
++count;
}
- if (iter % 500 == 0 && INTEGER(verbose)[0] != 0 ) {
+ if (iter % INTEGER(verbose)[0] == 0 && INTEGER(verbose)[0] > 0 ) {
Rprintf("MCMCmetrop1R iteration %i of %i \n", (iter+1), tot_iter);
Rprintf("theta = \n");
for (int i=0; i<npar; ++i)
@@ -148,8 +148,9 @@ extern "C" {
static_cast<double>(iter+1));
}
+
UNPROTECT(1);
- void R_CheckUserInterrupt(void); // allow user interrupts
+ R_CheckUserInterrupt(); // allow user interrupts
}
// put the sample into a SEXP and return it
diff --git a/src/MCMCmixfactanal.cc b/src/MCMCmixfactanal.cc
index dcb7136..7e34d46 100644
--- a/src/MCMCmixfactanal.cc
+++ b/src/MCMCmixfactanal.cc
@@ -271,7 +271,7 @@ mixfactanalpost (double* sampledata, const int* samplerow,
}
// print results to screen
- if (*verbose == 1 && iter % 500 == 0){
+ if (*verbose > 0 && iter % *verbose == 0){
Rprintf("\n\nMCMCmixfactanal iteration %i of %i \n", (iter+1), tot_iter);
Rprintf("Lambda = \n");
for (int i=0; i<K; ++i){
@@ -322,7 +322,8 @@ mixfactanalpost (double* sampledata, const int* samplerow,
}
// allow user interrupts
- void R_CheckUserInterrupt(void);
+ R_CheckUserInterrupt();
+
} // end Gibbs loop
delete stream; // clean up random number stream
diff --git a/src/MCMCmnlMH.cc b/src/MCMCmnlMH.cc
index 7def301..230d482 100644
--- a/src/MCMCmnlMH.cc
+++ b/src/MCMCmnlMH.cc
@@ -141,14 +141,16 @@ extern "C" {
// store values in matrices
if (iter >= *burnin && ((iter % *thin)==0)){
+
for (int j = 0; j < k; j++)
storemat(count, j) = beta[j];
++count;
}
// print output to stdout
- if(*verbose == 1 && iter % 500 == 0){
- Rprintf("\n\nMCMCmnl Metropolis iteration %i of %i \n", (iter+1), tot_iter);
+ if(*verbose > 0 && iter % *verbose == 0){
+ Rprintf("\n\nMCMCmnl Metropolis iteration %i of %i \n",
+ (iter+1), tot_iter);
Rprintf("beta = \n");
for (int j=0; j<k; ++j)
Rprintf("%10.5f\n", beta[j]);
@@ -157,7 +159,7 @@ extern "C" {
static_cast<double>(iter+1));
}
- void R_CheckUserInterrupt(void); // allow user interrupts
+ R_CheckUserInterrupt(); // allow user interrupts
} // end MCMC loop
delete stream; // clean up random number stream
diff --git a/src/MCMCmnlslice.cc b/src/MCMCmnlslice.cc
index 8b97e3a..3a34cdb 100644
--- a/src/MCMCmnlslice.cc
+++ b/src/MCMCmnlslice.cc
@@ -313,7 +313,7 @@ extern "C" {
}
// print output to stdout
- if(*verbose == 1 && iter % 50 == 0){
+ if(*verbose > 0 && iter % *verbose == 0){
Rprintf("\n\nMCMCmnl slice iteration %i of %i \n", (iter+1),
tot_iter);
Rprintf("beta = \n");
@@ -321,7 +321,8 @@ extern "C" {
Rprintf("%10.5f\n", beta[j]);
}
- void R_CheckUserInterrupt(void); // allow user interrupts
+ R_CheckUserInterrupt(); // allow user interrupts
+
} // end MCMC loop
delete stream; // clean up random number stream
diff --git a/src/MCMCoprobit.cc b/src/MCMCoprobit.cc
index 31e82a2..addb3b6 100644
--- a/src/MCMCoprobit.cc
+++ b/src/MCMCoprobit.cc
@@ -163,7 +163,7 @@ extern "C"{
}
// print output to stdout
- if(*verbose == 1 && iter % 500 == 0){
+ if(*verbose > 0 && iter % *verbose == 0){
Rprintf("\n\nMCMCoprobit iteration %i of %i \n", (iter+1), tot_iter);
Rprintf("beta = \n");
for (int j=0; j<k; ++j)
@@ -173,8 +173,7 @@ extern "C"{
static_cast<double>(iter+1));
}
-
- void R_CheckUserInterrupt(void); // allow user interrupts
+ R_CheckUserInterrupt(); // allow user interrupts
}
delete stream; // clean up random number stream
diff --git a/src/MCMCordfactanal.cc b/src/MCMCordfactanal.cc
index f531143..02a1fb9 100644
--- a/src/MCMCordfactanal.cc
+++ b/src/MCMCordfactanal.cc
@@ -224,7 +224,7 @@ ordfactanalpost (double* sampledata, const int* samplerow,
// print results to screen
- if (verbose[0] == 1 && iter % 500 == 0 && *outswitch == 1){
+ if (verbose[0] > 0 && iter % verbose[0] == 0 && *outswitch == 1){
Rprintf("\n\nMCMCordfactanal iteration %i of %i \n", (iter+1), tot_iter);
Rprintf("Lambda = \n");
for (int i=0; i<K; ++i){
@@ -240,16 +240,20 @@ ordfactanalpost (double* sampledata, const int* samplerow,
static_cast<double>((iter+1)));
}
}
- if (verbose[0] == 1 && iter % 500 == 0 && *outswitch == 2){
+ if (verbose[0] > 0 && iter % verbose[0] == 0 && *outswitch == 2){
Rprintf("\n\nMCMCirtKd iteration %i of %i \n", (iter+1), tot_iter);
}
// store results
if ((iter >= burnin[0]) && ((iter % thin[0]==0))) {
-
// store Lambda
if (storelambda[0]==1){
+ if (*outswitch==2){
+ for(int l=0; l<K; ++l){
+ Lambda(l,0) = Lambda(l,0) * -1.0;
+ }
+ }
Matrix<double> Lambda_store_vec = reshape(Lambda,1,K*D);
for (int l=0; l<K*D; ++l)
Lambda_store(count, l) = Lambda_store_vec[l];
@@ -268,9 +272,10 @@ ordfactanalpost (double* sampledata, const int* samplerow,
}
count++;
}
-
+
// allow user interrupts
- void R_CheckUserInterrupt(void);
+ R_CheckUserInterrupt();
+
} // end MCMC loop
delete stream; // clean up random number stream
diff --git a/src/MCMCpanel.cc b/src/MCMCpanel.cc
index 23fa7d9..e9d6d31 100644
--- a/src/MCMCpanel.cc
+++ b/src/MCMCpanel.cc
@@ -198,7 +198,7 @@ panelpost (double* sample, const int* samrow, const int* samcol,
}
// print output to stdout
- if(*verbose == 1 && g % 500 == 0) {
+ if(*verbose > 0 && g % *verbose == 0) {
Rprintf("\n\nMCMCpanel iteration %i of %i \n",
(g+1), Mtotiter);
Rprintf("beta = \n");
@@ -207,13 +207,11 @@ panelpost (double* sample, const int* samrow, const int* samcol,
Rprintf("sigma2 = %10.5f\n", sigma2_sim);
}
-
- void R_CheckUserInterrupt(void); // allow user interrupts
-
- }
-
- delete stream; // clean up random number stream
-
+ R_CheckUserInterrupt(); // allow user interrupt
+ }
+
+ delete stream; // clean up random number stream
+
// return posterior denisty sample to R
Matrix<double> storeagem = cbind(beta_holder, D_holder);
storeagem = cbind(storeagem, sigma2_holder);
diff --git a/src/MCMCpoisson.cc b/src/MCMCpoisson.cc
index ce516f6..4c5d26e 100644
--- a/src/MCMCpoisson.cc
+++ b/src/MCMCpoisson.cc
@@ -118,14 +118,14 @@ extern "C"{
}
// store values in matrices
- if (iter >= *burnin && ((iter % *thin)==0)){
+ if (iter >= *burnin && ((iter % *thin)==0)){
for (int j = 0; j < k; j++)
storemat(count, j) = beta[j];
++count;
}
// print output to stdout
- if(*verbose == 1 && iter % 500 == 0){
+ if(*verbose > 0 && iter % *verbose == 0){
Rprintf("\n\nMCMCpoisson iteration %i of %i \n", (iter+1), tot_iter);
Rprintf("beta = \n");
for (int j=0; j<k; ++j)
@@ -135,7 +135,8 @@ extern "C"{
static_cast<double>(iter+1));
}
- void R_CheckUserInterrupt(void); // allow user interrupts
+ R_CheckUserInterrupt(); // allow user interrupts
+
}// end MCMC loop
delete stream; // clean up random number stream
diff --git a/src/MCMCprobit.cc b/src/MCMCprobit.cc
index 45e70bb..3429253 100644
--- a/src/MCMCprobit.cc
+++ b/src/MCMCprobit.cc
@@ -98,14 +98,15 @@ extern "C"{
}
// print output to stdout
- if(*verbose == 1 && iter % 500 == 0){
+ if(*verbose > 0 && iter % *verbose == 0){
Rprintf("\n\nMCMCprobit iteration %i of %i \n", (iter+1), tot_iter);
Rprintf("beta = \n");
for (int j=0; j<k; ++j)
Rprintf("%10.5f\n", beta[j]);
}
- void R_CheckUserInterrupt(void); // allow user interrupts
+ R_CheckUserInterrupt(); // allow user interrupts
+
} // end MCMC loop
delete stream; // clean up random number stream
diff --git a/src/MCMCprobitres.cc b/src/MCMCprobitres.cc
index 2c1168c..072b664 100644
--- a/src/MCMCprobitres.cc
+++ b/src/MCMCprobitres.cc
@@ -107,14 +107,15 @@ extern "C"{
}
// print output to stdout
- if(*verbose == 1 && iter % 500 == 0){
+ if(*verbose > 0 && iter % *verbose == 0){
Rprintf("\n\nMCMCprobit iteration %i of %i \n", (iter+1), tot_iter);
Rprintf("beta = \n");
for (int j=0; j<k; ++j)
Rprintf("%10.5f\n", beta[j]);
}
- void R_CheckUserInterrupt(void); // allow user interrupts
+ R_CheckUserInterrupt(); // allow user interrupts
+
} // end MCMC loop
delete stream; // clean up random number stream
diff --git a/src/MCMCregress.cc b/src/MCMCregress.cc
index 5d732f6..1c277fa 100644
--- a/src/MCMCregress.cc
+++ b/src/MCMCregress.cc
@@ -101,7 +101,7 @@ extern "C" {
}
// print output to stdout
- if(*verbose == 1 && iter % 500 == 0) {
+ if(*verbose > 0 && iter % *verbose == 0) {
Rprintf("\n\nMCMCregress iteration %i of %i \n",
(iter+1), tot_iter);
Rprintf("beta = \n");
@@ -110,7 +110,8 @@ extern "C" {
Rprintf("sigma2 = %10.5f\n", sigma2);
}
- void R_CheckUserInterrupt(void); // allow user interrupts
+ R_CheckUserInterrupt(); // allow user interrupts
+
} // end MCMC loop
delete stream; // clean up random number stream
diff --git a/src/MCMCtobit.cc b/src/MCMCtobit.cc
index 4fca55b..4f1d2ad 100644
--- a/src/MCMCtobit.cc
+++ b/src/MCMCtobit.cc
@@ -112,6 +112,7 @@ extern "C" {
// store draws in storage matrix (or matrices)
if (iter >= *burnin && (iter % 1 == 0)) {
+
sigmamatrix (0, count) = sigma2;
for (int j = 0; j < k; j++)
betamatrix (j, count) = beta[j];
@@ -119,7 +120,7 @@ extern "C" {
}
// print output to stdout
- if(*verbose == 1 && iter % 1000 == 0) {
+ if(*verbose > 0 && iter % *verbose == 0) {
Rprintf("\n\nMCMCtobit iteration %i of %i \n",
(iter+1), tot_iter);
Rprintf("beta = \n");
@@ -127,8 +128,9 @@ extern "C" {
Rprintf("%10.5f\n", beta[j]);
Rprintf("sigma2 = %10.5f\n", sigma2);
}
+
+ R_CheckUserInterrupt(); // allow user interrupts
- void R_CheckUserInterrupt(void); // allow user interrupts
} // end MCMC loop
delete stream; // clean up random number stream
diff --git a/src/Makevars b/src/Makevars
index ae8cde5..3bff6ac 100644
--- a/src/Makevars
+++ b/src/Makevars
@@ -1,2 +1,2 @@
-PKG_CXXFLAGS = -DSCYTHE_COMPILE_DIRECT -DSCYTHE_NO_RANGE
+PKG_CXXFLAGS = -DSCYTHE_COMPILE_DIRECT -DSCYTHE_NO_RANGE -DHAVE_TRUNC
diff --git a/src/lecuyer.h b/src/lecuyer.h
index 806c8b0..f1cb79a 100644
--- a/src/lecuyer.h
+++ b/src/lecuyer.h
@@ -16,8 +16,11 @@
* Provides the class definition for the L'Ecuyer random number
* generator, a rng capable of generating many independent substreams.
* This class extends the abstract rng class by implementing runif().
- * Based on RngStream.cpp, which was released under the following
- * license:
+ * Based on RngStream.cpp.
+ *
+ * Pierre L'Ecuyer agreed to the following dual-licensing terms in an
+ * email received 7 August 2004. This dual-license was prompted by
+ * the Debian maintainers of R and MCMCpack.
*
* This software is Copyright (C) 2004 Pierre L'Ecuyer.
*
--
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