[r-cran-mcmcpack] 07/90: Imported Upstream version 0.6-5
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 a54539e2f7e9ee794c80176701dd379377cbd101
Author: Andreas Tille <tille at debian.org>
Date: Fri Dec 16 08:07:03 2016 +0100
Imported Upstream version 0.6-5
---
DESCRIPTION | 16 ++++++++--------
HISTORY | 23 +++++++++++++++++++++++
R/MCMCirt1d.R | 7 +++++--
R/MCMCmetrop1R.R | 35 +++++++++++++++++++++++++++++++++--
R/MCMCmixfactanal.R | 15 +++++++++++----
R/distn.R | 8 ++++++--
R/hidden.R | 3 ++-
R/zzz.R | 11 ++++++++++-
data/Senate.rda | Bin 294202 -> 568880 bytes
data/SupremeCourt.rda | Bin 6317 -> 3896 bytes
man/MCMCirt1d.Rd | 8 --------
man/MCMCirtKd.Rd | 10 ----------
man/MCMCmetrop1R.Rd | 21 ++++++++++++++++++++-
man/Senate.Rd | 2 +-
man/SupremeCourt.Rd | 4 ++--
src/MCMCfactanal.cc | 2 +-
src/MCMCmetrop1R.cc | 2 +-
src/MCMCpanel.cc | 4 ++--
18 files changed, 125 insertions(+), 46 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index 2446cdb..989de0e 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,20 +1,20 @@
Package: MCMCpack
-Version: 0.6-4
-Date: 2005-3-17
+Version: 0.6-5
+Date: 2005-6-29
Title: Markov chain Monte Carlo (MCMC) Package
Author: Andrew D. Martin <admartin at wustl.edu>, and
Kevin M. Quinn <kevin_quinn at harvard.edu>
Maintainer: Andrew D. Martin <admartin at wustl.edu>
Depends: R (>= 2.0.1), coda (>= 0.9-2), MASS
-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
- Library Version 1.0. All models return coda mcmc objects that
- can then be summarized using the coda package. MCMCpack
+Description: This package contains functions to perform Bayesian
+ inference using posterior simulation for a number of statistical
+ models. All simulation is done in compiled C++ written in the Scythe
+ Statistical Library Version 1.0. All models return coda mcmc objects
+ that can then be summarized using the coda package. MCMCpack
also contains some useful utility functions, including some
additional density functions and pseudo-random number generators
for statistical distributions, a general purpose Metropolis
sampling algorithm, and tools for visualization.
License: GPL version 2 or newer
URL: http://mcmcpack.wustl.edu
-Packaged: Thu Mar 17 09:42:26 2005; adm
+Packaged: Wed Jun 29 17:57:40 2005; adm
diff --git a/HISTORY b/HISTORY
index dca215b..271a231 100644
--- a/HISTORY
+++ b/HISTORY
@@ -2,6 +2,29 @@
// Changes and Bug Fixes
//
+0.6-4 to 0.6-5
+ * added a check so one can run MCMCirt1d without passing constraints
+ * fixed the Senate dataset so there are no duplicate names in the
+ member variables, and modified examples accordingly
+ * slightly edited DESCRIPTION file to get some important keywords
+ in the first sentence
+ * compute dinvgamma() on the log scale so it does not blow up
+ for large shape and scale parameters
+ * changed how MCMCmetrop1R handles a non-negative definite Hessian
+ from optim--sampling can now proceed if the user flips a switch
+ * fixed minor error in how std.mean was passed to MCMCmixfactanal
+ * fixed the error in the beta full conditional in MCMCpanel that resulted
+ from a typo in the Chib and Carlin paper
+ * changed factor.score.eigen.start() in hidden.R to use squared
+ distances rather than distances so as to be compatible with
+ Poole, 2005. _Spatial Models of Parliamentary Voting_
+ * fixed a minor bug in how the verbose argument was handled in the C++
+ code for MCMCmetrop1R and MCMCfactanal [thanks to Lingji Chen]
+ * fixed a minor documentation bug in the MCMCirtKd model [thanks to
+ Guillermo Rosas]
+ * fixed the text echoed at start-up so it does not have to be updated
+ year to year
+
0.6-3 to 0.6-4
* fixed bug with verbose in MCMCmetrop1R pointed out by a referee for
Rnews
diff --git a/R/MCMCirt1d.R b/R/MCMCirt1d.R
index 21b0d1e..972635e 100644
--- a/R/MCMCirt1d.R
+++ b/R/MCMCirt1d.R
@@ -35,8 +35,11 @@
subject.names <- rownames(datamatrix)
## setup constraints on theta
- for (i in 1:length(theta.constraints)){
- theta.constraints[[i]] <- list(as.integer(1), theta.constraints[[i]][1])
+ if(length(theta.constraints) != 0) {
+ for (i in 1:length(theta.constraints)){
+ theta.constraints[[i]] <-
+ list(as.integer(1), theta.constraints[[i]][1])
+ }
}
holder <- build.factor.constraints(theta.constraints, t(datamatrix), J, 1)
theta.eq.constraints <- holder[[1]]
diff --git a/R/MCMCmetrop1R.R b/R/MCMCmetrop1R.R
index 4bc4da8..9380859 100644
--- a/R/MCMCmetrop1R.R
+++ b/R/MCMCmetrop1R.R
@@ -2,12 +2,13 @@
## random walk Metropolis algorithm
##
## KQ 6/24/2004
+## modified to work with non-invertible Hessian KQ 6/28/2005
##
"MCMCmetrop1R" <- function(fun, theta.init,
burnin=500, mcmc=20000, thin=1,
tune=1, verbose=0, seed=NA, logfun=TRUE,
- optim.trace=0, optim.REPORT=10,
+ force.samp=FALSE, optim.trace=0, optim.REPORT=10,
optim.maxit=500, ...){
## error checking here
@@ -46,8 +47,38 @@
if(opt.out$convergence!=0){
warning("Mode and Hessian were not found with call to optim().\nSampling proceeded anyway. \n")
}
+
+
+ CC <- NULL
+ try(CC <- chol(-1*opt.out$hessian), silent=TRUE)
+ hess.new <- opt.out$hessian
+ hess.flag <- 0
+ if (force.samp==TRUE){
+ while (is.null(CC)){
+ hess.flag <- 1
+ hess.new <- hess.new - diag(diag(0.01 * abs(opt.out$hessian)))
+ try(CC <- chol(-1*hess.new), silent=TRUE)
+ }
+ }
+ else{
+ if (is.null(CC)){
+ hess.flag <- 2
+ }
+ }
+ if (hess.flag==1){
+ warning("Hessian from call to optim() not negative definite.\nSampling proceeded after enforcing negative definiteness. \n")
+ }
+ if (hess.flag==2){
+ cat("Hessian from call to optim() not negative definite.\n")
+ cat("Sampling (as specified) cannot proceed.\n")
+ stop("Check data and fun() and call MCMCmetrop1R() again. \n",
+ call.=FALSE)
+ }
- propvar <- tune %*% solve(-1*opt.out$hessian) %*% tune
+ propvar <- tune %*% solve(-1*hess.new) %*% tune
+ ## the old way that only worked when Hessian was neg. def.
+ ##propvar <- tune %*% solve(-1*opt.out$hessian) %*% tune
+
## Call the C++ function to do the MCMC sampling
sample <- .Call("MCMCmetrop1R_cc", fun, as.double(theta.init),
diff --git a/R/MCMCmixfactanal.R b/R/MCMCmixfactanal.R
index 2cb04dc..14b017f 100644
--- a/R/MCMCmixfactanal.R
+++ b/R/MCMCmixfactanal.R
@@ -22,6 +22,7 @@
##
## 12/2/2003
## Revised to accommodate new spec 7/20/2004
+## Minor bug fix regarding std.mean 6/25/2004
##
##########################################################################
@@ -43,7 +44,8 @@
mf$factors <- mf$lambda.constraints <- mf$burnin <- mf$mcmc <- NULL
mf$thin <- mf$tune <- mf$verbose <- mf$seed <- NULL
mf$lambda.start <- mf$l0 <- mf$L0 <- mf$a0 <- mf$b0 <- NULL
- mf$store.lambda <- mf$store.scores <- mf$std.var <- mf$... <- NULL
+ mf$store.lambda <- mf$store.scores <- mf$std.mean <- NULL
+ mf$std.var <- mf$... <- NULL
mf$drop.unused.levels <- TRUE
mf[[1]] <- as.name("model.frame")
mf$na.action <- 'na.pass'
@@ -61,8 +63,13 @@
ncat[i] <- -999
X[is.na(X[,i]), i] <- -999
}
- else if (is.ordered(X[, i])) {
ncat[i] <- nlevels(X[, i])
temp <- as.integer(X[,i])
temp <- ifelse(is.na(X[,i]) | (X[,i] == "<NA>"), -999, temp)
X[, i] <- temp
}
- else {
+ else if (is.ordered(X[, i])) {
+ ncat[i] <- nlevels(X[, i])
+ temp <- as.integer(X[,i])
+ temp <- ifelse(is.na(X[,i]) | (X[,i] == "<NA>"), -999, temp)
+ X[, i] <- temp
+ }
+ else {
stop("Manifest variable ", dimnames(X)[[2]][i],
" neither ordered factor nor numeric variable.\n")
}
@@ -87,7 +94,7 @@
if (std.var){
for (i in 1:K){
if (ncat[i] == -999){
- X[,i] <- (X[,i])/sd(X[,i])
+ X[,i] <- (X[,i] - mean(X[,i]))/sd(X[,i]) + mean(X[,i])
}
}
}
diff --git a/R/distn.R b/R/distn.R
index c69836d..bfa10f6 100644
--- a/R/distn.R
+++ b/R/distn.R
@@ -198,10 +198,14 @@
if(shape <= 0 | scale <=0) {
stop("Shape or scale parameter negative in dinvgamma().\n")
}
-
+
alpha <- shape
beta <- scale
- return(beta^alpha / gamma(alpha) * x^(-1*(alpha + 1)) * exp(-beta/x))
+
+ # done on log scale to allow for large alphas and betas
+ log.density <- alpha * log(beta) - lgamma(alpha) -
+ (alpha + 1) * log(x) - (beta/x)
+ return(exp(log.density))
}
## generate draws from the inverse gamma density (using
diff --git a/R/hidden.R b/R/hidden.R
index 53f4245..647507d 100644
--- a/R/hidden.R
+++ b/R/hidden.R
@@ -216,6 +216,7 @@
## takes a subject by subject agreement score matrix as input
"factor.score.eigen.start" <- function(A, factors){
+ A <- (1 - A)^2
AA <- A
arow <- matrix(NA, nrow(A), 1)
acol <- matrix(NA, ncol(A), 1)
@@ -231,7 +232,7 @@
for (i in 1:nrow(A)){
for (j in 1:ncol(A)){
- AA[i,j] <- (A[i,j]-arow[i]-acol[j]+matrixmean)/(2)
+ AA[i,j] <- (A[i,j]-arow[i]-acol[j]+matrixmean)/(-2)
}
}
diff --git a/R/zzz.R b/R/zzz.R
index 417fc49..7666a6c 100644
--- a/R/zzz.R
+++ b/R/zzz.R
@@ -1,6 +1,14 @@
.onAttach <- function(...) {
+
+ # figure out year automatically (probably could be done more elegantly)
+ date <- date()
+ x <- regexpr("[0-9]{4}", date)
+ this.year <- substr(date, x[1], x[1] + attr(x, "match.length") - 1)
+
+ # echo output to screen
cat("##\n## Markov Chain Monte Carlo Package (MCMCpack)\n")
- cat("## Copyright (C) 2003-2005, Andrew D. Martin and Kevin M. Quinn\n")
+ cat("## Copyright (C) 2003-", this.year,
+ " Andrew D. Martin and Kevin M. Quinn\n", sep="")
cat("##\n## Support provided by the U.S. National Science Foundation\n")
cat("## (Grants SES-0350646 and SES-0350613)\n##\n")
require(coda, quietly=TRUE)
@@ -10,3 +18,4 @@
.onUnload <- function(libpath) {
library.dynam.unload("MCMCpack", libpath)
}
+
diff --git a/data/Senate.rda b/data/Senate.rda
index 479c5e9..1a28cfe 100644
Binary files a/data/Senate.rda and b/data/Senate.rda differ
diff --git a/data/SupremeCourt.rda b/data/SupremeCourt.rda
index 107e7ee..1f1d1e7 100644
Binary files a/data/SupremeCourt.rda and b/data/SupremeCourt.rda differ
diff --git a/man/MCMCirt1d.Rd b/man/MCMCirt1d.Rd
index 7f0b222..3f1db55 100644
--- a/man/MCMCirt1d.Rd
+++ b/man/MCMCirt1d.Rd
@@ -200,15 +200,7 @@ 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]
-
posterior2 <- MCMCirt1d(Sen.rollcalls,
theta.constraints=list(KENNEDY=-2, HELMS=2),
burnin=2000, mcmc=100000, thin=20, verbose=500)
diff --git a/man/MCMCirtKd.Rd b/man/MCMCirtKd.Rd
index d218935..20de18c 100644
--- a/man/MCMCirtKd.Rd
+++ b/man/MCMCirtKd.Rd
@@ -220,17 +220,7 @@ MCMCirtKd(datamatrix, dimensions, item.constraints=list(),
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]
-
- # note that we need to transpose the data to get
- # the bills on the rows
posterior2 <- MCMCirtKd(Sen.rollcalls, dimensions=2,
burnin=5000, mcmc=50000, thin=10,
item.constraints=list(rc2=list(2,"-"), rc2=c(3,0),
diff --git a/man/MCMCmetrop1R.Rd b/man/MCMCmetrop1R.Rd
index 5b5f3c4..d1ceadd 100644
--- a/man/MCMCmetrop1R.Rd
+++ b/man/MCMCmetrop1R.Rd
@@ -9,7 +9,8 @@
\usage{
MCMCmetrop1R(fun, theta.init, burnin = 500, mcmc = 20000, thin = 1,
tune = 1, verbose = 0, seed=NA, logfun = TRUE,
- optim.trace = 0, optim.REPORT = 10, optim.maxit = 500, ...)
+ force.samp=FALSE, optim.trace = 0, optim.REPORT = 10,
+ optim.maxit = 500, ...)
}
\arguments{
\item{fun}{The (log)density from which to take a sample. This should
@@ -51,6 +52,24 @@ MCMCmetrop1R(fun, theta.init, burnin = 500, mcmc = 20000, thin = 1,
\item{logfun}{Logical indicating whether \code{fun} returns the natural log
of a density function (TRUE) or a density (FALSE).}
+ \item{force.samp}{Logical indicating whether the sampling should
+ proceed if the Hessian calculated from the initial call to
+ \code{optim} routine to maximize the (log)posterior
+ is not negative definite. If \code{force.samp==TRUE}
+ and the Hessian from \code{optim} is non-negative definite, the
+ Hessian is rescaled by subtracting small values from it's main
+ diagonal until it is negative definite. Sampling proceeds using
+ this rescaled Hessian in place of the original Hessian from
+ \code{optim}. By default, if \code{force.samp==FALSE} and the Hessian from
+ \code{optim} is non-negative definite, an error message is
+ printed and the call to \code{MCMCmetrop1R} is terminated.
+
+ \emph{Please note that a non-negative Hessian at the mode is a very
+ good indication that the model being fitted is not
+ identified. Thus, \code{force.samp} should only be set equal to
+ \code{TRUE} with great caution.}
+ }
+
\item{optim.trace}{The value of the \code{trace} parameter sent to
\code{optim} during an initial maximization of \code{fun}. }
diff --git a/man/Senate.Rd b/man/Senate.Rd
index 5b1ce9e..4c17dee 100644
--- a/man/Senate.Rd
+++ b/man/Senate.Rd
@@ -22,7 +22,7 @@ data(Senate)
}
\source{
- Keith Poole. 2003. \emph{106th Roll Call Vote Data}.
+ Keith Poole. 2005. \emph{106th Roll Call Vote Data}.
\url{http://voteview.uh.edu/}.
}
diff --git a/man/SupremeCourt.Rd b/man/SupremeCourt.Rd
index d5db85c..5520340 100644
--- a/man/SupremeCourt.Rd
+++ b/man/SupremeCourt.Rd
@@ -24,8 +24,8 @@ data(SupremeCourt)
}
\source{
- Harold J. Spaeth. 2003. \emph{Original United States Supreme Court Database:
- 1953-2001 Terms.} \url{http://polisci.msu.edu/pljp/}.
+ Harold J. Spaeth. 2005. \emph{Original United States Supreme Court Database:
+ 1953-2004 Terms.} \url{http://www.as.uky.edu/polisci/ulmerproject/sctdata.htm}.
}
\keyword{datasets}
diff --git a/src/MCMCfactanal.cc b/src/MCMCfactanal.cc
index d6a63ee..386018e 100644
--- a/src/MCMCfactanal.cc
+++ b/src/MCMCfactanal.cc
@@ -124,7 +124,7 @@ extern "C" {
// print results to screen
- if (iter % verbose[0] == 0 && verbose[0] > 0){
+ if (verbose[0] > 0 && iter % 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){
diff --git a/src/MCMCmetrop1R.cc b/src/MCMCmetrop1R.cc
index 07ff7a3..3db8cce 100644
--- a/src/MCMCmetrop1R.cc
+++ b/src/MCMCmetrop1R.cc
@@ -137,7 +137,7 @@ extern "C" {
++count;
}
- if (iter % INTEGER(verbose)[0] == 0 && INTEGER(verbose)[0] > 0 ) {
+ if (INTEGER(verbose)[0] > 0 && iter % 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)
diff --git a/src/MCMCpanel.cc b/src/MCMCpanel.cc
index e9d6d31..76fa6ed 100644
--- a/src/MCMCpanel.cc
+++ b/src/MCMCpanel.cc
@@ -135,8 +135,8 @@ panelpost (double* sample, const int* samrow, const int* samcol,
beta_var_sum = beta_var_sum + t(Xarr[i]) * invVi * Xarr[i];
beta_mean_sum = beta_mean_sum + t(Xarr[i]) * invVi * Yarr[i];
}
- Matrix<double> beta_sim_var = invpd(MB0 + (1/sigma2_sim) * beta_var_sum);
- Matrix<double> beta_sim_mean = beta_sim_var * (MB0 * Mb0 + (1/sigma2_sim) * beta_mean_sum);
+ Matrix<double> beta_sim_var = invpd(MB0 + beta_var_sum);
+ Matrix<double> beta_sim_mean = beta_sim_var * (MB0 * Mb0 + beta_mean_sum);
Matrix<double> beta_sim = beta_sim_mean + cholesky(beta_sim_var) *
stream->rnorm(Mp,1);
--
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