[r-cran-mcmcpack] 26/90: Imported Upstream version 0.9-2
Andreas Tille
tille at debian.org
Fri Dec 16 09:07:36 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 c04a2d6767a83beb1199267ea7f471429335b3e1
Author: Andreas Tille <tille at debian.org>
Date: Fri Dec 16 08:07:10 2016 +0100
Imported Upstream version 0.9-2
---
DESCRIPTION | 11 +-
HISTORY | 8 +
NAMESPACE | 4 +-
R/MCMCfactanal.R | 3 +-
R/MCMCirtKd.R | 74 ++--
R/MCMClogit.R | 2 +-
R/MCMCmixfactanal.R | 2 +
R/MCMCmnl.R | 111 ++++--
R/MCMCoprobit.R | 67 +++-
R/MCMCpoisson.R | 2 +-
R/MCMCpoissonChangepoint.R | 232 +++++-------
R/MCMCprobit.R | 1 +
R/MCMCregress.R | 3 +-
R/MCMCtobit.R | 131 +++----
R/btsutil.R | 297 +++++----------
R/zzz.R | 2 +-
README | 14 +-
man/MCMCdynamicEI.Rd | 4 +-
man/MCMChierEI.Rd | 4 +-
man/MCMCmnl.Rd | 13 +-
man/MCMCoprobit.Rd | 39 +-
man/MCMCpoissonChangepoint.Rd | 261 ++++++-------
man/MCMCtobit.Rd | 316 ++++++++--------
man/dtomog.Rd | 3 +-
man/{plotPostChangepoint.Rd => plotChangepoint.Rd} | 10 +-
man/{plotPostState.Rd => plotState.Rd} | 10 +-
man/tomog.Rd | 4 +-
src/MCMCmnl.h | 29 +-
src/MCMCmnlMH.cc | 277 +++++++++-----
src/MCMCoprobit.cc | 415 +++++++++++++++------
src/MCMCpoissonChangepoint.cc | 382 +++++++++++++++++++
src/Makevars | 3 +-
src/Makevars.in | 3 +-
33 files changed, 1675 insertions(+), 1062 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index 1aa6ce7..0a3ec6c 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,9 +1,10 @@
Package: MCMCpack
-Version: 0.9-1
-Date: 2007-08-22
+Version: 0.9-2
+Date: 2007-10-2
Title: Markov chain Monte Carlo (MCMC) Package
-Author: Andrew D. Martin <admartin at wustl.edu>, and
- Kevin M. Quinn <kevin_quinn at harvard.edu>
+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.5.0), coda (>= 0.11-3), MASS, stats
Description: This package contains functions to perform Bayesian
@@ -17,4 +18,4 @@ Description: This package contains functions to perform Bayesian
sampling algorithm, and tools for visualization.
License: GPL version 2
URL: http://mcmcpack.wustl.edu
-Packaged: Tue Aug 21 14:49:42 2007; adm
+Packaged: Tue Oct 2 19:42:13 2007; adm
diff --git a/HISTORY b/HISTORY
index 04b9ce7..b889b1a 100644
--- a/HISTORY
+++ b/HISTORY
@@ -2,6 +2,14 @@
// Changes and Bug Fixes
//
+0.9-1 to 0.9-2
+ * fixed Makevars per Ripley's email [PKG_CXXFLAGS to PKG_CPPFLAGS]
+ * fixed encoding on some documentation files [thanks to Kurt Hornik]
+ * added other estimation algorithms for MCMCoprobit
+ * added other estimation algorithms for MCMCmnl
+ * added MCMCpoissonChangepoint with estimation in compiled C++
+ * a variety of other minor fixes
+
0.8-2 to 0.9-1
* first release of JStatSoft version
* full port to Scythe 1.0.2 (nearly all C++ has been radically changed)
diff --git a/NAMESPACE b/NAMESPACE
index 135cfb5..67c7087 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -34,8 +34,8 @@ export(
MCMCregress,
MCMCSVDreg,
MCMCtobit,
- plotPostState,
- plotPostChangepoint,
+ plotState,
+ plotChangepoint,
PostProbMod,
procrustes,
rdirichlet,
diff --git a/R/MCMCfactanal.R b/R/MCMCfactanal.R
index 2c4e79c..27de0b5 100644
--- a/R/MCMCfactanal.R
+++ b/R/MCMCfactanal.R
@@ -110,7 +110,8 @@
else {
sample <- matrix(data=0, mcmc/thin, K*factors+K+factors*N)
}
-
+ posterior <- NULL
+
## call C++ code to do the sampling
auto.Scythe.call(output.object="posterior", cc.fun.name="MCMCfactanal",
sample.nonconst=sample, X=X, burnin=as.integer(burnin),
diff --git a/R/MCMCirtKd.R b/R/MCMCirtKd.R
index 855dc97..5869eef 100644
--- a/R/MCMCirtKd.R
+++ b/R/MCMCirtKd.R
@@ -1,37 +1,37 @@
-##########################################################################
-## sample from a K-dimensional two-parameter item response model with
-## probit link. This is just a wrapper function that calls
-## MCMCordfactanal.
-##
-## Andrew D. Martin
-## Washington University
-##
-## Kevin M. Quinn
-## Harvard University
-##
-## June 8, 2003
-##
-##########################################################################
-
-"MCMCirtKd" <-
- function(datamatrix, dimensions, item.constraints=list(),
- burnin = 1000, mcmc = 10000,
- thin=1, verbose = 0, seed = NA,
- alphabeta.start = NA, b0=0, B0=0,
- store.item=FALSE, store.ability=TRUE,
- drop.constant.items=TRUE, ... ) {
-
- datamatrix <- as.matrix(datamatrix)
-
- post <- MCMCordfactanal(x=datamatrix, factors=dimensions,
- lambda.constraints=item.constraints,
- burnin=burnin, mcmc=mcmc, thin=thin,
- tune=NA, verbose=verbose, seed=seed,
- lambda.start=alphabeta.start,
- l0=b0, L0=B0, store.lambda=store.item,
- store.scores=store.ability,
- drop.constantvars=drop.constant.items,
- model="MCMCirtKd")
- return(post)
- }
-
+##########################################################################
+## sample from a K-dimensional two-parameter item response model with
+## probit link. This is just a wrapper function that calls
+## MCMCordfactanal.
+##
+## Andrew D. Martin
+## Washington University
+##
+## Kevin M. Quinn
+## Harvard University
+##
+## June 8, 2003
+##
+##########################################################################
+
+"MCMCirtKd" <-
+ function(datamatrix, dimensions, item.constraints=list(),
+ burnin = 1000, mcmc = 10000,
+ thin=1, verbose = 0, seed = NA,
+ alphabeta.start = NA, b0=0, B0=0,
+ store.item=FALSE, store.ability=TRUE,
+ drop.constant.items=TRUE, ... ) {
+
+ datamatrix <- as.matrix(datamatrix)
+
+ post <- MCMCordfactanal(x=datamatrix, factors=dimensions,
+ lambda.constraints=item.constraints,
+ burnin=burnin, mcmc=mcmc, thin=thin,
+ tune=NA, verbose=verbose, seed=seed,
+ lambda.start=alphabeta.start,
+ l0=b0, L0=B0, store.lambda=store.item,
+ store.scores=store.ability,
+ drop.constantvars=drop.constant.items,
+ model="MCMCirtKd")
+ return(post)
+ }
+
diff --git a/R/MCMClogit.R b/R/MCMClogit.R
index 1d09c8d..84f3c25 100644
--- a/R/MCMClogit.R
+++ b/R/MCMClogit.R
@@ -116,7 +116,7 @@
propvar <- tune %*% V %*% tune
-
+ posterior <- NULL
## call C++ code to draw sample
if (is.null(user.prior.density)){
## define holder for posterior density sample
diff --git a/R/MCMCmixfactanal.R b/R/MCMCmixfactanal.R
index 35525a3..25f662c 100644
--- a/R/MCMCmixfactanal.R
+++ b/R/MCMCmixfactanal.R
@@ -36,6 +36,7 @@
std.mean=TRUE, std.var=TRUE, ... ) {
call <- match.call()
+ echo.name <- NULL
mt <- terms(x, data=data)
if (attr(mt, "response") > 0)
stop("Response not allowed in formula in MCMCmixfactanal().\n")
@@ -261,6 +262,7 @@
accepts <- matrix(0, K, 1)
# Call the C++ code to do the real work
+ posterior <- NULL
posterior <- .C("mixfactanalpost",
samdata = as.double(sample),
samrow = as.integer(nrow(sample)),
diff --git a/R/MCMCmnl.R b/R/MCMCmnl.R
index 1d41b06..316f286 100644
--- a/R/MCMCmnl.R
+++ b/R/MCMCmnl.R
@@ -167,20 +167,42 @@
}
+## MNL log-posterior function (used to get starting values)
+## vector Y without NAs
+"mnl.logpost.noNA" <- function(beta, new.Y, X, b0, B0){
+ nobs <- length(new.Y)
+ ncat <- nrow(X) / nobs
+ Xb <- X %*% beta
+ Xb <- matrix(Xb, byrow=TRUE, ncol=ncat)
+ indices <- cbind(1:nobs, new.Y)
+ Xb.reform <- Xb[indices]
+ eXb <- exp(Xb)
+ #denom <- log(apply(eXb, 1, sum))
+ z <- rep(1, ncat)
+ denom <- log(eXb %*% z)
-## MNL log-likelihood function (used to get starting values)
-"mnl.loglike" <- function(beta, Y, X){
- k <- ncol(X)
- numer <- exp(X %*% beta)
- numer[Y== -999] <- NA
- numer.mat <- matrix(numer, nrow(Y), ncol(Y), byrow=TRUE)
- denom <- apply(numer.mat, 1, sum, na.rm=TRUE)
- choice.probs <- numer.mat / denom
- Yna <- Y
- Yna[Y == -999] <- NA
- log.like.mat <- log(choice.probs) * Yna
- log.like <- sum(apply(log.like.mat, 1, sum, na.rm=TRUE))
- return(log.like)
+ log.prior <- 0.5 * t(beta - b0) %*% B0 %*% (beta - b0)
+
+ return(sum(Xb.reform - denom) + log.prior)
+}
+
+## MNL log-posterior function (used to get starting values)
+## matrix Y with NAs
+"mnl.logpost.NA" <- function(beta, Y, X, b0, B0){
+ k <- ncol(X)
+ numer <- exp(X %*% beta)
+ numer[Y== -999] <- NA
+ numer.mat <- matrix(numer, nrow(Y), ncol(Y), byrow=TRUE)
+ denom <- apply(numer.mat, 1, sum, na.rm=TRUE)
+ choice.probs <- numer.mat / denom
+ Yna <- Y
+ Yna[Y == -999] <- NA
+ log.like.mat <- log(choice.probs) * Yna
+ log.like <- sum(apply(log.like.mat, 1, sum, na.rm=TRUE))
+
+ log.prior <- 0.5 * t(beta - b0) %*% B0 %*% (beta - b0)
+
+ return(log.like + log.prior)
}
@@ -189,12 +211,16 @@
"MCMCmnl" <-
function(formula, baseline=NULL, data = parent.frame(),
burnin = 1000, mcmc = 10000, thin=1,
- mcmc.method = "MH", tune = 1.1, verbose = 0, seed = NA,
+ mcmc.method = c("IndMH", "RWM", "slice"),
+ tune = 1.0, tdf=6, verbose = 0, seed = NA,
beta.start = NA, b0 = 0, B0 = 0, ...) {
## checks
check.offset(list(...))
check.mcmc.parameters(burnin, mcmc, thin)
+ if (tdf <= 0){
+ stop("degrees of freedom for multivariate-t proposal must be positive.\n Respecify tdf and try again.\n")
+ }
## seeds
seeds <- form.seeds(seed)
@@ -225,15 +251,32 @@
## form the tuning parameter
tune <- vector.tune(tune, K)
+
+ ## priors and starting values
+ mvn.prior <- form.mvn.prior(b0, B0, K)
+ b0 <- mvn.prior[[1]]
+ B0 <- mvn.prior[[2]]
+
beta.init <- rep(0, K)
cat("Calculating MLEs and large sample var-cov matrix.\n")
cat("This may take a moment...\n")
- optim.out <- optim(beta.init, mnl.loglike, method="BFGS",
- control=list(fnscale=-1),
- hessian=TRUE, Y=Y, X=X)
- V <- solve(-1*optim.out$hessian)
-
- ## starting values and priors
+ if (max(is.na(Y))){
+ optim.out <- optim(beta.init, mnl.logpost.NA, method="BFGS",
+ control=list(fnscale=-1),
+ hessian=TRUE, Y=Y, X=X, b0=b0, B0=B0)
+ }
+ else{
+ new.Y <- apply(Y==1, 1, which)
+ optim.out <- optim(beta.init, mnl.logpost.noNA, method="BFGS",
+ control=list(fnscale=-1),
+ hessian=TRUE, new.Y=new.Y, X=X, b0=b0, B0=B0)
+ }
+ cat("Inverting Hessian to get large sample var-cov matrix.\n")
+ ##V <- solve(-1*optim.out$hessian)
+ V <- chol2inv(chol(-1*optim.out$hessian))
+ beta.mode <- matrix(optim.out$par, K, 1)
+
+
if (is.na(beta.start) || is.null(beta.start)){
beta.start <- matrix(optim.out$par, K, 1)
}
@@ -243,14 +286,12 @@
else if (length(beta.start != K)){
stop("beta.start not of appropriate dimension\n")
}
- mvn.prior <- form.mvn.prior(b0, B0, K)
- b0 <- mvn.prior[[1]]
- B0 <- mvn.prior[[2]]
## define holder for posterior sample
sample <- matrix(data=0, mcmc/thin, dim(X)[2] )
+ posterior <- NULL
- if (mcmc.method=="MH"){
+ if (mcmc.method=="RWM"){
## call C++ code to draw sample
auto.Scythe.call(output.object="posterior", cc.fun.name="MCMCmnlMH",
sample.nonconst=sample, Y=Y, X=X,
@@ -260,14 +301,32 @@
seedarray=as.integer(seed.array),
lecuyerstream=as.integer(lecuyer.stream),
verbose=as.integer(verbose),
- betastart=beta.start,
+ betastart=beta.start, betamode=beta.mode,
b0=b0, B0=B0,
- V=V)
+ V=V, RW=as.integer(1), tdf=as.double(tdf))
## put together matrix and build MCMC object to return
output <- form.mcmc.object(posterior, names=xnames,
title="MCMCmnl Posterior Sample")
}
+ else if (mcmc.method=="IndMH"){
+ auto.Scythe.call(output.object="posterior", cc.fun.name="MCMCmnlMH",
+ sample.nonconst=sample, Y=Y, X=X,
+ burnin=as.integer(burnin),
+ mcmc=as.integer(mcmc), thin=as.integer(thin),
+ tune=tune, lecuyer=as.integer(lecuyer),
+ seedarray=as.integer(seed.array),
+ lecuyerstream=as.integer(lecuyer.stream),
+ verbose=as.integer(verbose),
+ betastart=beta.start, betamode=beta.mode,
+ b0=b0, B0=B0,
+ V=V, RW=as.integer(0), tdf=as.double(tdf))
+
+ ## put together matrix and build MCMC object to return
+ output <- form.mcmc.object(posterior, names=xnames,
+ title="MCMCmnl Posterior Sample")
+
+ }
else if (mcmc.method=="slice"){
## call C++ code to draw sample
auto.Scythe.call(output.object="posterior", cc.fun.name="MCMCmnlslice",
diff --git a/R/MCMCoprobit.R b/R/MCMCoprobit.R
index 0e96596..33eb794 100644
--- a/R/MCMCoprobit.R
+++ b/R/MCMCoprobit.R
@@ -8,9 +8,8 @@
"MCMCoprobit" <-
function(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
- thin = 1, tune = NA, verbose = 0, seed = NA, beta.start = NA,
- b0 = 0, B0 = 0, ...) {
-
+ thin = 1, tune = NA, tdf = 1, verbose = 0, seed = NA, beta.start = NA,
+ b0 = 0, B0 = 0, a0 = 0, A0 = 0, mcmc.method = c("Cowles", "AC"), ...) {
## checks
check.offset(list(...))
@@ -21,16 +20,15 @@
lecuyer <- seeds[[1]]
seed.array <- seeds[[2]]
lecuyer.stream <- seeds[[3]]
-
- ## extract X, Y, and variable names from the model formula and frame
- call <- match.call()
- mt <- terms(formula, data=data)
- if(missing(data)) data <- sys.frame(sys.parent())
+ ## extract X, Y, and variable names from the model formula and frame
+ call <- match.call()
+ mt <- terms(formula, data=data)
+ if(missing(data)) data <- sys.frame(sys.parent())
mf <- match.call(expand.dots = FALSE)
mf$burnin <- mf$mcmc <- mf$b0 <- mf$B0 <- NULL
- mf$thin <- mf$... <- mf$tune <- mf$verbose <- mf$seed <- NULL
- mf$beta.start <- NULL
+ mf$thin <- mf$... <- mf$tune <- mf$tdf <- mf$verbose <- mf$seed <- NULL
+ mf$beta.start <- mf$mcmc.method <- NULL
mf$drop.unused.levels <- TRUE
mf[[1]] <- as.name("model.frame")
mf <- eval(mf, sys.frame(sys.parent()))
@@ -43,8 +41,8 @@
Y <- factor(Y, ordered=TRUE)
ncat <- nlevels(Y) # number of categories of y
cat <- levels(Y) # values of categories of y
- N <- nrow(X) # number of observations
- K <- ncol(X) # number of covariates
+ N <- nrow(X) # number of observations
+ K <- ncol(X) # number of covariates
if (length(Y) != N){
cat("X and Y do not have same number of rows.\n")
stop("Please respecify and call MCMCoprobit() again.\n")
@@ -102,32 +100,65 @@
cat("N(b0,B0) prior B0 not conformable.\n")
stop("Please respecify and call MCMCoprobit() again.\n")
}
-
+
+ ## prior for alpha error checking
+ if(is.null(dim(a0))) {
+ a0 <- a0 * matrix(1, ncat-1, 1)
+ }
+ if((dim(a0)[1] != ncat-1) || (dim(a0)[2] != 1)) {
+ cat("N(a0,A0) prior a0 not conformable.\n")
+ stop("Please respecify and call MCMCoprobit() again.\n")
+ }
+ if(is.null(dim(A0))) {
+ A0 <- A0 + diag(ncat - 1)
+ }
+ if((dim(A0)[1] != ncat - 1) || (dim(A0)[2] != ncat - 1)) {
+ cat("N(a0, A0) prior A0 not conformable.\n")
+ stop("Please respecify and call MCMCoprobit() again.\n")
+ }
+
## form gamma starting values (note: not changeable)
gamma <- matrix(NA,ncat+1,1)
gamma[1] <- -300
gamma[2] <- 0
gamma[3:ncat] <- (polr.out$zeta[2:(ncat-1)] - polr.out$zeta[1])*.588
gamma[ncat+1] <- 300
-
+
## posterior sample
sample <- matrix(data=0, mcmc/thin, K + ncat + 1)
## call C++ code to draw sample
+ nY <- as.matrix(as.numeric(Y))
+
+ ## mcmc.method
+ cowles <- as.integer(1)
+ if(mcmc.method[1]=="AC") {cowles <- as.integer(0)}
+
+ ## form the tuning parameter
+ tune <- vector.tune(tune, ncat-1)
+ posterior <- NULL
+
auto.Scythe.call(output.object="posterior", cc.fun.name="MCMCoprobit",
- sample.nonconst=sample, Y=as.integer(Y), X=X,
+ sample.nonconst=sample, Y=as.integer(Y), nY=nY, X=X,
burnin=as.integer(burnin),
mcmc=as.integer(mcmc), thin=as.integer(thin),
- tune=as.double(tune), lecuyer=as.integer(lecuyer),
+ tune=tune, tdf=as.double(tdf),
+ lecuyer=as.integer(lecuyer),
seedarray=as.integer(seed.array),
lecuyerstream=as.integer(lecuyer.stream),
verbose=as.integer(verbose), beta=beta.start,
- gamma=gamma, b0=b0, B0=B0)
+ gamma=gamma, b0=b0, B0=B0, a0=a0, A0=A0,
+ cowles=as.integer(cowles))
## put together matrix and build MCMC object to return
sample <- matrix(posterior$sampledata, posterior$samplerow,
posterior$samplecol, byrow=FALSE)
- sample <- sample[,c(1:K, (K+3):(K+ncat))]
+ if(mcmc.method[1]=="AC"){
+ sample[ , 1] <- sample[, 1] - sample[, K+2] ## post-MCMC normalization
+ sample[ , (K+2):(K+ncat)] <- sample[ , (K+2):(K+ncat)] - sample[, K+2] ## post-MCMC normalization
+ }
+ sample <- sample[, c(1:K, (K+3):(K+ncat))]
+
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
diff --git a/R/MCMCpoisson.R b/R/MCMCpoisson.R
index e86b9fa..bef353e 100644
--- a/R/MCMCpoisson.R
+++ b/R/MCMCpoisson.R
@@ -87,7 +87,7 @@
logpost.poisson(theta.tilde, Y, X, b0, B0)
}
-
+ posterior <- NULL
## call C++ code to draw sample
auto.Scythe.call(output.object="posterior", cc.fun.name="MCMCpoisson",
diff --git a/R/MCMCpoissonChangepoint.R b/R/MCMCpoissonChangepoint.R
old mode 100644
new mode 100755
index 4669cb7..348ee34
--- a/R/MCMCpoissonChangepoint.R
+++ b/R/MCMCpoissonChangepoint.R
@@ -1,148 +1,102 @@
-"MCMCpoissonChangepoint" <-
- function (data, m = 1, burnin = 1000, mcmc = 1000, thin = 1, verbose = 0,
- seed = NA, c0, d0, a = NULL, b = NULL, marginal.likelihood = c("none", "Chib95"), ...)
-{
+## sample from the posterior distribution
+## of a Poisson model with multiple changepoints
+## using linked C++ code in Scythe 1.0
+##
+## JHP 07/01/2007
+##
+## Revised on 09/12/2007 JHP
+
+"MCMCpoissonChangepoint"<-
+ function(data, m = 1, c0 = NA, d0 = NA, a = NULL, b = NULL,
+ burnin = 10000, mcmc = 10000, thin = 1, verbose = 0,
+ seed = NA, lambda.start = NA, P.start = NA,
+ marginal.likelihood = c("none", "Chib95"), ...) {
+
+ ## check iteration parameters
check.mcmc.parameters(burnin, mcmc, thin)
totiter <- mcmc + burnin
cl <- match.call()
- if (!is.na(seed))
- set.seed(seed)
- y <- data
- n <- length(y)
- A0 <- trans.mat.prior(m = m, n = n, a = a, b = b)
- marginal.likelihood <- match.arg(marginal.likelihood)
- lambda.store <- matrix(NA, mcmc/thin, m + 1)
- P.store <- matrix(NA, mcmc/thin, (m + 1)^2)
- ps.store <- matrix(0, n, m + 1)
- s1.store <- matrix(NA, mcmc/thin, n)
- py <- rep(0, m + 1)
- pdf.P.store <- matrix(NA, mcmc/thin, m + 1)
- lambda1 <- rep(NA, m + 1)
- P1 <- matrix(NA, m + 1, m + 1)
- lambda0 <- runif(m + 1)
- P0 <- trans.mat.prior(m = m, n = n, a = 0.9, b = 0.1)
- for (iter in 1:totiter) {
- state.out <- Poisson.state.sampler(m = m, y = y, lambda = lambda0,
- P = P0)
- s1 <- state.out$s1
- ps1 <- state.out$ps1
- for (j in 1:(m + 1)) {
- ej <- as.numeric(s1 == j)
- yj <- y[ej == 1]
- nj <- length(yj)
- c1 <- sum(yj) + c0
- d1 <- nj + d0
- lambda1[j] <- rgamma(1, c1, d1)
- }
- switch <- switchg(s1)
- for (j in 1:(m + 1)) {
- switch1 <- A0[j, ] + switch[j, ]
- pj <- rdirichlet.cp(1, switch1)
- P1[j, ] <- pj
- }
- lambda0 <- lambda1
- P0 <- P1
- if (iter > burnin && (iter%%thin == 0)) {
- lambda.store[(iter - burnin)/thin, ] <- lambda1
- P.store[(iter - burnin)/thin, ] <- as.vector(t(P1))
- s1.store[(iter - burnin)/thin, ] <- s1
- ps.store <- ps.store + ps1
- }
- if (verbose > 0 && iter%%verbose == 0) {
- cat("----------------------------------------------",
- "\n")
- cat("iteration = ", iter, "\n")
- cat("lambda = ", lambda1, "\n")
- cat("Transition Matrix", "\n")
- for (i in 1:(m + 1)) cat(paste("", P1[i, ]), fill = TRUE,
- labels = paste("{", i, "}:", sep = ""), sep = ",")
- }
- }
- if (marginal.likelihood == "Chib95") {
- lambda.st <- apply(lambda.store, 2, mean)
- P.vec.st <- apply(P.store, 2, mean)
- P.st <- t(matrix(P.vec.st, m + 1, m + 1))
- density.lambda <- matrix(NA, (mcmc/thin), m + 1)
- for (i in 1:(mcmc/thin)) {
- for (j in 1:(m + 1)) {
- ej <- as.numeric(s1.store[i, ] == j)
- yj <- y[ej == 1]
- nj <- length(yj)
- c1 <- sum(yj) + c0
- d1 <- nj + d0
- density.lambda[i, j] <- dgamma(lambda.st[j],
- c1, d1)
- }
- }
- pdf.lambda <- log(prod(apply(density.lambda, 2, mean)))
- for (g in 1:(mcmc/thin)) {
- state.out <- Poisson.state.sampler(m = m, y = y,
- lambda = lambda.st, P = P0)
- s1 <- state.out$s1
- ps1 <- state.out$ps1
- switch <- switchg(s1)
- for (j in 1:(m + 1)) {
- switch1 <- A0[j, ] + switch[j, ]
- pj <- rdirichlet.cp(1, switch1)
- P1[j, ] <- pj
- pdf.P.store[g, j] <- ddirichlet.cp(P.st[j, ],
- switch1)
- }
- P0 <- P1
- }
- pdf.P <- log(prod(apply(pdf.P.store, 2, mean)))
- F <- matrix(NA, n, m + 1)
- like <- rep(NA, n)
- pr1 <- c(1, rep(0, m))
- for (t in 1:n) {
- py <- sapply(c(1:(m + 1)), function(i) {
- poisson.pdf(y[t], lambda.st[i])
- })
- if (t == 1) {
- pstyt1 = pr1
- }
- else {
- pstyt1 <- F[t - 1, ] %*% P.st
- }
- unnorm.pstyt <- pstyt1 * py
- pstyt <- unnorm.pstyt/sum(unnorm.pstyt)
- F[t, ] <- pstyt
- like[t] <- sum(unnorm.pstyt)
- }
- loglik <- sum(log(like))
- nprior.lambda <- nprior.P <- rep(NA, m + 1)
- nprior.lambda <- sapply(c(1:(m + 1)), function(i) {
- dgamma(lambda.st[i], c0, d0, log = TRUE)
- })
- nprior.P <- sapply(c(1:(m + 1)), function(i) {
- log(ddirichlet.cp(P.st[i, ], A0[i, ]))
- })
- prior.lambda <- sum(nprior.lambda)
- prior.P <- sum(nprior.P)
- numerator <- loglik + prior.lambda + prior.P
- denominator <- pdf.lambda + pdf.P
- logmarglike <- numerator - denominator
- if (verbose > 0) {
- cat("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n")
- cat("Log Marginal Likelihood\n")
- cat("-------------------------------------------------",
- "\n")
- cat("log(marglike)= ", logmarglike, "\n")
- cat("log(likelihood)= ", loglik, "\n")
- cat("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n")
- }
+
+ ## ## seeds
+ seeds <- form.seeds(seed)
+ lecuyer <- seeds[[1]]
+ seed.array <- seeds[[2]]
+ lecuyer.stream <- seeds[[3]]
+
+ ## sample size
+ y <- as.matrix(data)
+ n <- nrow(y)
+ ns <- m+1
+
+ ## prior
+ A0 <- trans.mat.prior(m=m, n=n, a=a, b=b)
+ if (is.na(c0)||is.na(d0))
+ stop("Please specify prior for lambda (c0 and d0) and call MCMCpoissonChangepoint again.\n")
+
+ ## get marginal likelihood argument
+ marginal.likelihood <- match.arg(marginal.likelihood)
+
+ ## following MCMCregress, set chib as binary
+ logmarglike <- NULL
+ chib <- 0
+ if (marginal.likelihood == "Chib95"){
+ chib <- 1
}
- else {
- logmarglike <- marginal <- loglik <- NULL
+
+ Pstart <- check.P(P.start, m=m, n=n, a=a, b=b)
+ lambdastart <- check.theta(lambda.start, ns, y, min=range(y)[1], max=range(y)[2])
+
+ nstore <- mcmc/thin
+
+ ## call C++ code to draw sample
+ posterior <- .C("MCMCpoissonChangepoint",
+ lambdaout = as.double(rep(0.0, nstore*ns)),
+ Pout = as.double(rep(0.0, nstore*ns*ns)),
+ psout = as.double(rep(0.0, n*ns)),
+ sout = as.double(rep(0.0, nstore*n)),
+ Ydata = as.double(y),
+ Yrow = as.integer(nrow(y)),
+ Ycol = as.integer(ncol(y)),
+ m = as.integer(m),
+ 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),
+ lambdastart = as.double(lambdastart),
+ Pstart = as.double(Pstart),
+ a = as.double(a),
+ b = as.double(b),
+ c0 = as.double(c0),
+ d0 = as.double(d0),
+ A0data = as.double(A0),
+ logmarglikeholder = as.double(0.0),
+ chib = as.integer(chib))
+
+ ## get marginal likelihood if Chib95
+ if (marginal.likelihood == "Chib95"){
+ logmarglike <- posterior$logmarglikeholder
+ ##loglike <- posterior$loglikeholder
}
- output <- as.mcmc(lambda.store)
- varnames(output) <- paste("lambda.", 1:(m + 1), sep = "")
- attr(output, "title") <- "MCMCpoissonChangepoint Posterior Sample"
- attr(output, "y") <- data
+
+ ## pull together matrix and build MCMC object to return
+ lambda.holder <- matrix(posterior$lambdaout, mcmc/thin)
+ P.holder <- matrix(posterior$Pout, mcmc/thin)
+ s.holder <- matrix(posterior$sout, mcmc/thin)
+ ps.holder <- matrix(posterior$psout, n)
+
+ output <- mcmc(data=lambda.holder, start=burnin+1, end=burnin + mcmc, thin=thin)
+ varnames(output) <- paste("lambda.", 1:ns, sep = "")
+ attr(output,"title") <- "MCMCpoissonChangepoint Posterior Sample"
+ attr(output, "y") <- y
+ attr(output, "m") <- m
attr(output, "call") <- cl
attr(output, "logmarglike") <- logmarglike
- attr(output, "loglik") <- loglik
- attr(output, "prob.state") <- ps.store/(mcmc/thin)
+ attr(output, "prob.state") <- ps.holder/(mcmc/thin)
+ attr(output, "s.store") <- s.holder
return(output)
-}
+
+ }## end of MCMC function
diff --git a/R/MCMCprobit.R b/R/MCMCprobit.R
index 395b47c..a2edfd8 100644
--- a/R/MCMCprobit.R
+++ b/R/MCMCprobit.R
@@ -92,6 +92,7 @@
logpost.probit(theta.tilde, Y, X, b0, B0)
}
+ posterior <- NULL
if (is.null(resvec)){
## define holder for posterior density sample
diff --git a/R/MCMCregress.R b/R/MCMCregress.R
index fd92ef5..68bf281 100644
--- a/R/MCMCregress.R
+++ b/R/MCMCregress.R
@@ -64,7 +64,8 @@
## define holder for posterior sample
sample <- matrix(data=0, mcmc/thin, K+1)
-
+ posterior <- NULL
+
## call C++ code to draw sample
auto.Scythe.call(output.object="posterior", cc.fun.name="MCMCregress",
sample.nonconst=sample, Y=Y, X=X,
diff --git a/R/MCMCtobit.R b/R/MCMCtobit.R
index ba35697..e5735fd 100644
--- a/R/MCMCtobit.R
+++ b/R/MCMCtobit.R
@@ -1,65 +1,66 @@
-"MCMCtobit" <-
-function(formula, data=parent.frame(), 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, ...) {
-
- # checks
- check.offset(list(...))
- check.mcmc.parameters(burnin, mcmc, thin)
- if (!is.numeric(below) | !is.numeric(above)) {
- cat("Error: Censoring points must be numeric, which includes +-Inf.\n")
- stop("Please respecify and call ", calling.function(), " again.",
- call.=FALSE)
- }
- if (below >= above) {
- cat("Error: Truncation points are logically inconsistent.\n")
- stop("Please respecify and call ", calling.function(), " again.",
- call.=FALSE)
- }
-
- # convert infinite values to finite approximations
- if(is.infinite(below)) below <- .Machine$double.xmax*-1
- if(is.infinite(above)) above <- .Machine$double.xmax
-
- # seeds
- seeds <- form.seeds(seed)
- lecuyer <- seeds[[1]]
- seed.array <- seeds[[2]]
- lecuyer.stream <- seeds[[3]]
-
- # form response and model matrices
- holder <- parse.formula(formula, data)
- Y <- holder[[1]]
- X <- holder[[2]]
- xnames <- holder[[3]]
- K <- ncol(X) # number of covariates
-
- # starting values and priors
- beta.start <- coef.start(beta.start, K, formula, family=gaussian, data)
- mvn.prior <- form.mvn.prior(b0, B0, K)
- b0 <- mvn.prior[[1]]
- B0 <- mvn.prior[[2]]
- check.ig.prior(c0, d0)
-
- # define holder for posterior sample
- sample <- matrix(data=0, mcmc/thin, K+1)
-
- # call C++ code to draw sample
- auto.Scythe.call(output.object="posterior", cc.fun.name="MCMCtobit",
- sample.nonconst=sample, Y=Y, X=X, below=as.double(below),
- above=as.double(above),
- burnin=as.integer(burnin), mcmc=as.integer(mcmc),
- thin=as.integer(thin),
- lecuyer=as.integer(lecuyer),
- seedarray=as.integer(seed.array),
- lecuyerstream=as.integer(lecuyer.stream),
- verbose=as.integer(verbose), betastart=beta.start,
- b0=b0, B0=B0, c0=as.double(c0), d0=as.double(d0))
-
- # pull together matrix and build MCMC object to return
- output <- form.mcmc.object(posterior,
- names=c(xnames, "sigma2"),
- title="MCMCtobit Posterior Sample")
- return(output)
- }
+"MCMCtobit" <-
+function(formula, data=parent.frame(), 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, ...) {
+
+ # checks
+ check.offset(list(...))
+ check.mcmc.parameters(burnin, mcmc, thin)
+ if (!is.numeric(below) | !is.numeric(above)) {
+ cat("Error: Censoring points must be numeric, which includes +-Inf.\n")
+ stop("Please respecify and call ", calling.function(), " again.",
+ call.=FALSE)
+ }
+ if (below >= above) {
+ cat("Error: Truncation points are logically inconsistent.\n")
+ stop("Please respecify and call ", calling.function(), " again.",
+ call.=FALSE)
+ }
+
+ # convert infinite values to finite approximations
+ if(is.infinite(below)) below <- .Machine$double.xmax*-1
+ if(is.infinite(above)) above <- .Machine$double.xmax
+
+ # seeds
+ seeds <- form.seeds(seed)
+ lecuyer <- seeds[[1]]
+ seed.array <- seeds[[2]]
+ lecuyer.stream <- seeds[[3]]
+
+ # form response and model matrices
+ holder <- parse.formula(formula, data)
+ Y <- holder[[1]]
+ X <- holder[[2]]
+ xnames <- holder[[3]]
+ K <- ncol(X) # number of covariates
+
+ # starting values and priors
+ beta.start <- coef.start(beta.start, K, formula, family=gaussian, data)
+ mvn.prior <- form.mvn.prior(b0, B0, K)
+ b0 <- mvn.prior[[1]]
+ B0 <- mvn.prior[[2]]
+ check.ig.prior(c0, d0)
+
+ # define holder for posterior sample
+ sample <- matrix(data=0, mcmc/thin, K+1)
+ posterior <- NULL
+
+ # call C++ code to draw sample
+ auto.Scythe.call(output.object="posterior", cc.fun.name="MCMCtobit",
+ sample.nonconst=sample, Y=Y, X=X, below=as.double(below),
+ above=as.double(above),
+ burnin=as.integer(burnin), mcmc=as.integer(mcmc),
+ thin=as.integer(thin),
+ lecuyer=as.integer(lecuyer),
+ seedarray=as.integer(seed.array),
+ lecuyerstream=as.integer(lecuyer.stream),
+ verbose=as.integer(verbose), betastart=beta.start,
+ b0=b0, B0=B0, c0=as.double(c0), d0=as.double(d0))
+
+ # pull together matrix and build MCMC object to return
+ output <- form.mcmc.object(posterior,
+ names=c(xnames, "sigma2"),
+ title="MCMCtobit Posterior Sample")
+ return(output)
+ }
diff --git a/R/btsutil.R b/R/btsutil.R
index 6babf17..c49eecd 100644
--- a/R/btsutil.R
+++ b/R/btsutil.R
@@ -6,6 +6,8 @@
## 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.
@@ -13,19 +15,8 @@
## Helper functions for MCMCPoissonChangepoint()
##############################################################
- "switchg" <- function(s1){
-
- ## "switchg" computes the frequency of elements in state vector
- ## for one step ahead Markov process only
- ## example: s1 <- c(1,1,1,2,2,2,3,4,4,4,4,5,5)
- ## switchg(s1)
- ## [,1] [,2] [,3] [,4] [,5]
- ## [1,] 2 1 0 0 0
- ## [2,] 0 2 1 0 0
- ## [3,] 0 0 0 1 0
- ## [4,] 0 0 0 3 1
- ## [5,] 0 0 0 0 1
-
+ ## switch a state vector into a matrix containing the number of states
+ "switchg" <- function(s1){
s <- max(s1)
out <- matrix(0,s,s)
@@ -39,16 +30,6 @@
}
## "trans.mat.prior" makes a transition matrix
-
- ## a and b are beta prior for transition matrix
- ## By default, the expected duration is computed and
- ## corresponding a and b values are assigned. The expected
- ## duration is the sample period divided by the number of states.
- ## ex) 300 years and 2 changepoints (3 states)
- ## expected.duration <- 300/(2+1)
- ## b <- 0.1
- ## a <- b*expected.duration
-
"trans.mat.prior" <- function(m, n, a=NULL, b=NULL){
if (!is.null(a)|!is.null(b)){
a <- a
@@ -66,197 +47,101 @@
return(trans)
}
- "rdirichlet.cp" <- function(n, alpha){
- ## "rdirichlet.cp" picks n random deviates from the Dirichlet function
- ## with shape parameters alpha
- ## Note that alpha can contain zero to deal with transition matrix rowwise.
- ## It returns a matrix.
-
- col <- length(alpha)
- out <- matrix(NA, n, col)
- out[,which(alpha==0)]<- 0
- a <- alpha[alpha!=0]
- l <- length(a);
- x <- matrix(rgamma(l*n,a), ncol=l, byrow=TRUE);
- sm <- x%*%rep(1,l);
- dir.out <- x/as.vector(sm);
- ## combine zero and nonzero parts prob
- out[,which(alpha!=0)] <- dir.out
- return(out)
- }
-
- "ddirichlet.cp" <- function (x, alpha){
- ## "ddirichlet.cp" is a density function for the Dirichlet distribution
- ## alpha and x can contain zero
- ## but zero should be in the same place
- ## x=c(0.9, 0.1, 0) and alpha=c(6,1,0)
- ## It returns a matrix
- col <- length(alpha)
- out <- rep(NA, col)
- out[which(alpha==0)] <- 0
- a <- alpha[alpha!=0]
- x1 <- x[x!=0]
- dirichlet1 <- function(x1, a) {
- logD <- sum(lgamma(a)) - lgamma(sum(a))
- s <- sum((a - 1) * log(x1))
- exp(sum(s) - logD)
- }
- if (!is.matrix(x1))
- if (is.data.frame(x1))
- x1 <- as.matrix(x1)
- else x1 <- t(x1)
- if (!is.matrix(a))
- a <- matrix(a, ncol = length(a), nrow = nrow(x1),
- byrow = TRUE)
- if (any(dim(x1) != dim(a)))
- stop("Mismatch between dimensions of x and alpha in
- ddirichlet().\n")
- pd <- vector(length = nrow(x1))
- for (i in 1:nrow(x1)) pd[i] <- dirichlet1(x1[i, ], a[i,])
- pd[apply(x1, 1, function(z) any(z < 0 | z > 1))] <- 0
- pd[apply(x1, 1, function(z) all.equal(sum(z), 1) != TRUE)] <- 0
- return(pd)
- }
-
- "poisson.pdf" <- function(y, lambda){
- ## What this function does is same with "prod(dpois(y, lambda))"
- n <- length(y)
- log.like <- sum(sapply(c(1:n), function(i){-lambda + y[i]*log(lambda)- log(factorial(y[i]))}))
- return(exp(log.like))
- }
-
- "Poisson.state.sampler" <- function(m, y, lambda, P){
- ## "Poisson.state.sampler" samples state vector (s1) by forward and backward sampling
- ## P is a (m+1 by m+1) transition matrix
- ## F matrix that contains all the information of Pr(st|Yt)
-
- ## Forward sampling: update F matrix
- n <- length(y)
- F <- matrix(NA, n, m+1) # storage for the Filtered probabilities
- pr1 <- c(1,rep(0, m)) # initial probability Pr(s1=k|Y0, lambda)
- for (t in 1:n){
- py <- sapply(c(1:(m+1)), function(i){poisson.pdf(y[t], lambda[i])})
- if(t==1) pstyt1 = pr1
- else pstyt1 <- F[t-1,]%*%P
- unnorm.pstyt <- pstyt1*py
- pstyt <- unnorm.pstyt/sum(unnorm.pstyt) # Pr(st|Yt)
- F[t,] <- pstyt
- }
-
- ## Backward recursions: sample s1 using F matrix and a transition matrix (P)
- s1 <- matrix(NA, n, 1) ## holder for state variables
- ps1 <- matrix(NA, n, m+1) ## holder for state probabilities
- ps1[n,] <- F[n,] ## we know last elements of ps1 and s1
- s1[n,1] <- m+1
- t <- n-1
-
- while (t>=1){
- st1 <- s1[t+1]
- unnorm.pstyn <- F[t,]*P[,st1]
- ## normalize into a prob. density
- pstyn <- unnorm.pstyn/sum(unnorm.pstyn)
- if (st1==1) {s1[t]<-1}
- ## If this is not the first period,
- ## draw a state variable from a discrete prob. distribution("pstyn")
- ## using the inverse CDF method.
-
- ## What inverse CDF does is as follows
- ## In backward recursion, the current state (t) is always at st1, which is one of (1...m+1)
- ## The probability that transition happens at t-1, that is Pr(t-1=st1-1|t=st1) can be found
- ## by using (dicrete case) inverse CDF method for the probability vector we have (pstyn).
- ## draw u ~ U[0,1] then, if u<pstyn[st1-1], t-1=st1-1, otherwise t-1=st1
- ## This is because Pr(x<t)=Pr(F^-1(u)<t)=Pr(u<F(t))=F(t)
- ## Thus, when u ~ U[0,1], x=F^-1(u) ~ F
-
- else {
- pone <- pstyn[st1-1]
- s1[t] <- ifelse (runif(1) < pone, st1-1, st1)}
- ps1[t,] <- pstyn # probabilities pertaining to a certain state
- t <- t-1
- }
- ## name and report outputs
- out <- list(s1, ps1)
- names(out)<-c("s1","ps1")
- return(out)
- }
-
-
-##############################################################
-## Plot functions for MCMCPoissonChangepoint()
-##############################################################
-
- ## "plotPostState" draws a plot of posterior distribution of states
-
- "plotPostState" <-
- function(mcmcout, y, m, legend.control=NULL, cex=0.8, lwd=1.2){
-
- ## To use legend.control, a user has to specify xy coordinate for legend.control.
- ## ex) legend.control = c(x, y)
-
- ## extract posterior state probability vector
+ ## "plotState" draws a plot of posterior distribution of states
+ "plotState" <-
+ function (mcmcout, legend.control = NULL, cex = 0.8, lwd = 1.2)
+ {
out <- attr(mcmcout, "prob.state")
-
- ## y must be a ts object. Otherwise, it is coerced as a ts object.
- if (!is.ts(y)) y <- ts(y)
- time.frame <- as.vector(time(y))
-
- ## draw a plot for each state variable
- plot(1, 0, xlim=range(time.frame), ylim=c(0,1), type="n", main="", xlab="Time", cex=cex, lwd=lwd,
- ylab=expression(paste("Pr(", S[t], "= k |", Y[t], ")")))
- for (i in 1:(m+1))
- points(time.frame, out[,i], type="o", lty=i, lwd=lwd, col=i, cex=cex)
-
- ## put a legend to indicate each line for the state
- if (!is.null(legend.control)){
- if (length(legend.control)!=2)
+ y <- attr(mcmcout, "y")
+ m <- attr(mcmcout, "m")
+ if (!is.ts(y))
+ y <- ts(y)
+ time.frame <- as.vector(time(y))
+ plot(1, 0, xlim = range(time.frame), ylim = c(0, 1), type = "n",
+ main = "", xlab = "Time", cex = cex, lwd = lwd, ylab = expression(paste("Pr(",
+ S[t], "= k |", Y[t], ")")))
+ for (i in 1:(m + 1)) points(time.frame, out[, i], type = "o",
+ lty = i, lwd = lwd, col = i, cex = cex)
+ if (!is.null(legend.control)) {
+ if (length(legend.control) != 2)
stop("You should specify x and y coordinate for a legend.")
- else
- legend(legend.control[1], legend.control[2], legend = paste("State", 1:(m+1), sep=""),
- col=1:(m+1), lty=1:(m+1), lwd = rep(lwd, m+1), pch=rep(1, m+1), bty="n")
- }
+ else legend(legend.control[1], legend.control[2], legend = paste("State",
+ 1:(m + 1), sep = ""), col = 1:(m + 1), lty = 1:(m +
+ 1), lwd = rep(lwd, m + 1), pch = rep(1, m + 1), bty = "n")
}
-
- ## "plotPostChangepoint" draws a plot of posterior changepoint probability
- "plotPostChangepoint" <-
- function(mcmcout, y, m, xlab="Time", ylab=""){
+}
- ## extract posterior state probability vector
+ ## "plotChangepoint" draws a plot of posterior changepoint probability
+ "plotChangepoint" <-
+ function (mcmcout, xlab = "Time", ylab = "")
+ {
out <- attr(mcmcout, "prob.state")
-
- ## y must be a ts object. Otherwise, it is coerced as a ts object.
- if (!is.ts(y)) y <- ts(y)
- time.frame <- as.vector(time(y))
-
- ## divide space
- if (m == 1){
- pr.st <- c(0, diff(out[,(m+1)]))
- plot(time.frame, pr.st, type="h", main="", xlab=xlab, ylab=ylab)
- cp <- which(cumsum(pr.st) > 0.5)[1]-1
- abline(v=cp+ time.frame[1], lty=3, col="red")
- }
- else {
- par(mfrow=c(m,1));par(mar=c(2,4,1,1))
- cp <- rep(NA, m) # holder for expected changepoints
- for (i in 2:m){
- pr.st <- c(0, diff(out[,i]))
- pr.st <- ifelse(pr.st <0, 0, pr.st)
- plot(time.frame, pr.st, type="h", main="", xlab=xlab, ylab=ylab)
- cp[i-1] <- which(cumsum(pr.st)>0.5)[1]-1
- abline(v=cp[i-1]+ time.frame[1], lty=3, col="red")
-
- }
- pr.st <- c(0, diff(out[,(m+1)]))
- plot(time.frame, pr.st, type="h", main="", xlab=xlab, ylab=ylab)
- cp[m] <- which(cumsum(pr.st)>0.5)[1]-1
- abline(v=cp[m]+ time.frame[1], lty=3, col="red")
+ y <- attr(mcmcout, "y")
+ m <- attr(mcmcout, "m")
+ if (!is.ts(y))
+ y <- ts(y)
+ time.frame <- as.vector(time(y))
+ if (m == 1) {
+ pr.st <- c(0, diff(out[, (m + 1)]))
+ plot(time.frame, pr.st, type = "h", main = "", xlab = xlab,
+ ylab = ylab)
+ cp <- which(cumsum(pr.st) > 0.5)[1] - 1
+ abline(v = cp + time.frame[1], lty = 3, col = "red")
+ }
+ else {
+ par(mfrow = c(m, 1))
+ par(mar = c(2, 4, 1, 1))
+ cp <- rep(NA, m)
+ for (i in 2:m) {
+ pr.st <- c(0, diff(out[, i]))
+ pr.st <- ifelse(pr.st < 0, 0, pr.st)
+ plot(time.frame, pr.st, type = "h", main = "", xlab = xlab,
+ ylab = ylab)
+ cp[i - 1] <- which(cumsum(pr.st) > 0.5)[1] - 1
+ abline(v = cp[i - 1] + time.frame[1], lty = 3, col = "red")
}
-
- cp.means <- rep(NA, m+1)
+ pr.st <- c(0, diff(out[, (m + 1)]))
+ plot(time.frame, pr.st, type = "h", main = "", xlab = xlab,
+ ylab = ylab)
+ cp[m] <- which(cumsum(pr.st) > 0.5)[1] - 1
+ abline(v = cp[m] + time.frame[1], lty = 3, col = "red")
+ }
+ cp.means <- rep(NA, m + 1)
cp.start <- c(1, cp + 1)
cp.end <- c(cp, length(y))
cat("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n")
- cat("Expected changepoint(s) ", cp + time.frame[1], '\n')
- for (i in 1:(m+1)) cp.means[i] <- mean(y[cp.start[i]:cp.end[i]])
- cat("Local means for each regime are ", cp.means, '\n')
+ cat("Expected changepoint(s) ", cp + time.frame[1], "\n")
+ for (i in 1:(m + 1)) cp.means[i] <- mean(y[cp.start[i]:cp.end[i]])
+ cat("Local means for each regime are ", cp.means, "\n")
cat("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n")
+}
+ ## prior check for transition matrix
+ "check.P" <-
+ function(P.start = NA, m=m, n=n, a=a, b=b){
+ if (is.na(P.start)[1]){
+ P <- trans.mat.prior(m=m, n=n, a=0.9, b=0.1)}
+ else if ((dim(P.start)[1]==m+1)&(dim(P.start)[2]==m+1)){
+ if ((max(P.start)>1)||(min(P.start)<0)){
+ stop("Error: P starting values are out of unit range.\n")
+ }
+ else
+ P <- P.start
+ }
+ else {
+ stop("Error: P starting values are not conformable.\n")
+ }
+ return(P)
+ }
+
+ ## priro check for mean
+ "check.theta" <-
+ function(theta.start = NA, ns = ns, y = y, min, max){
+ if (is.na(theta.start)[1]){ # draw from uniform with range(y)
+ theta <- runif(ns, min=min, max=max)}
+ else if (length(theta.start)==ns)
+ theta <- theta.start
+ else if (length(theta.start)!=ns) {
+ stop("Error: theta starting values are not conformable.\n")
+ }
+ return(theta)
}
diff --git a/R/zzz.R b/R/zzz.R
index de34ba7..f078b95 100644
--- a/R/zzz.R
+++ b/R/zzz.R
@@ -8,7 +8,7 @@
# echo output to screen
cat("##\n## Markov Chain Monte Carlo Package (MCMCpack)\n")
cat("## Copyright (C) 2003-", this.year,
- " Andrew D. Martin and Kevin M. Quinn\n", sep="")
+ " Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park\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)
diff --git a/README b/README
index abcdb17..d5461ee 100644
--- a/README
+++ b/README
@@ -6,19 +6,17 @@
Andrew D. Martin <admartin at wustl.edu>
Kevin M. Quinn <kevin_quinn at harvard.edu>
+Jong Hee Park <jhp at uchicago.edu>
// Compilation
This package (along with Scythe) uses C++ and the Standard Template
Library (STL). We suggest using of the GCC compiler 4.0 or greater. The
-current package has been tested on GCC 3.1, 3.2, 3.3, 3.4, and 4.0 on Linux
-and MacOS X.
-
-We have also successfully compiled this package for Windows using a
-MinGW cross-compiler using the setup described by Yan and Rossini
-(2006). Many thanks to Dan Pemstein for helping with the
-cross-compilation, and to Kurt Hornik and Fritz Leisch for their help
-with debugging as well as their service to the R community.
+current package has been tested using GCC 4.0 on Linux and MacOS X.
+
+Many thanks to Dan Pemstein for helping with all sorts of C++ issues,
+and to Kurt Hornik and Fritz Leisch for their help with debugging as
+well as their service to the R community.
// Acknowledgments
diff --git a/man/MCMCdynamicEI.Rd b/man/MCMCdynamicEI.Rd
index 02c02d9..04a9a39 100644
--- a/man/MCMCdynamicEI.Rd
+++ b/man/MCMCdynamicEI.Rd
@@ -152,7 +152,7 @@ York: Cambridge University Press.
\url{http://www-fis.iarc.fr/coda/}.
Jonathan C. Wakefield. 2004. ``Ecological Inference for 2 x 2 Tables.'' \emph{Journal
- of the Royal Statistical Society, Series A}. 167(3): 385�445.
+ of the Royal Statistical Society, Series A}. 167(3): 385445.
}
@@ -227,3 +227,5 @@ pairs(cbind(p0.true, p0meanDyn, p0meanHier, p1.true, p1meanDyn, p1meanHier))
\seealso{\code{\link{MCMChierEI}},
\code{\link[coda]{plot.mcmc}},\code{\link[coda]{summary.mcmc}}}
+
+
diff --git a/man/MCMChierEI.Rd b/man/MCMChierEI.Rd
index 8728270..dc1d04b 100644
--- a/man/MCMChierEI.Rd
+++ b/man/MCMChierEI.Rd
@@ -140,7 +140,7 @@ MCMChierEI(r0, r1, c0, c1, burnin=5000, mcmc=50000, thin=1,
\references{
Jonathan C. Wakefield. 2004. ``Ecological Inference for 2 x 2 Tables.'' \emph{Journal
- of the Royal Statistical Society, Series A}. 167(3): 385�445.
+ of the Royal Statistical Society, Series A}. 167(3): 385445.
Radford Neal. 2003. ``Slice Sampling" (with discussion). \emph{Annals of
Statistics}, 31: 705-767.
@@ -187,3 +187,5 @@ pairs(cbind(p0.true, p0meanHier, p1.true, p1meanHier))
\seealso{\code{\link{MCMCdynamicEI}},
\code{\link[coda]{plot.mcmc}},\code{\link[coda]{summary.mcmc}}}
+
+
diff --git a/man/MCMCmnl.Rd b/man/MCMCmnl.Rd
index accb635..13eaee7 100644
--- a/man/MCMCmnl.Rd
+++ b/man/MCMCmnl.Rd
@@ -14,8 +14,8 @@
\usage{
MCMCmnl(formula, baseline = NULL, data = parent.frame(),
burnin = 1000, mcmc = 10000, thin = 1,
- mcmc.method = "MH", tune = 1.1, verbose = 0,
- seed = NA, beta.start = NA, b0 = 0, B0 = 0, ...) }
+ mcmc.method = c("IndMH", "RWM", "slice"), tune = 1, tdf=6,
+ verbose = 0, seed = NA, beta.start = NA, b0 = 0, B0 = 0, ...) }
\arguments{
\item{formula}{Model formula.
@@ -66,10 +66,15 @@ MCMCmnl(formula, baseline = NULL, data = parent.frame(),
\item{thin}{The thinning interval used in the simulation. The number of
mcmc iterations must be divisible by this value. }
- \item{mcmc.method}{Can be set to either "MH" (default) or "slice" to
- perform random walk Metropolis sampling or slice sampling
+ \item{mcmc.method}{Can be set to either "IndMH" (default), "RWM",
+ or "slice" to perform independent Metropolis-Hastings sampling,
+ random walk Metropolis sampling or slice sampling
respectively.}
+ \item{tdf}{Degrees of freedom for the multivariate-t proposal
+ distribution when \code{mcmc.method} is set to "IndMH". Must be
+ positive. }
+
\item{tune}{Metropolis tuning parameter. Can be either a positive
scalar or a \eqn{k}{k}-vector, where \eqn{k}{k} is the length of
\eqn{\beta}{beta}. Make sure that the
diff --git a/man/MCMCoprobit.Rd b/man/MCMCoprobit.Rd
index 02ad1b6..faceaab 100644
--- a/man/MCMCoprobit.Rd
+++ b/man/MCMCoprobit.Rd
@@ -12,8 +12,8 @@
\usage{
MCMCoprobit(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
- thin=1, tune = NA, verbose = 0, seed = NA, beta.start = NA,
- b0 = 0, B0 = 0, ...) }
+ thin=1, tune = NA, tdf = 1, verbose = 0, seed = NA, beta.start = NA,
+ b0 = 0, B0 = 0, a0 = 0, A0 = 0, mcmc.method = c("Cowles", "AC"), ...) }
\arguments{
\item{formula}{Model formula.}
@@ -31,6 +31,10 @@ MCMCoprobit(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
step. Default of NA corresponds to a choice of 0.05 divided by the
number of categories in the response variable.}
+ \item{tdf}{Degrees of freedom for the multivariate-t proposal
+ distribution when \code{mcmc.method} is set to "IndMH". Must be
+ positive. }
+
\item{verbose}{A switch which determines whether or not the progress of
the sampler is printed to the screen. If \code{verbose} is greater
than 0 the iteration
@@ -65,6 +69,21 @@ MCMCoprobit(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
identity matrix serves as the prior precision of \eqn{\beta}{beta}.
Default value of 0 is equivalent to an improper uniform prior on
\eqn{\beta}{beta}. }
+
+ \item{a0}{The prior mean of \eqn{\alpha}{alpha}. This can either be a
+ scalar or a column vector with dimension equal to the number of
+ betas. If this takes a scalar value, then that value will serve as
+ the prior mean for all of the betas.}
+
+ \item{A0}{The prior precision of \eqn{\alpha}{alpha}. This can either be a
+ scalar or a square matrix with dimensions equal to the number of
+ betas. If this takes a scalar value, then that value times an
+ identity matrix serves as the prior precision of \eqn{\beta}{beta}.
+ Default value of 0 is equivalent to an improper uniform prior on
+ \eqn{\beta}{beta}. }
+
+ \item{mcmc.method}{Can be set to either "Cowles" (default) or "AC" to perform
+posterior sampling of cutpoints based on Cowles (1996) or Albert and Chib (2001) respectively.}
\item{...}{further arguments to be passed}
}
@@ -98,13 +117,12 @@ that can be used to analyze the posterior sample.
These probabilities are used to form the multinomial distribution
that defines the likelihoods.
- The algorithm employed is discussed in depth by Cowles (1996). Note that
- the model does include a constant in the data matrix. Thus, the first
- element \eqn{\gamma_1}{gamma_1} is normalized to zero, and is not
- returned in the mcmc object.
+ \code{MCMCoprobit} provides two ways to sample the cutpoints. Cowles (1996) proposes a sampling scheme that groups sampling of a latent variable with cutpoints. In this case, for identification the first element \eqn{\gamma_1}{gamma_1} is normalized to zero. Albert and Chib (2001) show that we can sample cutpoints indirectly without constraints by transforming cutpoints into real-valued parameters (\eqn{\alpha}{alpha}).
}
\references{
+ Albert, James and Siddhartha Chib. 2001. ``Sequential Ordinal Modeling with Applications to Survival Data." \emph{Biometrics.} 57: 829-836.
+
M. K. Cowles. 1996. ``Accelerating Monte Carlo Markov Chain Convergence for
Cumulative-link Generalized Linear Models." \emph{Statistics and Computing.}
6: 101-110.
@@ -127,9 +145,12 @@ that can be used to analyze the posterior sample.
z <- 1.0 + x1*0.1 - x2*0.5 + rnorm(100);
y <- z; y[z < 0] <- 0; y[z >= 0 & z < 1] <- 1;
y[z >= 1 & z < 1.5] <- 2; y[z >= 1.5] <- 3;
- posterior <- MCMCoprobit(y ~ x1 + x2, tune=0.3, mcmc=20000)
- plot(posterior)
- summary(posterior)
+ out1 <- MCMCoprobit(y ~ x1 + x2, tune=0.3)
+ out2 <- MCMCoprobit(y ~ x1 + x2, tune=0.3, tdf=3, verbose=1000, mcmc.method="AC")
+ summary(out1)
+ summary(out2)
+ plot(out1)
+ plot(out2)
}
}
diff --git a/man/MCMCpoissonChangepoint.Rd b/man/MCMCpoissonChangepoint.Rd
old mode 100644
new mode 100755
index 7a63728..8e1c4f3
--- a/man/MCMCpoissonChangepoint.Rd
+++ b/man/MCMCpoissonChangepoint.Rd
@@ -1,137 +1,124 @@
-\name{MCMCpoissonChangepoint}
-\alias{MCMCpoissonChangepoint}
-
-\title{Markov Chain Monte Carlo for a Poisson Multiple Changepoint Model}
-\description{
- This function generates a sample from the posterior distribution
- of a Poisson model with multiple changepoints. The function uses
- the Markov chain Monte Carlo method of Chib (1998).
- The user supplies data and priors, and
- a sample from the posterior distribution is returned as an mcmc
- object, which can be subsequently analyzed with functions
- provided in the coda package.
-}
-
-\usage{MCMCpoissonChangepoint(data, m = 1, burnin = 1000, mcmc = 1000,
- thin = 1, verbose = 0, seed = NA, c0, d0, a = NULL, b = NULL,
- marginal.likelihood = c("none", "Chib95"), ...) }
-
-\arguments{
- \item{data}{The data.}
-
- \item{m}{The number of changepoints.}
-
- \item{burnin}{The number of burn-in iterations for the sampler.}
-
- \item{mcmc}{The number of MCMC iterations after burn-in.}
-
- \item{thin}{The thinning interval used in the simulation. The number of
- MCMC iterations must be divisible by this value.}
-
- \item{verbose}{A switch which determines whether or not the progress of
- the sampler is printed to the screen. If \code{verbose} is greater
- than 0, the iteration number and the posterior density samples are printed to the screen every \code{verbose}th iteration.}
-
- \item{seed}{The seed for the random number generator. If NA, current R
- system seed is used.}
-
- \item{c0}{\eqn{c_0}{c0} is the shape parameter for Gamma prior on \eqn{\lambda}{lambda}
- (the mean).}
-
- \item{d0}{\eqn{d_0}{d0} is the scale parameter for Gamma prior on \eqn{\lambda}{lambda}
- (the mean).}
-
- \item{a}{\eqn{a}{a} is the shape1 beta prior for transition probabilities. By default,
- the expected duration is computed and corresponding a and b values are assigned. The expected
- duration is the sample period divided by the number of states.}
-
- \item{b}{\eqn{b}{b} is the shape2 beta prior for transition probabilities. By default,
- the expected duration is computed and corresponding a and b values are assigned. The expected
- duration is the sample period divided by the number of states.}
-
- \item{marginal.likelihood}{How should the marginal likelihood be
- calculated? Options are: \code{none} in which case the marginal
- likelihood will not be calculated, and
- \code{Chib95} in which case the method of Chib (1995) is used.}
-
- \item{...}{further arguments to be passed}
-}
-
-\value{
- An mcmc object that contains the posterior sample. This
- object can be summarized by functions provided by the coda package. The object contains an attribute \code{prob.state} storage matrix that contains the probability of \eqn{state_i}{state_i} for each period, and the
- log-likelihood of the model (\code{log.like}).
-}
-
-\details{
- \code{MCMCpoissonChangepoint} simulates from the posterior distribution of
- a Poisson model with multiple changepoints.
-
- The model takes the following form:
- \deqn{Y_t \sim \mathcal{P}oisson(\lambda_i),\;\; i = 1, \ldots, k}{Y_t ~ Poisson(lambda_i), i = 1,...,k.}
- Where \eqn{k}{k} is the number of states.
-
- We assume Gamma priors for \eqn{\lambda_{i}}{lambda_i} and Beta priors for transition probabilities:
- \deqn{\lambda_i \sim \mathcal{G}amma(c_0, d_0)}{lambda_i ~ Gamma(c0, d0)}
- \deqn{p_{ii} \sim \mathcal{B}eta{a}{b},\;\; i = 1, \ldots, k}{p_ii ~ Beta(a, b), i = 1,...,k.}
- Where \eqn{k}{k} is the number of states.
-
- Note that no default value is provided for Gamma priors.
-
- The simulation in this model-fitting is performed in R.
- }
-
-\author{Jong Hee Park, \email{jhp at uchicago.edu}, \url{http://home.uchicago.edu/~jhp/}.}
-
-\references{
- Siddhartha Chib. 1995. "Marginal Likelihood from the Gibbs Output."
- \emph{Journal of the American Statistical Association}. 90:
- 1313-1321.
-
- Siddhartha Chib. 1998. "Estimation and comparison of multiple change-point models."
- \emph{Journal of Econometrics}. 86: 221-241.
-
-
- Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2002.
- \emph{Output Analysis and Diagnostics for MCMC (CODA)}.
- \url{http://www-fis.iarc.fr/coda/}.
-}
-
-\examples{
- \dontrun{
- ## generate event count time series data with two breakpoints at 100 and 200
- set.seed(1973)
- n <- 300
- true.lambda <- c(3, 2, 4)
- y1 <- rpois(100, true.lambda[1])
- y2 <- rpois(100, true.lambda[2])
- y3 <- rpois(100, true.lambda[3])
- y <- c(y1, y2, y3)
-
- ## run the example
- model1 <- MCMCpoissonChangepoint(y, m=1, burnin = 100, mcmc = 100,
- thin=1, verbose=5, c0=3, d0=1, marginal.likelihood =
- c("Chib95"))
- model2 <- MCMCpoissonChangepoint(y, m=2, burnin = 100, mcmc = 100,
- thin=1, verbose=5, c0=3, d0=1, marginal.likelihood =
- c("Chib95"))
- model3 <- MCMCpoissonChangepoint(y, m=3, burnin = 100, mcmc = 100,
- thin=1, verbose=5, c0=3, d0=1, marginal.likelihood =
- c("Chib95"))
- model4 <- MCMCpoissonChangepoint(y, m=4, burnin = 100, mcmc = 100,
- thin=1, verbose=5, c0=3, d0=1, marginal.likelihood =
- c("Chib95"))
-
- ## perform model comparison
- BayesFactor(model1, model2, model3, model4)
-
- ## Draw plots using the "right" model
- plotPostState(model2, y, m=2, legend.control=NULL)
- plotPostChangepoint(model2, y, m=2)
- }
-}
-
-\keyword{models}
-
-\seealso{\code{\link{plotPostState}}, \code{\link{plotPostChangepoint}}}
-
+\name{MCMCpoissonChangepoint}
+\alias{MCMCpoissonChangepoint}
+
+\title{Markov Chain Monte Carlo for a Poisson Multiple Changepoint Model}
+\description{
+ This function generates a sample from the posterior distribution
+ of a Poisson model with multiple changepoints. The function uses
+ the Markov chain Monte Carlo method of Chib (1998).
+ The user supplies data and priors, and
+ a sample from the posterior distribution is returned as an mcmc
+ object, which can be subsequently analyzed with functions
+ provided in the coda package.
+}
+
+\usage{MCMCpoissonChangepoint(data, m = 1, c0 = NA, d0 = NA, a = NULL, b = NULL,
+ burnin = 10000, mcmc = 10000, thin = 1, verbose = 0, seed = NA,
+ lambda.start = NA, P.start = NA,
+ marginal.likelihood = c("none", "Chib95"), ...) }
+
+\arguments{
+ \item{data}{The data.}
+
+ \item{m}{The number of changepoints.}
+
+ \item{c0}{\eqn{c_0}{c0} is the shape parameter for Gamma prior on \eqn{\lambda}{lambda}
+ (the mean). No default value is provided.}
+
+ \item{d0}{\eqn{d_0}{d0} is the scale parameter for Gamma prior on \eqn{\lambda}{lambda}
+ (the mean). No default value is provided.}
+
+ \item{a}{\eqn{a}{a} is the shape1 beta prior for transition probabilities. By default,
+ the expected duration is computed and corresponding a and b values are assigned. The expected
+ duration is the sample period divided by the number of states.}
+
+ \item{b}{\eqn{b}{b} is the shape2 beta prior for transition probabilities. By default,
+ the expected duration is computed and corresponding a and b values are assigned. The expected
+ duration is the sample period divided by the number of states.}
+
+ \item{burnin}{The number of burn-in iterations for the sampler.}
+
+ \item{mcmc}{The number of MCMC iterations after burn-in.}
+
+ \item{thin}{The thinning interval used in the simulation. The number of
+ MCMC iterations must be divisible by this value.}
+
+ \item{verbose}{A switch which determines whether or not the progress of
+ the sampler is printed to the screen. If \code{verbose} is greater
+ than 0, the iteration number and the posterior density samples are printed to the screen every \code{verbose}th iteration.}
+
+ \item{seed}{The seed for the random number generator. If NA, current R
+ system seed is used.}
+
+ \item{lambda.start}{The starting values for the lambda vector. This can either be a scalar or a column vector with dimension equal to the number of lambdas. The default value of NA will use draws from the Uniform distribution with the same boundary with the data as the starting value. If this is a scalar, that value will serve as the starting value mean for all of the lambdas.}
+
+ \item{P.start}{The starting values for the transition matrix. A user should provide a square matrix with dimension equal to the number of states. By default, draws from the \code{Beta(0.9, 0.1)} are used to construct a proper transition matrix for each raw except the last raw.}
+
+ \item{marginal.likelihood}{How should the marginal likelihood be
+ calculated? Options are: \code{none} in which case the marginal
+ likelihood will not be calculated, and
+ \code{Chib95} in which case the method of Chib (1995) is used.}
+
+ \item{...}{further arguments to be passed}
+}
+
+\value{
+ An mcmc object that contains the posterior sample. This
+ object can be summarized by functions provided by the coda package.
+ The object contains an attribute \code{prob.state} storage matrix that contains the probability of \eqn{state_i}{state_i} for each period, and
+ the log-marginal likelihood of the model (\code{logmarglike}).
+}
+
+\details{
+ \code{MCMCpoissonChangepoint} simulates from the posterior distribution of
+ a Poisson model with multiple changepoints.
+
+ The model takes the following form:
+ \deqn{Y_t \sim \mathcal{P}oisson(\lambda_i),\;\; i = 1, \ldots, k}{Y_t ~ Poisson(lambda_i), i = 1,...,k.}
+ Where \eqn{k}{k} is the number of states.
+
+ We assume Gamma priors for \eqn{\lambda_{i}}{lambda_i} and Beta priors for transition probabilities:
+ \deqn{\lambda_i \sim \mathcal{G}amma(c_0, d_0)}{lambda_i ~ Gamma(c0, d0)}
+ \deqn{p_{ii} \sim \mathcal{B}eta{a}{b},\;\; i = 1, \ldots, k}{p_ii ~ Beta(a, b), i = 1,...,k.}
+ Where \eqn{k}{k} is the number of states.
+
+ The simulation is done in compiled C++ code to maximize efficiency.
+ }
+
+\author{Jong Hee Park, \email{jhp at uchicago.edu}, \url{http://home.uchicago.edu/~jhp/}.}
+
+\references{
+ Siddhartha Chib. 1995. "Marginal Likelihood from the Gibbs Output."
+ \emph{Journal of the American Statistical Association}. 90:
+ 1313-1321.
+
+ Siddhartha Chib. 1998. "Estimation and comparison of multiple change-point models."
+ \emph{Journal of Econometrics}. 86: 221-241.
+}
+
+\examples{
+ \dontrun{
+ set.seed(1973)
+ n <- 100
+ true.lambda <- c(5, 8)
+ y1 <- rpois(n/2, true.lambda[1])
+ y2 <- rpois(n/2, true.lambda[2])
+ y <- c(y1, y2)
+
+ model1 <- MCMCpoissonChangepoint(y, m=1, c0=6.85, d0=1, verbose = 10000, marginal.likelihood = "Chib95")
+ model2 <- MCMCpoissonChangepoint(y, m=2, c0=6.85, d0=1, verbose = 10000, marginal.likelihood = "Chib95")
+ model3 <- MCMCpoissonChangepoint(y, m=3, c0=6.85, d0=1, verbose = 10000, marginal.likelihood = "Chib95")
+ model4 <- MCMCpoissonChangepoint(y, m=4, c0=6.85, d0=1, verbose = 10000, marginal.likelihood = "Chib95")
+ model5 <- MCMCpoissonChangepoint(y, m=5, c0=6.85, d0=1, verbose = 10000, marginal.likelihood = "Chib95")
+
+ print(BayesFactor(model1, model2, model3, model4, model5))
+
+ ## Draw plots using the "right" model
+ plotState(model1)
+ plotChangepoint(model1)
+ }
+}
+
+\keyword{models}
+
+\seealso{\code{\link{plotState}}, \code{\link{plotChangepoint}}}
diff --git a/man/MCMCtobit.Rd b/man/MCMCtobit.Rd
index b0ceea2..0db4277 100644
--- a/man/MCMCtobit.Rd
+++ b/man/MCMCtobit.Rd
@@ -1,160 +1,160 @@
-\name{MCMCtobit}
-\alias{MCMCtobit}
-\title{Markov Chain Monte Carlo for Gaussian Linear Regression with a Censored Dependent Variable}
-\description{
- This function generates a sample from the posterior distribution
- of a linear regression model with Gaussian errors using
- Gibbs sampling (with a multivariate Gaussian prior on the
- beta vector, and an inverse Gamma prior on the conditional
- error variance). The dependent variable may be censored
- from below, from above, or both. The user supplies data
- and priors, and a sample from the posterior distribution is
- returned as an mcmc object, which can be subsequently
- analyzed with functions provided in the coda package.
- }
-\usage{
-MCMCtobit(formula, data = parent.frame(), below = 0, above = Inf,
- burnin = 1000, mcmc = 10000, thin = 1, verbose = 0, seed = NA,
- beta.start = NA, b0 = 0, B0 = 0, c0 = 0.001, d0 = 0.001, ...)
-}
-\arguments{
- \item{formula}{A model formula.}
-
- \item{data}{A dataframe.}
-
- \item{below}{The point at which the dependent variable is
- censored from below. The default is zero. To censor from
- above only, specify that below = -Inf.}
-
- \item{above}{The point at which the dependent variable is
- censored from above. To censor from below only, use the
- default value of Inf.}
-
- \item{burnin}{The number of burn-in iterations for the sampler.}
-
- \item{mcmc}{The number of MCMC iterations after burnin.}
-
- \item{thin}{The thinning interval used in the simulation. The number of
- MCMC iterations must be divisible by this value.}
-
- \item{verbose}{A switch which determines whether or not the progress of
- the sampler is printed to the screen. If \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
- passed it is used to seed the Mersenne twister. The user can also
- pass a list of length two to use the L'Ecuyer random number generator,
- which is suitable for parallel computation. The first element of the
- list is the L'Ecuyer seed, which is a vector of length six or NA (if NA
- a default seed of \code{rep(12345,6)} is used). The second element of
- list is a positive substream number. See the MCMCpack specification
- for more details.}
-
- \item{beta.start}{The starting values for the \eqn{\beta}{beta} vector.
- This can either be a scalar or a column vector with dimension equal to the
- number of betas. The default value of of NA will use the OLS estimate of
- \eqn{\beta}{beta} as the starting value. If this is a scalar, that value
- will serve as the starting value mean for all of the betas.}
-
- \item{b0}{The prior mean of \eqn{\beta}{beta}. This can either be a
- scalar or a column vector with dimension equal to the number of betas. If
- this takes a scalar value, then that value will serve as the prior
- mean for all of the betas.}
-
- \item{B0}{The prior precision of \eqn{\beta}{beta}. This can either be a
- scalar or a square matrix with dimensions equal to the number of
- betas.
- If this takes a scalar value, then that value times an
- identity matrix serves as the prior precision of beta. Default value
- of 0 is equivalent to an improper uniform prior for beta.}
-
- \item{c0}{\eqn{c_0/2}{c0/2} is the shape parameter for the inverse
- Gamma prior on \eqn{\sigma^2}{sigma^2} (the variance of the
- disturbances). The amount of information in the inverse Gamma prior
- is something like that from \eqn{c_0}{c0} pseudo-observations.}
-
- \item{d0}{\eqn{d_0/2}{d0/2} is the scale parameter for the
- inverse Gamma prior on \eqn{\sigma^2}{sigma^2} (the variance of the
- disturbances). In constructing the inverse Gamma prior,
- \eqn{d_0}{d0} acts like the sum of squared errors from the
- \eqn{c_0}{c0} pseudo-observations.}
-
- \item{...}{further arguments to be passed}
-
-}
-
-\details{
- \code{MCMCtobit} simulates from the posterior distribution using standard
- Gibbs sampling (a multivariate Normal draw for the betas, and an
- inverse Gamma draw for the conditional error variance). \code{MCMCtobit}
- differs from \code{MCMCregress} in that the dependent variable may be
- censored from below, from above, or both. The simulation
- proper is done in compiled C++ code to maximize efficiency. Please consult
- the coda documentation for a comprehensive list of functions that can be
- used to analyze the posterior sample.
-
- The model takes the following form:
- \deqn{y_i = x_i ' \beta + \varepsilon_{i},}{y_i = x_i'beta + epsilon_i,}
- where the errors are assumed to be Gaussian:
- \deqn{\varepsilon_{i} \sim \mathcal{N}(0, \sigma^2).}{epsilon_i ~ N(0,
- sigma^2).}
- Let \eqn{c_1} and \eqn{c_2} be the two censoring points, and let
- \eqn{y_i^\ast}{y_i^star} be the partially observed dependent variable. Then,
- \deqn{y_i = y_i^{\ast} \texttt{ if } c_1 < y_i^{\ast} < c_2,}{y_i = y_i^star if c_1 < y_i^star < c_2,}
- \deqn{y_i = c_1 \texttt{ if } c_1 \geq y_i^{\ast},}{y_i = c_1 if c_1 >= y_i^star,}
- \deqn{y_i = c_2 \texttt{ if } c_2 \leq y_i^{\ast}.}{y_i = c_2 if c_1 <= y_i^star.}
- We assume standard, semi-conjugate priors:
- \deqn{\beta \sim \mathcal{N}(b_0,B_0^{-1}),}{beta ~ N(b0,B0^(-1)),}
- and:
- \deqn{\sigma^{-2} \sim \mathcal{G}amma(c_0/2, d_0/2),}{sigma^(-2) ~
- Gamma(c0/2, d0/2),}
- where \eqn{\beta}{beta} and \eqn{\sigma^{-2}}{sigma^(-2)} are assumed
- \emph{a priori} independent. Note that only starting values for
- \eqn{\beta}{beta} are allowed because simulation is done using
- Gibbs sampling with the conditional error variance
- as the first block in the sampler.
- }
-
-\value{
- An mcmc object that contains the posterior sample. This
- object can be summarized by functions provided by the coda package.
-}
-
-\references{
- 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/}.
-
- James Tobin. 1958. ``Estimation of relationships for limited dependent
- variables." \emph{Econometrica.} 26:24-36.
-}
-
-\author{Ben Goodrich, \email{goodrich.ben at gmail.com},
- \url{http://www.people.fas.harvard.edu/~goodrich/}}
-
-\seealso{
- \code{\link[coda]{plot.mcmc}},
- \code{\link[coda]{summary.mcmc}},
- \code{\link[survival]{survreg}},
- \code{\link[MCMCpack]{MCMCregress}}}
-
+\name{MCMCtobit}
+\alias{MCMCtobit}
+\title{Markov Chain Monte Carlo for Gaussian Linear Regression with a Censored Dependent Variable}
+\description{
+ This function generates a sample from the posterior distribution
+ of a linear regression model with Gaussian errors using
+ Gibbs sampling (with a multivariate Gaussian prior on the
+ beta vector, and an inverse Gamma prior on the conditional
+ error variance). The dependent variable may be censored
+ from below, from above, or both. The user supplies data
+ and priors, and a sample from the posterior distribution is
+ returned as an mcmc object, which can be subsequently
+ analyzed with functions provided in the coda package.
+ }
+\usage{
+MCMCtobit(formula, data = parent.frame(), below = 0, above = Inf,
+ burnin = 1000, mcmc = 10000, thin = 1, verbose = 0, seed = NA,
+ beta.start = NA, b0 = 0, B0 = 0, c0 = 0.001, d0 = 0.001, ...)
+}
+\arguments{
+ \item{formula}{A model formula.}
+
+ \item{data}{A dataframe.}
+
+ \item{below}{The point at which the dependent variable is
+ censored from below. The default is zero. To censor from
+ above only, specify that below = -Inf.}
+
+ \item{above}{The point at which the dependent variable is
+ censored from above. To censor from below only, use the
+ default value of Inf.}
+
+ \item{burnin}{The number of burn-in iterations for the sampler.}
+
+ \item{mcmc}{The number of MCMC iterations after burnin.}
+
+ \item{thin}{The thinning interval used in the simulation. The number of
+ MCMC iterations must be divisible by this value.}
+
+ \item{verbose}{A switch which determines whether or not the progress of
+ the sampler is printed to the screen. If \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
+ passed it is used to seed the Mersenne twister. The user can also
+ pass a list of length two to use the L'Ecuyer random number generator,
+ which is suitable for parallel computation. The first element of the
+ list is the L'Ecuyer seed, which is a vector of length six or NA (if NA
+ a default seed of \code{rep(12345,6)} is used). The second element of
+ list is a positive substream number. See the MCMCpack specification
+ for more details.}
+
+ \item{beta.start}{The starting values for the \eqn{\beta}{beta} vector.
+ This can either be a scalar or a column vector with dimension equal to the
+ number of betas. The default value of of NA will use the OLS estimate of
+ \eqn{\beta}{beta} as the starting value. If this is a scalar, that value
+ will serve as the starting value mean for all of the betas.}
+
+ \item{b0}{The prior mean of \eqn{\beta}{beta}. This can either be a
+ scalar or a column vector with dimension equal to the number of betas. If
+ this takes a scalar value, then that value will serve as the prior
+ mean for all of the betas.}
+
+ \item{B0}{The prior precision of \eqn{\beta}{beta}. This can either be a
+ scalar or a square matrix with dimensions equal to the number of
+ betas.
+ If this takes a scalar value, then that value times an
+ identity matrix serves as the prior precision of beta. Default value
+ of 0 is equivalent to an improper uniform prior for beta.}
+
+ \item{c0}{\eqn{c_0/2}{c0/2} is the shape parameter for the inverse
+ Gamma prior on \eqn{\sigma^2}{sigma^2} (the variance of the
+ disturbances). The amount of information in the inverse Gamma prior
+ is something like that from \eqn{c_0}{c0} pseudo-observations.}
+
+ \item{d0}{\eqn{d_0/2}{d0/2} is the scale parameter for the
+ inverse Gamma prior on \eqn{\sigma^2}{sigma^2} (the variance of the
+ disturbances). In constructing the inverse Gamma prior,
+ \eqn{d_0}{d0} acts like the sum of squared errors from the
+ \eqn{c_0}{c0} pseudo-observations.}
+
+ \item{...}{further arguments to be passed}
+
+}
+
+\details{
+ \code{MCMCtobit} simulates from the posterior distribution using standard
+ Gibbs sampling (a multivariate Normal draw for the betas, and an
+ inverse Gamma draw for the conditional error variance). \code{MCMCtobit}
+ differs from \code{MCMCregress} in that the dependent variable may be
+ censored from below, from above, or both. The simulation
+ proper is done in compiled C++ code to maximize efficiency. Please consult
+ the coda documentation for a comprehensive list of functions that can be
+ used to analyze the posterior sample.
+
+ The model takes the following form:
+ \deqn{y_i = x_i ' \beta + \varepsilon_{i},}{y_i = x_i'beta + epsilon_i,}
+ where the errors are assumed to be Gaussian:
+ \deqn{\varepsilon_{i} \sim \mathcal{N}(0, \sigma^2).}{epsilon_i ~ N(0,
+ sigma^2).}
+ Let \eqn{c_1} and \eqn{c_2} be the two censoring points, and let
+ \eqn{y_i^\ast}{y_i^star} be the partially observed dependent variable. Then,
+ \deqn{y_i = y_i^{\ast} \texttt{ if } c_1 < y_i^{\ast} < c_2,}{y_i = y_i^star if c_1 < y_i^star < c_2,}
+ \deqn{y_i = c_1 \texttt{ if } c_1 \geq y_i^{\ast},}{y_i = c_1 if c_1 >= y_i^star,}
+ \deqn{y_i = c_2 \texttt{ if } c_2 \leq y_i^{\ast}.}{y_i = c_2 if c_1 <= y_i^star.}
+ We assume standard, semi-conjugate priors:
+ \deqn{\beta \sim \mathcal{N}(b_0,B_0^{-1}),}{beta ~ N(b0,B0^(-1)),}
+ and:
+ \deqn{\sigma^{-2} \sim \mathcal{G}amma(c_0/2, d_0/2),}{sigma^(-2) ~
+ Gamma(c0/2, d0/2),}
+ where \eqn{\beta}{beta} and \eqn{\sigma^{-2}}{sigma^(-2)} are assumed
+ \emph{a priori} independent. Note that only starting values for
+ \eqn{\beta}{beta} are allowed because simulation is done using
+ Gibbs sampling with the conditional error variance
+ as the first block in the sampler.
+ }
+
+\value{
+ An mcmc object that contains the posterior sample. This
+ object can be summarized by functions provided by the coda package.
+}
+
+\references{
+ 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/}.
+
+ James Tobin. 1958. ``Estimation of relationships for limited dependent
+ variables." \emph{Econometrica.} 26:24-36.
+}
+
+\author{Ben Goodrich, \email{goodrich.ben at gmail.com},
+ \url{http://www.people.fas.harvard.edu/~goodrich/}}
+
+\seealso{
+ \code{\link[coda]{plot.mcmc}},
+ \code{\link[coda]{summary.mcmc}},
+ \code{\link[survival]{survreg}},
+ \code{\link[MCMCpack]{MCMCregress}}}
+
\examples{
-\dontrun{
-library(survival)
-example(tobin)
-summary(tfit)
-tfit.mcmc <- MCMCtobit(durable ~ age + quant, data=tobin, mcmc=30000,
- verbose=1000)
-plot(tfit.mcmc)
-raftery.diag(tfit.mcmc)
-summary(tfit.mcmc)
+\dontrun{
+library(survival)
+example(tobin)
+summary(tfit)
+tfit.mcmc <- MCMCtobit(durable ~ age + quant, data=tobin, mcmc=30000,
+ verbose=1000)
+plot(tfit.mcmc)
+raftery.diag(tfit.mcmc)
+summary(tfit.mcmc)
+}
}
-}
-
-\keyword{models}
+
+\keyword{models}
diff --git a/man/dtomog.Rd b/man/dtomog.Rd
index 5454a27..4fe99e9 100644
--- a/man/dtomog.Rd
+++ b/man/dtomog.Rd
@@ -75,7 +75,7 @@ dtomogplot(r0, r1, c0, c1, time.vec=NA, delay=0,
Princeton: Princeton University Press.
Jonathan C. Wakefield. 2004. ``Ecological Inference for 2 x 2 Tables.''
- \emph{Journal of the Royal Statistical Society, Series A}. 167(3): 385�445.
+ \emph{Journal of the Royal Statistical Society, Series A}. 167(3): 385445.
Kevin Quinn. 2004. ``Ecological Inference in the Presence of Temporal
Dependence." In \emph{Ecological Inference: New Methodological
@@ -121,3 +121,4 @@ dtomogplot(r0, r1, c0, c1, delay=0.1)
\seealso{\code{\link{MCMChierEI}},
\code{\link{MCMCdynamicEI}},\code{\link{tomogplot}}
}
+
diff --git a/man/plotPostChangepoint.Rd b/man/plotChangepoint.Rd
old mode 100644
new mode 100755
similarity index 75%
rename from man/plotPostChangepoint.Rd
rename to man/plotChangepoint.Rd
index 9ba2ffd..6bd6c7b
--- a/man/plotPostChangepoint.Rd
+++ b/man/plotChangepoint.Rd
@@ -1,20 +1,16 @@
-\name{plotPostChangepoint}
-\alias{plotPostChangepoint}
+\name{plotChangepoint}
+\alias{plotChangepoint}
\title{Changepoint Location Plots}
\description{Plot the posterior distribution over the location of the changepoints.}
\usage{
- plotPostChangepoint(mcmcout, y, m, xlab="Time", ylab="")
+ plotChangepoint(mcmcout, xlab="Time", ylab="")
}
\arguments{
\item{mcmcout}{The \code{mcmc} object containing the posterior density sample from a changepoint model. Note that this must have a \code{prob.state} attribute.}
-\item{y}{The data.}
-
-\item{m}{The number of changepoints.}
-
\item{xlab}{Label for the x-axis.}
\item{ylab}{Label for the y-axis.}
diff --git a/man/plotPostState.Rd b/man/plotState.Rd
old mode 100644
new mode 100755
similarity index 79%
rename from man/plotPostState.Rd
rename to man/plotState.Rd
index 5dc078d..1d21369
--- a/man/plotPostState.Rd
+++ b/man/plotState.Rd
@@ -1,20 +1,16 @@
-\name{plotPostState}
-\alias{plotPostState}
+\name{plotState}
+\alias{plotState}
\title{Changepoint State Plot}
\description{Plot the posterior probability that each time point is in each state.}
\usage{
- plotPostState(mcmcout, y, m, legend.control=NULL, cex=0.8, lwd=1.2)
+ plotState(mcmcout, legend.control=NULL, cex=0.8, lwd=1.2)
}
\arguments{
\item{mcmcout}{The \code{mcmc} object containing the posterior density sample from a changepoint model. Note that this must have a \code{prob.state} attribute.}
-\item{y}{The data.}
-
-\item{m}{The number of changepoints.}
-
\item{legend.control}{Control the location of the legend. It is necessary to pass both the x and y locations; i.e., \code{c(x,y)}.}
\item{cex}{Control point size.}
diff --git a/man/tomog.Rd b/man/tomog.Rd
index d3b1725..9a7ad03 100644
--- a/man/tomog.Rd
+++ b/man/tomog.Rd
@@ -62,7 +62,7 @@ tomogplot(r0, r1, c0, c1, xlab="fraction of r0 in c0 (p0)",
Princeton: Princeton University Press.
Jonathan C. Wakefield. 2004. ``Ecological Inference for 2 x 2 Tables.''
- \emph{Journal of the Royal Statistical Society, Series A}. 167(3): 385�445.
+ \emph{Journal of the Royal Statistical Society, Series A}. 167(3): 385445.
}
\examples{
@@ -78,3 +78,5 @@ tomogplot(r0, r1, c0, c1)
}
+
+
diff --git a/src/MCMCmnl.h b/src/MCMCmnl.h
index 44cb79a..d5a2ec5 100644
--- a/src/MCMCmnl.h
+++ b/src/MCMCmnl.h
@@ -34,13 +34,13 @@ using namespace scythe;
inline double
mnl_logpost(const Matrix<>& Y, const Matrix<>& X, const Matrix<>& beta,
- const Matrix<>& beta_prior_mean,
- const Matrix<>& beta_prior_prec)
+ const Matrix<>& beta_prior_mean,
+ const Matrix<>& beta_prior_prec)
{
-
+
// likelihood
double loglike = 0.0;
- Matrix<double,Row> numera = exp(X * beta);
+ const Matrix<double,Row> numera = exp(X * beta);
//numer = reshape(numer, Y.rows(), Y.cols());
//numer.resize(Y.rows(), Y.cols(), true);
Matrix<double,Row> numer(Y.rows(), Y.cols(), false);
@@ -50,23 +50,30 @@ mnl_logpost(const Matrix<>& Y, const Matrix<>& X, const Matrix<>& beta,
denom[i] = 0.0;
for (unsigned int j = 0; j < Y.cols(); ++j) {
if (Y(i,j) != -999){
- denom[i] += numer(i,j);
+ denom[i] += numer(i,j);
}
}
for (unsigned int j = 0; j < Y.cols(); ++j) {
if (Y(i,j) == 1.0){
- loglike += std::log(numer(i,j) / denom[i]);
+ loglike += std::log(numer(i,j) / denom[i]);
}
}
}
delete [] denom;
-
+
// prior
- double logprior = 0.0;
- if (beta_prior_prec(0,0) != 0) {
- logprior = lndmvn(beta, beta_prior_mean, invpd(beta_prior_prec));
- }
+ // double logprior = 0.0;
+ //if (beta_prior_prec(0,0) != 0) {
+ // logprior = lndmvn(beta, beta_prior_mean, invpd(beta_prior_prec));
+ // }
+ //
+ // the following is only up to proportionality
+ const double logprior = -0.5 *(t(beta - beta_prior_mean) *
+ beta_prior_prec *
+ (beta - beta_prior_mean))(0);
+
+
return (loglike + logprior);
}
diff --git a/src/MCMCmnlMH.cc b/src/MCMCmnlMH.cc
index 0f7a984..b4bf93e 100644
--- a/src/MCMCmnlMH.cc
+++ b/src/MCMCmnlMH.cc
@@ -44,113 +44,190 @@
using namespace std;
using namespace scythe;
-template <typename RNGTYPE>
-void MCMCmnlMH_impl(rng<RNGTYPE>& stream, const Matrix<>& Y,
- const Matrix<>& X, const Matrix<>& b0,
- const Matrix<>& B0, const Matrix<>& V,
- Matrix<>& beta, const Matrix<>& tune,
- unsigned int burnin, unsigned int mcmc,
- unsigned int thin, unsigned int verbose,
- Matrix<>& storemat)
-{
- // define constants
- const unsigned int tot_iter = burnin + mcmc; // total iterations
- const unsigned int nstore = mcmc / thin; // # of draws to store
- const unsigned int k = X.cols();
-
- // Initialize storage matrix
- storemat = Matrix<>(nstore, k, false);
+// natural log of the multivariate-t density (up to a constant of
+// proportionality)
+// theta: eval point
+// mu: mode
+// C: Cholesky factor of inverse scale matrix
+// df: degrees of freedom
+static inline double lnmulttdens(const Matrix<>& theta,
+ const Matrix<>& mu,
+ const Matrix<>& C,
+ const double& df){
+
+ const int d = theta.size();
+ //const Matrix<> z = t(theta - mu) * C;
+ // C is now C' if VC mat is C C'
+ const Matrix<> z = C * (theta - mu);
+ double zsumsq = 0;
+ for (int i=0; i<d; ++i){
+ zsumsq += std::pow(z[i], 2);
+ }
+
+ return ( (-(df + d)/2) * std::log(1 + (1/df) * zsumsq) );
+}
- // proposal parameters
- const Matrix<> propV = tune * invpd(B0 + invpd(V)) * tune;
- const Matrix<> propC = cholesky(propV) ;
-
- double logpost_cur = mnl_logpost(Y, X, beta, b0, B0);
-
- int count = 0;
- int accepts = 0;
- ///// MCMC SAMPLING OCCURS IN THIS FOR LOOP
- for (unsigned int iter = 0; iter < tot_iter; ++iter) {
-
- // sample beta
- const
- Matrix<> beta_can = gaxpy(propC, stream.rnorm(k,1,0,1), beta);
-
- const double logpost_can = mnl_logpost(Y, X, beta_can, b0, B0);
- const double ratio = std::exp(logpost_can - logpost_cur);
-
- if (stream() < ratio) {
- beta = beta_can;
- logpost_cur = logpost_can;
- ++accepts;
- }
-
- // store values in matrices
- if (iter >= burnin && ((iter % thin) == 0)) {
- for (unsigned int j = 0; j < k; j++)
- storemat(count, j) = beta[j];
- ++count;
- }
-
- // print output to stdout
- if (verbose > 0 && iter % verbose == 0) {
- Rprintf("\n\nMCMCmnl Metropolis iteration %i of %i \n",
- (iter+1), tot_iter);
- Rprintf("beta = \n");
- for (unsigned int j=0; j<k; ++j)
- Rprintf("%10.5f\n", beta[j]);
-
- Rprintf("Metropolis acceptance rate for beta = %3.5f\n\n",
- static_cast<double>(accepts) / static_cast<double>(iter+1));
- }
-
- R_CheckUserInterrupt(); // allow user interrupts
- } // end MCMC loop
+template <typename RNGTYPE>
+void MCMCmnlMH_impl(rng<RNGTYPE>& stream, const Matrix<>& Y,
+ const Matrix<>& X, const Matrix<>& b0,
+ const Matrix<>& B0, const Matrix<>& V,
+ Matrix<>& beta, const Matrix<>& beta_hat,
+ const Matrix<>& tune,
+ const unsigned int burnin, const unsigned int mcmc,
+ const unsigned int thin, const unsigned int verbose,
+ const unsigned int RW, const double tdf,
+ Matrix<>& storemat)
+{
+ // define constants
+ const unsigned int tot_iter = burnin + mcmc; // total iterations
+ const unsigned int nstore = mcmc / thin; // # of draws to store
+ const unsigned int k = X.cols();
+
+ // Initialize storage matrix
+ storemat = Matrix<>(nstore, k, false);
+
+ // proposal parameters
+ const Matrix<> propV = tune * V * tune;
+ const Matrix<> propC = cholesky(propV);
+ const Matrix<> propCinvT = t(cholesky(invpd(propV)));
+
+
+ double logpost_cur = mnl_logpost(Y, X, beta, b0, B0);
+ double logjump_cur = lnmulttdens(beta, beta_hat, propCinvT, tdf);
+
+ int count = 0;
+ int accepts = 0;
+
+ ///// MCMC SAMPLING OCCURS IN THIS FOR LOOP
+ for (unsigned int iter = 0; iter < tot_iter; ++iter) {
+
+ // sample beta
+ if (RW == 0){ // Independent Metropolis-Hastings
+
+ const double u = stream();
+
+ if (u < 0.75){
+ const Matrix<> beta_can = beta_hat + stream.rmvt(propV, tdf);
+ const double logpost_can = mnl_logpost(Y, X, beta_can, b0, B0);
+ const double logjump_can = lnmulttdens(beta_can, beta_hat,
+ propCinvT, tdf);
+
+ const double ratio = std::exp( logpost_can - logjump_can -
+ logpost_cur + logjump_cur );
+
+ if (stream() < ratio) {
+ beta = beta_can;
+ logpost_cur = logpost_can;
+ logjump_cur = logjump_can;
+ ++accepts;
+ }
+ }
+ else{
+ const Matrix<> beta_can = beta_hat + beta_hat - beta;
+ const double logpost_can = mnl_logpost(Y, X, beta_can, b0, B0);
+ const double logjump_can = lnmulttdens(beta_can, beta_hat,
+ propCinvT, tdf);
+
+ const double ratio = std::exp( logpost_can -
+ logpost_cur );
+
+ if (stream() < ratio) {
+ beta = beta_can;
+ logpost_cur = logpost_can;
+ logjump_cur = logjump_can;
+ ++accepts;
+ }
+
+ }
+
+
+ }
+ else{ // Random Walk Metropolis
+ const
+ Matrix<> beta_can = gaxpy(propC, stream.rnorm(k,1,0,1), beta);
+ const double logpost_can = mnl_logpost(Y, X, beta_can, b0, B0);
+ const double ratio = std::exp(logpost_can - logpost_cur);
+
+ if (stream() < ratio) {
+ beta = beta_can;
+ logpost_cur = logpost_can;
+ ++accepts;
+ }
+
+ }
+
+ // store values in matrices
+ if (iter >= burnin && ((iter % thin) == 0)) {
+ for (unsigned int j = 0; j < k; j++)
+ storemat(count, j) = beta[j];
+
+ ++count;
+ }
+
+ // print output to stdout
+ if (verbose > 0 && iter % verbose == 0) {
+ Rprintf("\n\nMCMCmnl Metropolis iteration %i of %i \n",
+ (iter+1), tot_iter);
+ Rprintf("beta = \n");
+
+ for (unsigned int j=0; j<k; ++j)
+ Rprintf("%10.5f\n", beta[j]);
+
+ Rprintf("Metropolis acceptance rate for beta = %3.5f\n\n",
+ static_cast<double>(accepts) / static_cast<double>(iter+1));
+ }
+
+ R_CheckUserInterrupt(); // allow user interrupts
+ } // end MCMC loop
}
extern "C" {
-
- // MCMC sampling for multinomial logit via Metropolis-Hastings
- void MCMCmnlMH(double *sampledata, const int *samplerow,
- const int *samplecol, const double *Ydata,
- const int *Yrow, const int *Ycol,
- const double *Xdata, const int *Xrow, const int *Xcol,
- const int *burnin, const int *mcmc, const int *thin,
- const double *tunedata, const int *tunerow,
- const int *tunecol, const int *uselecuyer,
- const int *seedarray, const int *lecuyerstream,
- const int *verbose, const double *betastartdata,
- const int *betastartrow, const int *betastartcol,
- const double *b0data, const int *b0row,
- const int *b0col, const double *B0data,
- const int *B0row, const int *B0col,
- const double *Vdata, const int *Vrow,
- const int *Vcol)
- {
-
- // pull together Matrix objects
- // REMEMBER TO ACCESS PASSED ints AND doubles PROPERLY
- const Matrix<> Y(*Yrow, *Ycol, Ydata);
- const Matrix<> X(*Xrow, *Xcol, Xdata);
- const Matrix<> tune(*tunerow, *tunecol, tunedata);
- Matrix<> beta(*betastartrow, *betastartcol, betastartdata);
- const Matrix<> b0(*b0row, *b0col, b0data);
- const Matrix<> B0(*B0row, *B0col, B0data);
- const Matrix<> V(*Vrow, *Vcol, Vdata);
-
- // storage matrix or matrices
- Matrix<> storemat;
- MCMCPACK_PASSRNG2MODEL(MCMCmnlMH_impl, Y, X, b0, B0, V, beta,
- tune, *burnin, *mcmc, *thin, *verbose, storemat);
-
- // load draws into sample array
- for(unsigned int i = 0; i < storemat.size(); ++i)
- sampledata[i] = storemat(i);
-
- } // end MCMCmnlMH
+
+ // MCMC sampling for multinomial logit via Metropolis-Hastings
+ void MCMCmnlMH(double *sampledata, const int *samplerow,
+ const int *samplecol, const double *Ydata,
+ const int *Yrow, const int *Ycol,
+ const double *Xdata, const int *Xrow, const int *Xcol,
+ const int *burnin, const int *mcmc, const int *thin,
+ const double *tunedata, const int *tunerow,
+ const int *tunecol, const int *uselecuyer,
+ const int *seedarray, const int *lecuyerstream,
+ const int *verbose, const double *betastartdata,
+ const int *betastartrow, const int *betastartcol,
+ const double *betamodedata,
+ const int *betamoderow, const int *betamodecol,
+ const double *b0data, const int *b0row,
+ const int *b0col, const double *B0data,
+ const int *B0row, const int *B0col,
+ const double *Vdata, const int *Vrow,
+ const int *Vcol, const int* RW, const double* tdf)
+ {
+
+ // pull together Matrix objects
+ // REMEMBER TO ACCESS PASSED ints AND doubles PROPERLY
+ const Matrix<> Y(*Yrow, *Ycol, Ydata);
+ const Matrix<> X(*Xrow, *Xcol, Xdata);
+ const Matrix<> tune(*tunerow, *tunecol, tunedata);
+ Matrix<> beta(*betastartrow, *betastartcol, betastartdata);
+ const Matrix<> betamode(*betamoderow, *betamodecol, betamodedata);
+ const Matrix<> b0(*b0row, *b0col, b0data);
+ const Matrix<> B0(*B0row, *B0col, B0data);
+ const Matrix<> V(*Vrow, *Vcol, Vdata);
+
+ // storage matrix or matrices
+ Matrix<> storemat;
+ MCMCPACK_PASSRNG2MODEL(MCMCmnlMH_impl, Y, X, b0, B0, V, beta, betamode,
+ tune, *burnin, *mcmc, *thin, *verbose, *RW,
+ *tdf, storemat);
+
+ // load draws into sample array
+ for(unsigned int i = 0; i < storemat.size(); ++i)
+ sampledata[i] = storemat(i);
+
+ } // end MCMCmnlMH
} // end extern "C"
#endif
diff --git a/src/MCMCoprobit.cc b/src/MCMCoprobit.cc
index 7327daf..6b10749 100644
--- a/src/MCMCoprobit.cc
+++ b/src/MCMCoprobit.cc
@@ -20,7 +20,8 @@
// 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
+
#ifndef MCMCOPROBIT_CC
#define MCMCOPROBIT_CC
@@ -35,6 +36,7 @@
#include "rng.h"
#include "mersenne.h"
#include "lecuyer.h"
+#include "optimize.h"
#include <R.h> // needed to use Rprintf()
#include <R_ext/Utils.h> // needed to allow user interrupts
@@ -42,15 +44,160 @@
using namespace std;
using namespace scythe;
+static inline double lnmulttdens(const Matrix<>& theta,
+ const Matrix<>& mu,
+ const Matrix<>& C,
+ const double df){
+
+ const int d = theta.size();
+ //const Matrix<> z = t(theta - mu) * C;
+ // C is now C' if VC mat is C C'
+ const Matrix<> z = C * (theta - mu);
+ double zsumsq = 0;
+ for (int i=0; i<d; ++i){
+ zsumsq += std::pow(z[i], 2);
+ }
+
+ return ( (-(df + d)/2) * std::log(1 + (1/df) * zsumsq) );
+}
+
+// function that transforms alpha to gamma
+Matrix<> gamma2alpha(const Matrix<>& gamma){
+ const int m = gamma.rows() - 2;
+ Matrix<> alpha(m, 1);
+ alpha[0] = std::log(gamma[1]);
+ for (int j=1; j< m ; ++j){
+ alpha[j] = std::log(gamma[j+1] - gamma[j]);
+ }
+ return alpha;
+}
+
+// function that transforms alpha to gamma
+Matrix<> alpha2gamma(const Matrix<>& alpha){
+ const int m = alpha.rows();
+ Matrix<> gamma(m+2, 1);
+ gamma[0] = -300;
+ gamma[m+1] = 300;
+ for (int j=1; j<m+1 ; ++j){
+ double gamma_sum = 0.0;
+ for(int i=0;i<j; ++i){
+ gamma_sum +=exp(alpha[i]);
+ }
+ gamma[j] = gamma_sum;
+ }
+ return gamma;
+}
+
+double lndmvn_jhp(const Matrix<double> &x, const Matrix<double> &mu,
+ const Matrix<double> &Sigma)
+{
+ int k = Sigma.cols();
+ double first = (-k/2.0) * ::log(2 * M_PI) -0.5 * ::log(det(Sigma));
+ Matrix< > second = ::t(x-mu)*invpd(Sigma)*(x-mu);
+ return (first - 0.5*second[0]);
+}
+
+
+
+// orpobit_logpost
+static inline double oprobit_logpost(const Matrix<double>& nY, const Matrix<double>& X,
+ const Matrix<double>& alpha,
+ const Matrix<double>& alpha_prior_mean,
+ const Matrix<double>& alpha_prior_var,
+ const Matrix<double>& beta){
+
+ // likelihood
+ double loglike = 0.0;
+ const int n = nY.rows();
+ const int ncat = alpha.rows() + 1;
+
+ // the linear predictor
+ Matrix<> mu = X*beta;
+ Matrix<> gamma = alpha2gamma(alpha);
+
+ // compute prob
+ Matrix<> cat_prob(n, ncat-1);
+ //cat_prob: CATegory specific PROBability
+ //the first col of cat_prob = pnorm(gamma_1 - mu)
+ //thus, the col number is ncat - 1
+ Matrix<> prob(n, ncat);
+ for (int j=0; j< ncat-1; ++j){
+ for (int i=0; i<n ; ++i){
+ cat_prob(i, j) = pnorm(gamma[j+1] - mu[i], 0.0, 1.0);
+ }
+ }
+ prob(_, ncat-1) = 1 - cat_prob(_, ncat-2); // top category
+ prob(_, 0) = cat_prob(_, 0); // bottom category
+ for (int j=1; j<(ncat-1); ++j){ // middle categories are actually differences of cumulatives
+ prob(_, j) = cat_prob(_,j) - cat_prob(_, j-1);
+ }
+
+ // the loglikelihood
+
+ for (int i=0; i<n; ++i){
+ int ind = (int) nY[i];
+ loglike += std::log(prob(i,ind-1));
+ }
+
+ // prior
+ double logprior = 0.0;
+ logprior = lndmvn_jhp(alpha, alpha_prior_mean, alpha_prior_var);
+
+ return (loglike + logprior);
+
+}// end of statis oprobit_logpost
+
+// The oprobitModel class stores additional data to be passed to BFGS
+class oprobitModel {
+public:
+ double operator() (const Matrix<double> alpha){
+ const int n = y_.rows();
+ const int ncat = alpha.rows() + 1;
+
+ // the linear predictor
+ Matrix<> mu = X_ * beta_;
+ Matrix<> gamma = alpha2gamma(alpha);
+
+ Matrix<> cat_prob(n, ncat-1);
+ //cat_prob: CATegory specific PROBability
+ //the first col of cat_prob = pnorm(gamma_1 - mu)
+ //thus, the col number is ncat - 1
+ Matrix<> prob(n, ncat);
+ for (int j=0; j< ncat-1; ++j){
+ for (int i=0; i<n ; ++i){
+ cat_prob(i, j) = pnorm(gamma[j+1] - mu[i], 0.0, 1.0);
+ }
+ }
+ prob(_, ncat-1) = 1 - cat_prob(_, ncat-2); // top category
+ prob(_, 0) = cat_prob(_, 0); // bottom category
+ for (int j=1; j<(ncat-1); ++j){ // middle categories are actually differences of cumulatives
+ prob(_, j) = cat_prob(_,j) - cat_prob(_, j-1);
+ }
+
+ // the loglikelihood
+ double loglike = 0.0;
+ for (int i=0; i<n; ++i){
+ int ind = y_[i];
+ loglike += std::log(prob(i,ind-1));
+ }
+ return -1 * loglike;
+ }
+ Matrix<double> y_;
+ Matrix<double> X_;
+ Matrix<double> beta_;
+};
+
/* MCMCoprobit implementation. Takes Matrix<> reference which it
* fills with the posterior.
*/
template <typename RNGTYPE>
void MCMCoprobit_impl (rng<RNGTYPE>& stream, const int * Y,
- const Matrix<>& X, Matrix<>& beta, Matrix<>& gamma,
+ const Matrix<>& nY, const Matrix<>& X, Matrix<>& beta, Matrix<>& gamma,
const Matrix<>& b0, const Matrix<>& B0,
+ const Matrix<>& alpha_prior_mean, const Matrix<>& alpha_prior_var,
const unsigned int burnin, const unsigned int mcmc, const unsigned int thin,
- const unsigned int verbose, const double tune, Matrix<>& result) {
+ const unsigned int verbose, const Matrix<>& tune, const double tdf,
+ const unsigned int cowles, Matrix<>& result) {
// define constants and from cross-product matrices
const unsigned int tot_iter = burnin + mcmc; // total number of mcmc iterations
@@ -62,114 +209,149 @@ void MCMCoprobit_impl (rng<RNGTYPE>& stream, const int * Y,
// storage matrix or matrices
Matrix<> storemat(nstore, k+ncat+1);
-
- // initialize Z
- Matrix<> gamma_p = gamma;
+
+ // initialize Z
Matrix<> Z(N,1);
Matrix<> Xbeta = X * beta;
+ // Gibbs loop
+ int count = 0;
+ int accepts = 0;
+ Matrix<> gamma_p = gamma;
+ Matrix<> gamma_new = gamma + 1;
+ Matrix<> alpha = gamma2alpha(gamma_new);
+ Matrix<> alpha_hat = alpha;
- // Gibbs loop
- int count = 0;
- int accepts = 0;
- for (unsigned int iter = 0; iter < tot_iter; ++iter) {
+ // initialize current value stuff
+ Matrix<> propV = tune * alpha_prior_var * tune;
+ Matrix<> propCinvT = t(cholesky(invpd(propV)));
+ double logpost_cur = oprobit_logpost(nY, X, alpha, alpha_prior_mean, alpha_prior_var, beta);
+ double logjump_cur = lnmulttdens(alpha_prior_mean, alpha_hat, propCinvT, tdf);
+
+ double tune_double = tune[0];
+ for (unsigned int iter = 0; iter < tot_iter; ++iter) {
+
+ //////////////////
+ if (cowles == 1){
+ //////////////////
+ // Cowles method [gamma | Z, beta]
+ for (int i=2; i<(ncat); ++i){
+ if (i==(ncat-1)){
+ gamma_p[i] = stream.rtbnorm_combo(gamma[i], ::pow(tune_double, 2.0),
+ gamma_p[i-1]);
+ }
+ else {
+ gamma_p[i] = stream.rtnorm_combo(gamma[i], ::pow(tune_double, 2.0),
+ gamma_p[i-1],
+ gamma[i+1]);
+ }
+ }
+ double loglikerat = 0.0;
+ double loggendenrat = 0.0;
+
+ // loop over observations and construct the acceptance ratio
+ for (unsigned int i=0; i<N; ++i){
+ if (Y[i] == ncat){
+ loglikerat = loglikerat
+ + log(1.0 - pnorm(gamma_p[Y[i]-1] - Xbeta[i], 0.0, 1.0) )
+ - log(1.0 - pnorm(gamma[Y[i]-1] - Xbeta[i], 0.0, 1.0) );
+ }
+ else if (Y[i] == 1){
+ loglikerat = loglikerat + log(pnorm(gamma_p[Y[i]] - Xbeta[i], 0.0, 1.0) )
+ - log(pnorm(gamma[Y[i]] - Xbeta[i], 0.0, 1.0) );
+ }
+ else{
+ loglikerat = loglikerat
+ + log(pnorm(gamma_p[Y[i]] - Xbeta[i], 0.0, 1.0) -
+ pnorm(gamma_p[Y[i]-1] - Xbeta[i], 0.0, 1.0) )
+ - log(pnorm(gamma[Y[i]] - Xbeta[i], 0.0, 1.0) -
+ pnorm(gamma[Y[i]-1] - Xbeta[i], 0.0, 1.0) );
+ }
+ }
+ double logacceptrat = loglikerat + loggendenrat;
+ if (stream.runif() <= exp(logacceptrat)){
+ gamma = gamma_p;
+ ++accepts;
+ }
+
+ }// end of Cowles
- // [gamma | Z, beta]
- for (int i=2; i<(ncat); ++i){
- if (i==(ncat-1)){
- gamma_p[i] = stream.rtbnorm_combo(gamma[i], ::pow(tune, 2.0),
- gamma_p[i-1]);
- }
+ //////////////////
else {
- gamma_p[i] = stream.rtnorm_combo(gamma[i], ::pow(tune, 2.0),
- gamma_p[i-1],
- gamma[i+1]);
- }
- }
- double loglikerat = 0.0;
- double loggendenrat = 0.0;
-
- // loop over observations and construct the acceptance ratio
- for (unsigned int i=0; i<N; ++i){
- if (Y[i] == ncat){
- loglikerat = loglikerat
- + log(1.0 -
- pnorm(gamma_p[Y[i]-1] - Xbeta[i], 0.0, 1.0) )
- - log(1.0 -
- pnorm(gamma[Y[i]-1] - Xbeta[i], 0.0, 1.0) );
- }
- else if (Y[i] == 1){
- loglikerat = loglikerat
- + log(pnorm(gamma_p[Y[i]] - Xbeta[i], 0.0, 1.0) )
- - log(pnorm(gamma[Y[i]] - Xbeta[i], 0.0, 1.0) );
- }
- else{
- loglikerat = loglikerat
- + log(pnorm(gamma_p[Y[i]] - Xbeta[i], 0.0, 1.0) -
- pnorm(gamma_p[Y[i]-1] - Xbeta[i], 0.0, 1.0) )
- - log(pnorm(gamma[Y[i]] - Xbeta[i], 0.0, 1.0) -
- pnorm(gamma[Y[i]-1] - Xbeta[i], 0.0, 1.0) );
- }
- }
- // for (int j=2; j<ncat; ++j){
- // loggendenrat = loggendenrat
- // + log(pnorm(gamma[j+1], gamma[j], tune) -
- // pnorm(gamma_p[j-1], gamma[j], tune) )
- // - log(pnorm(gamma_p[j+1], gamma_p[j], tune) -
- // pnorm(gamma[j-1], gamma_p[j], tune) );
- // }
- double logacceptrat = loglikerat + loggendenrat;
- if (stream.runif() <= exp(logacceptrat)){
- gamma = gamma_p;
- ++accepts;
- }
-
+ //////////////////
+ // Albert and Chib (2001) method
+ // Step 1: [gamma | Z, beta]
+ oprobitModel oprobit_model;
+ oprobit_model.y_ = nY;
+ oprobit_model.X_ = X;
+ oprobit_model.beta_ = beta;
+
+ // tailored proposal MH
+ Matrix<> alpha_hat = BFGS(oprobit_model, alpha, stream, 100, 1e-5, false);
+ Matrix<> alpha_V = invpd(hesscdif(oprobit_model, alpha_hat)); //note that oprobit_model contains the multiplication by -1
+ Matrix<> propV = tune * alpha_V * tune;
+ Matrix<> propCinvT = ::t(cholesky(invpd(propV)));
+
+ // Draw alpha_can from multivariate t
+ Matrix<> alpha_can = alpha_hat + stream.rmvt(propV, tdf);
+
+ // compute components of transition kernel
+ double logpost_can = oprobit_logpost(nY, X, alpha_can, alpha_prior_mean, alpha_prior_var, beta);
+ double logjump_can = lnmulttdens(alpha_can, alpha_hat, propCinvT, tdf);
+ double logpost_cur = oprobit_logpost(nY, X, alpha, alpha_prior_mean, alpha_prior_var, beta);
+ double logjump_cur = lnmulttdens(alpha, alpha_hat, propCinvT, tdf);
+ double ratio = exp(logpost_can - logjump_can - logpost_cur + logjump_cur);
+ const double u = stream();
+ if (u < ratio) {
+ alpha = alpha_can;
+ gamma = alpha2gamma(alpha);
+ logpost_cur = logpost_can;
+ logjump_cur = logjump_can;
+ ++accepts;
+ }
+ }// end of AC method
+
+ // Step 2: [Z| gamma, beta, y]
+ Xbeta = X * beta;
+ // Matrix<> Z_noconst(N,1);
+ for (unsigned int i=0; i<N; ++i){
+ Z[i] = stream.rtnorm_combo(Xbeta[i], 1.0, gamma[Y[i]-1], gamma[Y[i]]);
+ }
+
+ // Step 3: [beta|Z, gamma]
+ const Matrix<> XpZ = t(X) * Z;
+ beta = NormNormregress_beta_draw(XpX, XpZ, b0, B0, 1.0, stream);
+
+ // store values in matrices
+ if (iter >= burnin && ((iter % thin)==0)){
+ for (unsigned int j=0; j<k; ++j)
+ storemat(count, j) = beta[j];
+ for (int j=0; j<(ncat+1); ++j)
+ storemat(count, j+k) = gamma[j];
+ ++count;
+ }
+
+ // print output to stdout
+ if(verbose > 0 && iter % verbose == 0){
+ Rprintf("\n\nMCMCoprobit iteration %i of %i \n", (iter+1), tot_iter);
+ Rprintf("beta = \n");
+ Rprintf("%10.5f\n", beta[0]-gamma[1]);
+ for (unsigned int j=1; j<k; ++j)
+ Rprintf("%10.5f\n", beta[j]);
+ Rprintf("Metropolis acceptance rate for gamma = %3.5f\n\n",
+ static_cast<double>(accepts)/static_cast<double>(iter+1));
+ }
+
+ R_CheckUserInterrupt(); // allow user interrupts
+ }// end of MCMC
+
+ Rprintf("\n\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n");
+ Rprintf("The Metropolis acceptance rate for beta was %3.5f",
+ static_cast<double>(accepts) / static_cast<double>(tot_iter));
+ Rprintf("\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n");
- // [Z| gamma, beta, y]
- //Matrix<double> Z_mean = X * beta;
- for (unsigned int i=0; i<N; ++i){
- Z[i] = stream.rtnorm_combo(Xbeta[i], 1.0, gamma[Y[i]-1], gamma[Y[i]]);
- }
-
- // [beta|Z, gamma]
- const Matrix<> XpZ = t(X) * Z;
- beta = NormNormregress_beta_draw(XpX, XpZ, b0, B0, 1.0, stream);
- Xbeta = X * beta;
-
-
- // store values in matrices
- if (iter >= burnin && ((iter % thin)==0)){
- for (unsigned int j=0; j<k; ++j)
- storemat(count, j) = beta[j];
- for (int j=0; j<(ncat+1); ++j)
- storemat(count, j+k) = gamma[j];
- ++count;
- }
-
-
-
-
- // print output to stdout
- if(verbose > 0 && iter % verbose == 0){
- Rprintf("\n\nMCMCoprobit iteration %i of %i \n", (iter+1), tot_iter);
- Rprintf("beta = \n");
- for (unsigned int j=0; j<k; ++j)
- Rprintf("%10.5f\n", beta[j]);
- Rprintf("Metropolis acceptance rate for gamma = %3.5f\n\n",
- static_cast<double>(accepts) /
- static_cast<double>(iter+1));
- }
-
- R_CheckUserInterrupt(); // allow user interrupts
- }
result = storemat;
-
-
- Rprintf("\n\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n");
- Rprintf("The Metropolis acceptance rate for beta was %3.5f",
- static_cast<double>(accepts) / static_cast<double>(tot_iter));
- Rprintf("\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n");
}
@@ -178,30 +360,39 @@ extern "C"{
void MCMCoprobit(double *sampledata, const int *samplerow,
const int *samplecol, const int *Y,
+ const double *nYdata, const int *nYrow, const int *nYcol,
const double *Xdata,
const int *Xrow, const int *Xcol, const int *burnin,
- const int *mcmc, const int *thin, const double* tune,
+ const int *mcmc, const int *thin, const double *tunedata,
+ const int *tunerow, const int *tunecol,
+ const double* tdf,
const int *uselecuyer, const int *seedarray,
const int *lecuyerstream, const int *verbose,
const double *betadata, const int *betarow,
const int *betacol, const double* gammadata,
const int* gammarow, const int* gammacol,
const double *b0data, const int *b0row, const int *b0col,
- const double *B0data, const int *B0row, const int *B0col) {
+ const double *B0data, const int *B0row, const int *B0col,
+ const double *a0data, const int *a0row, const int *a0col,
+ const double *A0data, const int *A0row, const int *A0col,
+ const int *cowles) {
// pull together Matrix objects
+ const Matrix <> nY(*nYrow, *nYcol, nYdata);
const Matrix <> X(*Xrow, *Xcol, Xdata);
- Matrix <> beta(*betarow, *betacol,
- betadata);
- Matrix <> gamma(*gammarow, *gammacol,
- gammadata);
+ Matrix <> beta(*betarow, *betacol, betadata);
+ Matrix <> gamma(*gammarow, *gammacol, gammadata);
const Matrix <> b0(*b0row, *b0col, b0data);
const Matrix <> B0(*B0row, *B0col, B0data);
-
+ const Matrix <> alpha_prior_mean(*a0row, *a0col, a0data);
+ const Matrix <> alpha_prior_var(*A0row, *A0col, A0data);
+ const Matrix<> tune(*tunerow, *tunecol, tunedata);
+
Matrix<> storagematrix;
- MCMCPACK_PASSRNG2MODEL(MCMCoprobit_impl, Y, X, beta, gamma, b0, B0,
- *burnin, *mcmc, *thin, *verbose, *tune,
- storagematrix);
+ MCMCPACK_PASSRNG2MODEL(MCMCoprobit_impl, Y, nY, X, beta, gamma, b0, B0,
+ alpha_prior_mean, alpha_prior_var,
+ *burnin, *mcmc, *thin, *verbose, tune, *tdf,
+ *cowles, storagematrix);
const unsigned int size = *samplerow * *samplecol;
for (unsigned int i = 0; i < size; ++i)
diff --git a/src/MCMCpoissonChangepoint.cc b/src/MCMCpoissonChangepoint.cc
new file mode 100644
index 0000000..da68951
--- /dev/null
+++ b/src/MCMCpoissonChangepoint.cc
@@ -0,0 +1,382 @@
+// MCMCpoissonChange.cc is C++ code to estimate a Poisson changepoint model
+// with a gamma prior
+//
+// Jong Hee Park
+// Dept. of Political Science
+// University of Chicago
+// jhp at uchicago.edu
+//
+// This software is distributed under the terms of the GNU GENERAL
+// PUBLIC LICENSE Version 2, June 1991. See the package LICENSE
+// file for more information.
+//
+// Copyright (C) 2004 Andrew D. Martin and Kevin M. Quinn
+
+#ifndef MCMCPOISSONCHANGEPOINT_CC
+#define MCMCPOISSONCHANGEPOINT_CC
+
+#include "MCMCrng.h"
+#include "MCMCfcds.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
+
+using namespace std;
+using namespace scythe;
+
+
+template <typename RNGTYPE>
+Matrix<double> poisson_state_sampler(rng<RNGTYPE>& stream, const int& m,
+ const Matrix<double>& Y,
+ const Matrix<double>& lambda,
+ const Matrix<double>& P){
+ const int ns = m + 1;
+ const int n = Y.rows();
+ Matrix<> F(n, ns);
+ Matrix<> pr1(ns, 1);
+ pr1[0] = 1;
+ Matrix<> py(ns, 1);
+ Matrix<> pstyt1(ns, 1);
+
+ //
+ // Forward sampling: update F matrix
+ //
+ for (int t=0; t<n ; ++t){
+ int yt = (int) Y[t];
+ for (int j=0; j<ns ; ++j){
+ py[j] = dpois(yt, lambda[j]);
+ }
+ if (t==0) pstyt1 = pr1;
+ else {
+ pstyt1 = ::t(F(t-1,_)*P);
+ }
+ Matrix<> unnorm_pstyt = pstyt1%py;
+ Matrix<> pstyt = unnorm_pstyt/sum(unnorm_pstyt); // pstyt = Pr(st|Yt)
+ for (int j=0; j<ns ; ++j){
+ F(t,j) = pstyt(j);
+ }
+ }// end of F matrix filtering
+
+ //
+ // Backward recursions
+ //
+ Matrix<int> s(n, 1);
+ Matrix<> ps(n, ns);
+ ps(n-1,_) = F(n-1,_);
+ s(n-1) = ns;
+
+ Matrix<> pstyn(ns, 1);
+ double pone = 0.0;
+ int t = n-2;
+ while (t >= 0){
+ int st = s(t+1);
+ Matrix<> Pst_1 = ::t(P(_,st-1));
+ Matrix<> unnorm_pstyn = F(t,_)%Pst_1;
+ pstyn = unnorm_pstyn/sum(unnorm_pstyn);
+ if (st==1) s(t) = 1;
+ else{
+ pone = pstyn(st-2);
+ if(stream.runif () < pone) s(t) = st-1;
+ else s(t) = st;
+ }
+ ps(t,_) = pstyn;
+ --t;
+ }// end of while loop
+
+ Matrix<> Sout(n, ns+1);
+ Sout(_, 0) = s(_,0);
+ for (int j = 0; j<ns; ++j){
+ Sout(_,j+1) = ps(_, j);
+ }
+
+ return Sout;
+
+ } // end of state sampler
+
+////////////////////////////////////////////
+// MCMCpoissonChangepoint implementation.
+////////////////////////////////////////////
+template <typename RNGTYPE>
+void MCMCpoissonChangepoint_impl(rng<RNGTYPE>& stream, const Matrix<>& Y,
+ Matrix<>& lambda, Matrix<>& P, const Matrix<>& A0,
+ double m, double c0, double d0,
+ unsigned int burnin, unsigned int mcmc, unsigned int thin,
+ unsigned int verbose, bool chib,
+ Matrix<>& lambda_store, Matrix<>& P_store,
+ Matrix<>& ps_store, Matrix<>& s_store,
+ double& logmarglike)
+{
+ // define constants and form cross-product matrices
+ const unsigned int tot_iter = burnin + mcmc; //total iterations
+ const unsigned int nstore = mcmc / thin; // number of draws to store
+ const int n = Y.rows();
+ const int ns = m + 1; // number of states
+
+ //MCMC loop
+ unsigned int count = 0;
+ for (int iter = 0; iter < tot_iter; ++iter){
+
+ //////////////////////
+ // 1. Sample s
+ //////////////////////
+ Matrix<> Sout = poisson_state_sampler(stream, m, Y, lambda, P);
+ Matrix<> s = Sout(_, 0);
+ Matrix<> ps(n, ns);
+ for (int j = 0; j<ns; ++j){
+ ps(_,j) = Sout(_,j+1);
+ }
+
+ //////////////////////
+ // 2. Sample lambda
+ //////////////////////
+ Matrix<> addY(ns, 1);
+ Matrix<> addN(ns, 1);
+
+ for (int j = 0; j<ns; ++j){
+ for (int i = 0; i<n; ++i){
+ if (s[i] == (j+1)) { // since j starts from 0
+ addN[j] = addN[j] + 1;
+ addY[j] = addY[j] + Y[i];
+ }// end of if
+ }
+ double c1 = addY[j] + c0;
+ double d1 = addN[j] + 1/ d0;
+ lambda[j] = stream.rgamma(c1, d1);
+ }
+
+ //////////////////////
+ // 3. Sample P
+ //////////////////////
+ double shape1 = 0;
+ double shape2 = 0;
+ P(ns-1, ns-1) = 1; //no jump at the last state
+
+ for (int j =0; j< (ns-1); ++j){
+ shape1 = A0(j,j) + addN[j] - 1;
+ shape2 = A0(j,j+1) + 1; // SS(j,j+1);
+ P(j,j) = stream.rbeta(shape1, shape2);
+ P(j,j+1) = 1 - P(j,j);
+ }
+
+ // load draws into sample array
+ if (iter >= burnin && ((iter % thin)==0)){
+ for (int i=0; i<ns; ++i)
+ lambda_store(count,i) = lambda[i];
+ for (int j=0; j<ns*ns; ++j)
+ P_store(count,j)= P[j];
+ s_store(count,_) = s;
+ for (int l=0; l<n ; ++l)
+ ps_store(l,_) = ps_store(l,_) + ps(l,_); // add two matrices
+
+ ++count; // count when (iter % *thin)==0
+
+ } // end of if(iter...)
+
+
+ // print output to stdout
+ if(verbose > 0 && iter % verbose == 0){
+ Rprintf("\n\nMCMCpoissonChangepoint iteration %i of %i \n", (iter+1), tot_iter);
+ Rprintf("lambda = \n");
+ for (int j=0; j<ns; ++j)
+ Rprintf("%10.5f\n", lambda[j]);
+ }
+
+ R_CheckUserInterrupt(); // allow user interrupts
+
+ }// end MCMC loop
+
+if(chib ==1){
+
+ Matrix<> lambda_st = meanc(lambda_store);
+ Matrix<> P_vec_st = meanc(P_store);
+ const Matrix<> P_st(ns, ns);
+ for (int j = 0; j< ns*ns; ++j){
+ P_st[j] = P_vec_st[j];
+ }
+
+ //////////////////////
+ // 1. pdf.lambda
+ //////////////////////
+ Matrix<> density_lambda(nstore, ns);
+ for (int iter = 0; iter<nstore; ++iter){
+
+ Matrix<> addY(ns, 1);
+ Matrix<> addN(ns, 1);
+
+ for (int j = 0; j<ns; ++j){
+ for (int i = 0; i<n; ++i){
+ if (s_store(iter, i) == (j+1)) {
+ addN[j] = addN[j] + 1;
+ addY[j] = addY[j] + Y[i];
+ }// end of if
+ } // end of for i
+ double c1 = addY[j] + c0;
+ double d1 = addN[j] + d0;
+ density_lambda(iter, j) = dgamma(lambda_st[j], c1, 1/d1);
+ } // end of for j
+
+ }// end of pdf.lambda MCMC run
+ double pdf_lambda = log(prod(meanc(density_lambda)));
+
+ //////////////////////
+ // 2. pdf.P
+ //////////////////////
+ Matrix<> density_P(nstore, ns);
+ for (int iter = 0; iter < nstore; ++iter){
+ Matrix<> Sout = poisson_state_sampler(stream, m, Y, lambda_st, P);
+ Matrix <> s = Sout(_, 0);
+ Matrix <> ps(n, ns);
+ for (int j = 0; j<ns; ++j){
+ ps(_,j) = Sout(_,j+1);
+ }
+
+ double shape1 = 0;
+ double shape2 = 0;
+ P(ns-1, ns-1) = 1;
+ Matrix <> P_addN(ns, 1);
+ for (int j = 0; j<ns; ++j){
+ for (int i = 0; i<n; ++i){
+ if (s[i] == (j+1)) { // since j starts from 0
+ P_addN[j] = P_addN[j] + 1;
+ }// end of if
+ } // end of for i
+ } // end of for j
+
+ for (int j =0; j< (ns-1); ++j){
+ shape1 = A0(j,j) + P_addN[j] - 1;
+ shape2 = A0(j,j+1) + 1; //
+ P(j,j) = stream.rbeta(shape1, shape2);
+ P(j,j+1) = 1 - P(j,j);
+ density_P(iter, j) = dbeta(P_st(j,j), shape1, shape2);
+ } // end of for j
+ density_P(iter, ns-1) = 1; //
+
+ }// end of pdf.P MCMC run
+ double pdf_P = log(prod(meanc(density_P)));
+
+ //////////////////////
+ // likelihood
+ //////////////////////
+ Matrix<> F(n, ns);
+ Matrix<> like(n, 1);
+ Matrix<> pr1(ns, 1);
+ pr1[0] = 1;
+ Matrix<> py(ns, 1);
+ Matrix<> pstyt1(ns, 1);
+
+ for (int t=0; t<n ; ++t){
+ int yt = (int) Y[t];
+ for (int j=0; j<ns ; ++j){
+ py[j] = dpois(yt, lambda_st[j]);
+ }
+ if (t==0) pstyt1 = pr1;
+ else {
+ pstyt1 = ::t(F(t-1,_)*P_st); // make it an ns by 1 matrix
+ }
+ Matrix<double> unnorm_pstyt = pstyt1%py;
+ Matrix<double> pstyt = unnorm_pstyt/sum(unnorm_pstyt); // pstyt = Pr(st|Yt)
+ for (int j=0; j<ns ; ++j){
+ F(t,j) = pstyt(j);
+ }
+ like[t] = sum(unnorm_pstyt);
+ }// end of likelihood computation
+
+ double loglike = sum(log(like));
+
+ Rprintf("loglike is \n");
+ Rprintf("%10.5f\n", loglike);
+
+ //////////////////////
+ // log prior ordinate
+ //////////////////////
+ Matrix<> density_lambda_prior(ns, 1);
+ Matrix<> density_P_prior(ns, 1);
+ density_P[ns-1] = 1; //
+
+ for (int j=0; j<ns ; ++j){
+ density_lambda_prior[j] = log(dgamma(lambda_st[j], c0, d0));
+ }
+
+ for (int j =0; j< (ns-1); ++j){
+ density_P_prior[j] = log(dbeta(P_st(j,j), A0(j,j), A0(j,j+1)));
+ }
+
+ // compute marginal likelihood
+ double logprior = sum(density_lambda_prior) + sum(density_P_prior);
+ logmarglike = (loglike + logprior) - (pdf_lambda + pdf_P);
+ Rprintf("logmarglike is \n");
+ Rprintf("%10.5f\n", logmarglike);
+ Rprintf("-------------------------------------------\n");
+
+
+}// end of marginal likelihood computation
+
+}
+////////////////////////////////////////////
+// Start MCMCpoissonChangepoint function
+///////////////////////////////////////////
+extern "C"{
+ void MCMCpoissonChangepoint(double *lambdaout, double *Pout, double *psout, double *sout,
+ const double *Ydata, const int *Yrow, const int *Ycol, const int *m,
+ const int *burnin, const int *mcmc, const int *thin, const int *uselecuyer, const int *seedarray,
+ const int *lecuyerstream, const int *verbose,
+ const double *lambdastart, const double *Pstart,
+ const double *a, const double *b, const double *c0, const double *d0,
+ const double *A0data, double *logmarglikeholder,
+ const int *chib){
+
+ // pull together Matrix objects
+ const Matrix <double> Y(*Yrow, *Ycol, Ydata);
+ const unsigned int tot_iter = *burnin + *mcmc; //total iterations
+ const unsigned int nstore = *mcmc / *thin; // number of draws to store
+ const int n = Y.rows();
+ const int ns = *m + 1; // number of states
+
+ // generate starting values
+ Matrix <> lambda(ns, 1, lambdastart);
+ const Matrix <> A0(ns, ns, A0data);
+ Matrix <> P(ns, ns, Pstart);
+
+ double logmarglike;
+ // double loglike;
+
+ // storage matrices
+ Matrix<> lambda_store(nstore, ns);
+ Matrix<> P_store(nstore, ns*ns);
+ Matrix<> ps_store(n, ns);
+ Matrix<> s_store(nstore, n);
+
+ MCMCPACK_PASSRNG2MODEL(MCMCpoissonChangepoint_impl, Y,
+ lambda, P, A0, *m, *c0, *d0, *burnin, *mcmc, *thin, *verbose, *chib,
+ lambda_store, P_store, ps_store, s_store, logmarglike)
+
+ logmarglikeholder[0] = logmarglike;
+ // loglikeholder[0] = loglike;
+
+ // return output
+ for (int i = 0; i<(nstore*ns); ++i){
+ lambdaout[i] = lambda_store[i];
+ }
+ for (int i = 0; i<(nstore*ns*ns); ++i){
+ Pout[i] = P_store[i];
+ }
+ for (int i = 0; i<(n*ns); ++i){
+ psout[i] = ps_store[i];
+ }
+ for (int i = 0; i<(nstore*n); ++i){
+ sout[i] = s_store[i];
+ }
+
+ }// end of MCMCpoissonChange
+}//end extern "C"
+
+#endif
diff --git a/src/Makevars b/src/Makevars
index 8f3a8df..8e1d77d 100644
--- a/src/Makevars
+++ b/src/Makevars
@@ -1 +1,2 @@
-PKG_CXXFLAGS = -DSCYTHE_COMPILE_DIRECT -DSCYTHE_DEBUG=0 -DHAVE_TRUNC -DHAVE_TRUNC
+PKG_CPPFLAGS = -DSCYTHE_COMPILE_DIRECT -DSCYTHE_DEBUG=0 -DHAVE_TRUNC -DHAVE_TRUNC
+
diff --git a/src/Makevars.in b/src/Makevars.in
index afaee3c..ad21240 100644
--- a/src/Makevars.in
+++ b/src/Makevars.in
@@ -1 +1,2 @@
-PKG_CXXFLAGS = -DSCYTHE_COMPILE_DIRECT -DSCYTHE_DEBUG=0 -DHAVE_TRUNC @MV_HAVE_IEEEFP_H@ @MV_HAVE_TRUNC@
+PKG_CPPFLAGS = -DSCYTHE_COMPILE_DIRECT -DSCYTHE_DEBUG=0 -DHAVE_TRUNC @MV_HAVE_IEEEFP_H@ @MV_HAVE_TRUNC@
+
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/r-cran-mcmcpack.git
More information about the debian-science-commits
mailing list