[r-cran-mcmcpack] 52/90: Imported Upstream version 1.0-9
Andreas Tille
tille at debian.org
Fri Dec 16 09:07:46 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 a7517a1938eaaed18ad696e5f627faf93f75eba6
Author: Andreas Tille <tille at debian.org>
Date: Fri Dec 16 08:07:19 2016 +0100
Imported Upstream version 1.0-9
---
DESCRIPTION | 8 +-
HISTORY | 10 +-
NAMESPACE | 2 +
R/MCMCoprobitChange.R | 222 ++++++++++++++
R/MCMCpoissonChange.R | 3 +
R/MCMCprobit.R | 79 +++--
R/MCMCprobitChange.R | 123 ++++++++
man/MCMCmetrop1R.Rd | 3 +
man/MCMCmnl.Rd | 3 +
man/MCMCoprobit.Rd | 11 +-
man/MCMCoprobitChange.Rd | 183 ++++++++++++
man/MCMCpoissonChange.Rd | 22 +-
man/MCMCprobit.Rd | 31 +-
man/MCMCprobitChange.Rd | 167 +++++++++++
man/MCMCtobit.Rd | 5 +-
src/MCMCoprobitChange.cc | 747 +++++++++++++++++++++++++++++++++++++++++++++++
src/MCMCpoissonChange.cc | 11 +-
src/MCMCprobit.cc | 150 ++++++----
src/MCMCprobitChange.cc | 415 ++++++++++++++++++++++++++
src/MCMCprobitres.cc | 67 ++++-
20 files changed, 2144 insertions(+), 118 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index afdd562..5a07c7b 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
Package: MCMCpack
-Version: 1.0-8
-Date: 2010-10-17
+Version: 1.0-9
+Date: 2011-1-31
Title: Markov chain Monte Carlo (MCMC) Package
Author: Andrew D. Martin <admartin at wustl.edu>, Kevin M. Quinn
<kquinn at law.berkeley.edu>, Jong Hee Park <jhp at uchicago.edu>
@@ -19,6 +19,6 @@ Description: This package contains functions to perform Bayesian
License: GPL-2
SystemRequirements: gcc (>= 4.0)
URL: http://mcmcpack.wustl.edu
-Packaged: 2010-10-17 18:39:16 UTC; adm
+Packaged: 2011-01-31 17:28:11 UTC; adm
Repository: CRAN
-Date/Publication: 2010-10-19 16:29:53
+Date/Publication: 2011-01-31 18:50:39
diff --git a/HISTORY b/HISTORY
index d0b4df8..17a9085 100644
--- a/HISTORY
+++ b/HISTORY
@@ -2,6 +2,10 @@
// Changes and Bug Fixes
//
+1.0-8 to 1.0-9
+ * added MCMCprobitChange()
+ * added MCMCoprobitChange()
+ * modified MCMCprobit and MCMCprobitres for marginal likelihood estimation
1.0-7 to 1.0-8
* fixed some NAMESPACE issues [thanks to Shawn Treier]
@@ -31,12 +35,12 @@
* added Bayesian quantile regression [contributed by Craig Reed]
1.0-1 to 1.0-2
- * Added Poisson regression changepoint analysis MCMCpoissonChange() [JHP]
+ * added Poisson regression changepoint analysis MCMCpoissonChange() [JHP]
* Old MCMCpoissonChangepoint() is deprecated [JHP]
- * Added binary data changepoint model MCMCbinaryChange() [JHP]
+ * added binary data changepoint model MCMCbinaryChange() [JHP]
* plotState() function modified to support new models [JHP]
* plotChangepoint() function modified to support new models
- * Added a number of new helper functions in btsutil.R [JHP]
+ * added a number of new helper functions in btsutil.R [JHP]
0.9-6 to 1.0-1
* added one-dimensional dynamic IRT model
diff --git a/NAMESPACE b/NAMESPACE
index 50ae65c..811315c 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -31,10 +31,12 @@ export(
MCMCmixfactanal,
MCMCmnl,
MCMCoprobit,
+ MCMCoprobitChange,
MCMCordfactanal,
MCMCpoisson,
MCMCpoissonChange,
MCMCprobit,
+ MCMCprobitChange,
MCMCquantreg,
MCMCregress,
MCMCSVDreg,
diff --git a/R/MCMCoprobitChange.R b/R/MCMCoprobitChange.R
new file mode 100644
index 0000000..864e0e7
--- /dev/null
+++ b/R/MCMCoprobitChange.R
@@ -0,0 +1,222 @@
+#########################################################
+##
+## sample from the posterior distribution
+## of ordinal probit changepoint regression model
+## using a linear Gaussian approximation
+##
+## JHP 07/01/2007
+## JHP 03/03/2009
+## JHP 09/08/2010
+#########################################################
+
+"MCMCoprobitChange"<-
+ function(formula, data=parent.frame(), m = 1,
+ burnin = 1000, mcmc = 1000, thin = 1, tune = NA, verbose = 0,
+ seed = NA, beta.start = NA, gamma.start = NA, P.start = NA,
+ b0 = NULL, B0 = NULL, a = NULL, b = NULL,
+ marginal.likelihood = c("none", "Chib95"),
+ gamma.fixed=0, ...){
+
+ ## checks
+ check.offset(list(...))
+ check.mcmc.parameters(burnin, mcmc, thin)
+ cl <- match.call()
+ nstore <- mcmc/thin
+
+ ## seeds
+ seeds <- form.seeds(seed)
+ lecuyer <- seeds[[1]]
+ seed.array <- seeds[[2]]
+ lecuyer.stream <- seeds[[3]]
+
+ totiter <- mcmc+burnin
+ holder <- parse.formula(formula, data=data)
+ y <- holder[[1]]
+ X <- holder[[2]]
+ xnames <- holder[[3]]
+ K <- ncol(X)
+ Y <- factor(y, ordered = TRUE)
+ ncat <- nlevels(Y)
+ cat <- levels(Y)
+ ns <- m + 1
+ N <- nrow(X)
+ gk <- ncat + 1
+
+ if(sum(is.na(tune))==1) {
+ stop("Please specify a tune parameter and call MCMCoprobitChange() again.\n")
+ }
+ else if (length(tune)==1){
+ tune <- rep(tune, ns)
+ }
+ else if(length(tune)>1&length(tune)<ns){
+ tune <- rep(tune[1], ns)
+ cat("The first element of tune is repeated to make it conformable to the number of states.\n")
+ }
+ else{
+
+ }
+
+ xint <- match("(Intercept)", colnames(X), nomatch = 0)
+ if (xint > 0) {
+ new.X <- X[, -xint, drop = FALSE]
+ }
+ else
+ warning("An intercept is needed and assumed in MCMCoprobitChange()\n.")
+ if (ncol(new.X) == 0) {
+ polr.out <- polr(ordered(Y) ~ 1)
+ }
+ else {
+ polr.out <- polr(ordered(Y) ~ new.X)
+ }
+
+ ## prior for transition matrix
+ A0 <- trans.mat.prior(m=m, n=N, a=a, b=b)
+
+ ## prior for beta error checking
+ if(is.null(dim(b0))) {
+ b0 <- b0 * matrix(1,K,1)
+ }
+ if((dim(b0)[1] != K) || (dim(b0)[2] != 1)) {
+ cat("N(b0,B0) prior b0 not conformable.\n")
+ stop("Please respecify and call MCMCoprobitChange() again.\n")
+ }
+ if(is.null(dim(B0))) {
+ B0 <- B0 * diag(K)
+ }
+ if((dim(B0)[1] != K) || (dim(B0)[2] != K)) {
+ cat("N(b0,B0) prior B0 not conformable.\n")
+ stop("Please respecify and call MCMCoprobitChange() again.\n")
+ }
+ marginal.likelihood <- match.arg(marginal.likelihood)
+ B0.eigenvalues <- eigen(B0)$values
+ if (isTRUE(all.equal(min(B0.eigenvalues), 0))){
+ if (marginal.likelihood != "none"){
+ warning("Cannot calculate marginal likelihood with improper prior\n")
+ marginal.likelihood <- "none"
+ }
+ }
+ chib <- 0
+ if (marginal.likelihood == "Chib95"){
+ chib <- 1
+ }
+
+ ## to save time
+ B0inv <- solve(B0)
+ gamma.start <- matrix(NA, ncat + 1, 1)
+ gamma.start[1] <- -300
+ gamma.start[2] <- 0
+ gamma.start[3:ncat] <- (polr.out$zeta[2:(ncat - 1)] - polr.out$zeta[1]) * 0.588
+ gamma.start[ncat + 1] <- 300
+
+ ## initial values
+ mle <- polr(Y ~ X[,-1])
+ beta <- matrix(rep(c(mle$zeta[1], coef(mle)), ns), ns, , byrow=TRUE)
+ ols <- lm(as.double(Y) ~ X-1)
+ betalinearstart <- matrix(rep(coef(ols), ns), ns, , byrow=TRUE)
+ P <- trans.mat.prior(m=m, n=N, a=0.9, b=0.1)
+ Sigmastart <- summary(ols)$sigma
+ if (gamma.fixed==1){
+ gamma <- gamma.start
+ gamma.storage <-rep(0.0, nstore*gk)
+ }
+ else {
+ gamma <- matrix(rep(gamma.start, ns), ns, , byrow=T)
+ gamma.storage <- rep(0.0, nstore*ns*gk)
+ }
+
+ ## call C++ code to draw sample
+ posterior <- .C("MCMCoprobitChange",
+ betaout = as.double(rep(0.0, nstore*ns*K)),
+ betalinearout = as.double(rep(0.0, nstore*ns*K)),
+ gammaout = as.double(gamma.storage),
+ 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),
+ Xdata = as.double(X),
+ Xrow = as.integer(nrow(X)),
+ Xcol = as.integer(ncol(X)),
+
+ m = as.integer(m),
+ ncat = as.integer(ncat),
+
+ burnin = as.integer(burnin),
+ mcmc = as.integer(mcmc),
+ thin = as.integer(thin),
+ verbose = as.integer(verbose),
+
+ tunedata = as.double(tune),
+ lecuyer=as.integer(lecuyer),
+ seedarray=as.integer(seed.array),
+ lecuyerstream=as.integer(lecuyer.stream),
+
+ betastart = as.double(beta),
+ betalinearstart = as.double(betalinearstart),
+ gammastart = as.double(gamma),
+ Pstart = as.double(P),
+ sigmastart = as.double(Sigmastart),
+
+ a = as.double(a),
+ b = as.double(b),
+ b0data = as.double(b0),
+ B0data = as.double(B0),
+ A0data = as.double(A0),
+ logmarglikeholder = as.double(0.0),
+ loglikeholder = as.double(0.0),
+ chib = as.integer(chib),
+ gammafixed= as.integer(gamma.fixed))
+
+ ## get marginal likelihood if Chib95
+ if (chib==1){
+ logmarglike <- posterior$logmarglikeholder
+ loglike <- posterior$loglikeholder
+ }
+ else{
+ logmarglike <- loglike <- 0
+ }
+
+ ## pull together matrix and build MCMC object to return
+ beta.holder <- mcmc(matrix(posterior$betaout, nstore, ns*K))
+ if (gamma.fixed==1){
+ gamma.holder <- mcmc(matrix(posterior$gammaout, nstore, gk))
+ }
+ else {
+ gamma.holder <- mcmc(matrix(posterior$gammaout, nstore, ns*gk))
+ }
+ P.holder <- matrix(posterior$Pout, nstore, )
+ s.holder <- matrix(posterior$sout, nstore, )
+ ps.holder <- matrix(posterior$psout, N, )
+
+ varnames(beta.holder) <- sapply(c(1:ns),
+ function(i){
+ paste(c(xnames), "_regime", i, sep = "")
+ })
+ ## betalinear
+ betalinear.holder <- mcmc(matrix(posterior$betalinearout, nstore, ns*K))
+ varnames(betalinear.holder) <- sapply(c(1:ns),
+ function(i){
+ paste(c(xnames), "_regime", i, sep = "")
+ })
+ gamma.holder <- gamma.holder[, as.vector(sapply(1:ns, function(i){gk*(i-1) + (3:(gk-1))}))]
+ gamma.names <- paste("gamma", 3:(gk-1), sep="")
+ varnames(gamma.holder) <- sapply(c(1:ns),
+ function(i){
+ paste(gamma.names, "_regime", i, sep = "")
+ })
+
+ output <- mcmc(cbind(beta.holder, gamma.holder))
+ attr(output, "title") <- "MCMCoprobitChange Posterior Sample"
+ ## attr(output, "betalinear") <- mcmc(betalinear.holder)
+ attr(output, "formula") <- formula
+ attr(output, "y") <- Y
+ attr(output, "X") <- X
+ attr(output, "m") <- m
+ attr(output, "call") <- cl
+ attr(output, "logmarglike") <- logmarglike
+ attr(output, "loglike") <- loglike
+ attr(output, "prob.state") <- ps.holder/nstore
+ attr(output, "s.store") <- s.holder
+ return(output)
+
+ }## end of MCMC function
diff --git a/R/MCMCpoissonChange.R b/R/MCMCpoissonChange.R
index 4dc90e5..e525757 100644
--- a/R/MCMCpoissonChange.R
+++ b/R/MCMCpoissonChange.R
@@ -104,6 +104,7 @@
B0data = as.double(B0),
A0data = as.double(A0),
logmarglikeholder = as.double(0.0),
+ loglikeholder = as.double(0.0),
wrin = as.double(wr),
mrin = as.double(mr),
srin = as.double(sr),
@@ -114,6 +115,7 @@
## get marginal likelihood if Chib95
if (marginal.likelihood == "Chib95"){
logmarglike <- posterior$logmarglikeholder
+ loglike <- posterior$loglikeholder
}
## pull together matrix and build MCMC object to return
@@ -131,6 +133,7 @@
attr(output, "m") <- m
attr(output, "call") <- cl
attr(output, "logmarglike") <- logmarglike
+ attr(output, "loglike") <- loglike
attr(output, "prob.state") <- ps.holder/nstore
attr(output, "s.store") <- s.holder
attr(output, "P.store") <- P.holder
diff --git a/R/MCMCprobit.R b/R/MCMCprobit.R
index edb261c..a43b685 100644
--- a/R/MCMCprobit.R
+++ b/R/MCMCprobit.R
@@ -21,7 +21,7 @@
function(formula, data=NULL, burnin = 1000, mcmc = 10000,
thin = 1, verbose = 0, seed = NA, beta.start = NA,
b0 = 0, B0 = 0, bayes.resid=FALSE,
- marginal.likelihood = c("none", "Laplace"), ...) {
+ marginal.likelihood = c("none", "Laplace", "Chib95"), ...) {
## checks
check.offset(list(...))
@@ -84,24 +84,12 @@
stop("Check data and call MCMCprobit() again.\n")
}
-
- ## marginal likelihood calculation if Laplace
- if (marginal.likelihood == "Laplace"){
- theta.start <- beta.start
- optim.out <- optim(theta.start, logpost.probit, method="BFGS",
- control=list(fnscale=-1),
- hessian=TRUE, y=Y, X=X, b0=b0, B0=B0)
-
- theta.tilde <- optim.out$par
- beta.tilde <- theta.tilde[1:K]
-
- Sigma.tilde <- solve(-1*optim.out$hessian)
-
- logmarglike <- (length(theta.tilde)/2)*log(2*pi) +
- log(sqrt(det(Sigma.tilde))) +
- logpost.probit(theta.tilde, Y, X, b0, B0)
-
+ ## if Chib95 is true
+ chib <- 0
+ if (marginal.likelihood == "Chib95"){
+ chib <- 1
}
+
posterior <- NULL
if (is.null(resvec)){
@@ -117,18 +105,42 @@
seedarray=as.integer(seed.array),
lecuyerstream=as.integer(lecuyer.stream),
verbose=as.integer(verbose), betastart=beta.start,
- b0=b0, B0=B0)
+ b0=b0, B0=B0, logmarglikeholder.nonconst = as.double(0.0),
+ chib = as.integer(chib))
+
+ ## get marginal likelihood if Chib95
+ if (marginal.likelihood == "Chib95"){
+ logmarglike <- posterior$logmarglikeholder
+ }
+ ## marginal likelihood calculation if Laplace
+ if (marginal.likelihood == "Laplace"){
+ theta.start <- beta.start
+ optim.out <- optim(theta.start, logpost.probit, method="BFGS",
+ control=list(fnscale=-1),
+ hessian=TRUE, y=Y, X=X, b0=b0, B0=B0)
+
+ theta.tilde <- optim.out$par
+ beta.tilde <- theta.tilde[1:K]
+
+ Sigma.tilde <- solve(-1*optim.out$hessian)
+
+ logmarglike <- (length(theta.tilde)/2)*log(2*pi) +
+ log(sqrt(det(Sigma.tilde))) +
+ logpost.probit(theta.tilde, Y, X, b0, B0)
+
+ }
## put together matrix and build MCMC object to return
output <- form.mcmc.object(posterior, names=xnames,
title="MCMCprobit Posterior Sample",
- y=Y, call=cl, logmarglike=logmarglike)
+ y=Y, call=cl,
+ logmarglike=logmarglike)
}
else{
# define holder for posterior density sample
sample <- matrix(data=0, mcmc/thin, dim(X)[2]+length(resvec) )
-
+
## call C++ code to draw sample
auto.Scythe.call(output.object="posterior", cc.fun.name="MCMCprobitres",
sample.nonconst=sample, Y=Y, X=X, resvec=resvec,
@@ -138,7 +150,30 @@
seedarray=as.integer(seed.array),
lecuyerstream=as.integer(lecuyer.stream),
verbose=as.integer(verbose), betastart=beta.start,
- b0=b0, B0=B0)
+ b0=b0, B0=B0, logmarglikeholder.nonconst= as.double(0.0),
+ chib = as.integer(chib))
+
+ ## get marginal likelihood if Chib95
+ if (marginal.likelihood == "Chib95"){
+ logmarglike <- posterior$logmarglikeholder
+ }
+ ## marginal likelihood calculation if Laplace
+ if (marginal.likelihood == "Laplace"){
+ theta.start <- beta.start
+ optim.out <- optim(theta.start, logpost.probit, method="BFGS",
+ control=list(fnscale=-1),
+ hessian=TRUE, y=Y, X=X, b0=b0, B0=B0)
+
+ theta.tilde <- optim.out$par
+ beta.tilde <- theta.tilde[1:K]
+
+ Sigma.tilde <- solve(-1*optim.out$hessian)
+
+ logmarglike <- (length(theta.tilde)/2)*log(2*pi) +
+ log(sqrt(det(Sigma.tilde))) +
+ logpost.probit(theta.tilde, Y, X, b0, B0)
+
+ }
## put together matrix and build MCMC object to return
xnames <- c(xnames, paste("epsilonstar", as.character(resvec), sep="") )
diff --git a/R/MCMCprobitChange.R b/R/MCMCprobitChange.R
new file mode 100644
index 0000000..0eb8fdf
--- /dev/null
+++ b/R/MCMCprobitChange.R
@@ -0,0 +1,123 @@
+#########################################################
+##
+## sample from the posterior distribution
+## of a probit regression model with multiple changepoints
+##
+## JHP 07/01/2007
+## JHP 03/03/2009
+#########################################################
+
+"MCMCprobitChange"<-
+ function(formula, data=parent.frame(), m = 1,
+ burnin = 10000, mcmc = 10000, thin = 1, verbose = 0,
+ seed = NA, beta.start = NA, P.start = NA,
+ b0 = NULL, B0 = NULL, a = NULL, b = NULL,
+ marginal.likelihood = c("none", "Chib95"),
+ ...){
+
+ ## form response and model matrices
+ holder <- parse.formula(formula, data)
+ y <- holder[[1]]
+ X <- holder[[2]]
+ xnames <- holder[[3]]
+ k <- ncol(X)
+ n <- length(y)
+ ns <- m + 1
+
+ ## check iteration parameters
+ check.mcmc.parameters(burnin, mcmc, thin)
+ totiter <- mcmc + burnin
+ cl <- match.call()
+ nstore <- mcmc/thin
+
+ ## seeds
+ seeds <- form.seeds(seed)
+ lecuyer <- seeds[[1]]
+ seed.array <- seeds[[2]]
+ lecuyer.stream <- seeds[[3]]
+ if(!is.na(seed)) set.seed(seed)
+
+ ## prior
+ mvn.prior <- form.mvn.prior(b0, B0, k)
+ b0 <- mvn.prior[[1]]
+ B0 <- mvn.prior[[2]]
+ A0 <- trans.mat.prior(m=m, n=n, a=a, b=b)
+
+
+ ## initial values
+ Pstart <- check.P(P.start, m, a=a, b=b)
+ betastart <- beta.change.start(beta.start, ns, k, formula, family=binomial, data)
+
+ chib <- 0
+ if (marginal.likelihood == "Chib95"){
+ chib <- 1
+ }
+
+ ## call C++ code to draw sample
+ posterior <- .C("MCMCprobitChange",
+ betaout = as.double(rep(0.0, nstore*ns*k)),
+ 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)),
+ Xdata = as.double(X),
+ Xrow = as.integer(nrow(X)),
+ Xcol = as.integer(ncol(X)),
+ m = as.integer(m),
+ burnin = as.integer(burnin),
+ mcmc = as.integer(mcmc),
+ thin = as.integer(thin),
+ verbose = as.integer(verbose),
+
+ lecuyer=as.integer(lecuyer),
+ seedarray=as.integer(seed.array),
+ lecuyerstream=as.integer(lecuyer.stream),
+
+ betastart = as.double(betastart),
+ Pstart = as.double(Pstart),
+
+ a = as.double(a),
+ b = as.double(b),
+ b0data = as.double(b0),
+ B0data = as.double(B0),
+ A0data = as.double(A0),
+ logmarglikeholder = as.double(0.0),
+ loglikeholder = as.double(0.0),
+ chib = as.integer(chib))
+
+ ## get marginal likelihood if Chib95
+ if (chib==1){
+ logmarglike <- posterior$logmarglikeholder
+ loglike <- posterior$loglikeholder
+ }
+ else{
+ logmarglike <- loglike <- 0
+ }
+
+ ## pull together matrix and build MCMC object to return
+ beta.holder <- matrix(posterior$betaout, nstore, ns*k)
+ P.holder <- matrix(posterior$Pout, nstore, )
+ s.holder <- matrix(posterior$sout, nstore, )
+ ps.holder <- matrix(posterior$psout, n, )
+
+ output <- mcmc(data=beta.holder, start=burnin+1, end=burnin + mcmc, thin=thin)
+ varnames(output) <- sapply(c(1:ns),
+ function(i){
+ paste(c(xnames), "_regime", i, sep = "")
+ })
+ attr(output, "title") <- "MCMCprobitChange Posterior Sample"
+ attr(output, "formula") <- formula
+ attr(output, "y") <- y
+ attr(output, "X") <- X
+ attr(output, "m") <- m
+ attr(output, "call") <- cl
+ attr(output, "logmarglike") <- logmarglike
+ attr(output, "loglike") <- loglike
+ attr(output, "prob.state") <- ps.holder/nstore
+ attr(output, "s.store") <- s.holder
+ return(output)
+
+}## end of MCMC function
+
diff --git a/man/MCMCmetrop1R.Rd b/man/MCMCmetrop1R.Rd
index a933b3e..30d5f5d 100644
--- a/man/MCMCmetrop1R.Rd
+++ b/man/MCMCmetrop1R.Rd
@@ -202,6 +202,9 @@ mode. This last calculation is done via an initial call to
}
}
\references{
+ Siddhartha Chib; Edward Greenberg. 1995. ``Understanding the Metropolis-Hastings Algorithm."
+ \emph{The American Statistician}, 49, 327-335.
+
Andrew Gelman, John B. Carlin, Hal S. Stern, and Donald
B. Rubin. 2003. \emph{Bayesian Data Analysis}. 2nd Edition. Boca
Raton: Chapman & Hall/CRC.
diff --git a/man/MCMCmnl.Rd b/man/MCMCmnl.Rd
index 9a8eec8..418a40f 100644
--- a/man/MCMCmnl.Rd
+++ b/man/MCMCmnl.Rd
@@ -161,6 +161,9 @@ that can be used to analyze the posterior sample.
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/}.
+
+ Siddhartha Chib, Edward Greenberg, and Yuxin Chen. 1998.
+ ``MCMC Methods for Fitting and Comparing Multinomial Response Models."
}
\seealso{\code{\link[coda]{plot.mcmc}},\code{\link[coda]{summary.mcmc}},
diff --git a/man/MCMCoprobit.Rd b/man/MCMCoprobit.Rd
index faceaab..1f15640 100644
--- a/man/MCMCoprobit.Rd
+++ b/man/MCMCoprobit.Rd
@@ -9,7 +9,9 @@
object, which can be subsequently analyzed with functions
provided in the coda package.
}
-
+
+
+
\usage{
MCMCoprobit(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
thin=1, tune = NA, tdf = 1, verbose = 0, seed = NA, beta.start = NA,
@@ -62,7 +64,7 @@ MCMCoprobit(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
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
@@ -121,7 +123,8 @@ that can be used to analyze the posterior sample.
}
\references{
- Albert, James and Siddhartha Chib. 2001. ``Sequential Ordinal Modeling with Applications to Survival Data." \emph{Biometrics.} 57: 829-836.
+ Albert, J. H. and S. Chib. 1993. ``Bayesian Analysis of Binary and
+ Polychotomous Response Data.'' \emph{J. Amer. Statist. Assoc.} 88, 669-679
M. K. Cowles. 1996. ``Accelerating Monte Carlo Markov Chain Convergence for
Cumulative-link Generalized Linear Models." \emph{Statistics and Computing.}
@@ -129,6 +132,8 @@ that can be used to analyze the posterior sample.
Valen E. Johnson and James H. Albert. 1999. \emph{Ordinal Data Modeling}.
Springer: New York.
+
+ Albert, James and Siddhartha Chib. 2001. ``Sequential Ordinal Modeling with Applications to Survival Data." \emph{Biometrics.} 57: 829-836.
Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007.
\emph{Scythe Statistical Library 1.0.} \url{http://scythe.wustl.edu}.
diff --git a/man/MCMCoprobitChange.Rd b/man/MCMCoprobitChange.Rd
new file mode 100644
index 0000000..9bdedec
--- /dev/null
+++ b/man/MCMCoprobitChange.Rd
@@ -0,0 +1,183 @@
+\name{MCMCoprobitChange}
+\alias{MCMCoprobitChange}
+
+\title{Markov Chain Monte Carlo for Ordered Probit Changepoint Regression Model}
+\description{
+ This function generates a sample from the posterior distribution
+ of an ordered probit regression model with multiple parameter breaks. 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{MCMCoprobitChange(formula, data=parent.frame(), m = 1,
+ burnin = 1000, mcmc = 1000, thin = 1, tune = NA, verbose = 0,
+ seed = NA, beta.start = NA, gamma.start=NA, P.start = NA,
+ b0 = NULL, B0 = NULL, a = NULL, b = NULL,
+ marginal.likelihood = c("none", "Chib95"), gamma.fixed=0, ...)}
+
+\arguments{
+ \item{formula}{Model formula.}
+
+ \item{data}{Data frame.}
+
+ \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 burnin.}
+
+ \item{thin}{The thinning interval used in the simulation. The number of
+ MCMC iterations must be divisible by this value.}
+
+ \item{tune}{The tuning parameter for the Metropolis-Hastings
+ step. Default of NA corresponds to a choice of 0.05 divided by the
+ number of categories in the response variable.}
+
+ \item{verbose}{A switch which determines whether or not the progress of
+ the sampler is printed to the screen. If \code{verbose} is greater
+ than 0 the iteration number, the \eqn{\beta}{beta} vector, and the error variance are printed to
+ the screen every \code{verbose}th iteration.}
+
+ \item{seed}{The seed for the random number generator. If NA, the Mersenne
+ Twister generator is used with default seed 12345; if an integer is
+ 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 MLE
+ 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{gamma.start}{The starting values for the \eqn{\gamma}{gamma} vector.
+ This can either be a scalar or a
+ column vector with dimension equal to the number of gammas.
+ The default value of of NA will use the MLE
+ estimate of \eqn{\gamma}{gamma} as the starting value. If this is a
+ scalar, that value will serve as the starting value
+ mean for all of the gammas.}
+
+ \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{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{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{gamma.fixed}{1 if users want to constrain \eqn{\gamma}{gamma} values to be constant. By default,
+ \eqn{\gamma}{gamma} values are allowed to vary across regimes.}
+
+ \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, the log-likelihood of the model (\code{loglike}), and
+ the log-marginal likelihood of the model (\code{logmarglike}).
+}
+
+\details{
+ \code{MCMCoprobitChange} simulates from the posterior distribution of
+ an ordinal probit regression model with multiple parameter breaks. The simulation of latent states is based on
+ the linear approximation method discussed in Park (2011).
+
+ The model takes the following form:
+ \deqn{\Pr(y_t = 1) = \Phi(\gamma_{c, m} - x_i'\beta_m) - \Phi(\gamma_{c-1, m} - x_i'\beta_m)\;\; m = 1, \ldots, M}{
+ Pr(y_t = 1) = Phi(gamma_(c, m) - x_i'beta_m) - Phi(gamma_(c-1, m) - x_i'beta)
+ }
+ Where \eqn{M}{M} is the number of states, and \eqn{\gamma_{c, m}}{gamma_(c, m)} and \eqn{\beta_m}{beta_m}
+ are paramters when a state is \eqn{m}{m} at \eqn{t}{t}.
+
+ Note that when the fitted changepoint model has very few observations in any of states,
+ the marginal likelihood outcome can be ``nan," which indicates that too many breaks are assumed
+ given the model and data.
+
+}
+
+\references{
+
+ Jong Hee Park. 2011. ``Changepoint Analysis of Binary and Ordinal Probit Models:
+ An Application to Bank Rate Policy Under the Interwar Gold Standard." \emph{Political Analysis}.
+
+ Siddhartha Chib. 1998. "Estimation and comparison of multiple change-point models."
+ \emph{Journal of Econometrics}. 86: 221-241.
+
+}
+
+\examples{
+set.seed(1909)
+N <- 200
+x1 <- rnorm(N, 1, .5);
+
+## set a true break at 100
+z1 <- 1 + x1[1:100] + rnorm(100);
+z2 <- 1 -0.2*x1[101:200] + rnorm(100);
+z <- c(z1, z2);
+y <- z
+
+## generate y
+y[z < 1] <- 1;
+y[z >= 1 & z < 2] <- 2;
+y[z >= 2] <- 3;
+
+## inputs
+formula <- y ~ x1
+
+## fit multiple models with a varying number of breaks
+out1 <- MCMCoprobitChange(formula, m=1,
+ mcmc=1000, burnin=1000, thin=1, tune=c(.5, .5), verbose=1000,
+ b0=0, B0=10, marginal.likelihood = "Chib95")
+out2 <- MCMCoprobitChange(formula, m=2,
+ mcmc=1000, burnin=1000, thin=1, tune=c(.5, .5, .5), verbose=1000,
+ b0=0, B0=10, marginal.likelihood = "Chib95")
+out3 <- MCMCoprobitChange(formula, m=3,
+ mcmc=1000, burnin=1000, thin=1, tune=c(.5, .5, .5, .5), verbose=1000,
+ b0=0, B0=10, marginal.likelihood = "Chib95")
+
+## find the most reasonable one
+BayesFactor(out1, out2, out3)
+
+## draw plots using the "right" model
+plotState(out1)
+plotChangepoint(out1)
+}
+
+\keyword{models}
+
+\seealso{\code{\link{plotState}}, \code{\link{plotChangepoint}}}
diff --git a/man/MCMCpoissonChange.Rd b/man/MCMCpoissonChange.Rd
index 697fedd..63f8843 100644
--- a/man/MCMCpoissonChange.Rd
+++ b/man/MCMCpoissonChange.Rd
@@ -87,7 +87,8 @@
\details{
\code{MCMCpoissonChange} simulates from the posterior distribution of
- a Poisson regression model with multiple changepoints.
+ a Poisson regression model with multiple changepoints using the methods of Chib (1998) and Fruhwirth-Schnatter and Wagner (2006).
+ The details of the model are discussed in Park (2010).
The model takes the following form:
\deqn{y_t \sim \mathcal{P}oisson(\mu_t)}{y_t ~ Poisson(mu_t)}
@@ -104,13 +105,18 @@
\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.
+ Jong Hee Park. 2010. "Structural Change in the U.S. Presidents' Use of Force Abroad."
+ \emph{American Journal of Political Science} 54: 766-782.
+
+ Sylvia Fruhwirth-Schnatter and Helga Wagner 2006. "Auxiliary Mixture Sampling for Parameter-driven Models of
+ Time Series of Counts with Applications to State Space Modelling." \emph{Biometrika}. 93:827--841.
Siddhartha Chib. 1998. "Estimation and comparison of multiple change-point models."
\emph{Journal of Econometrics}. 86: 221-241.
+ Siddhartha Chib. 1995. "Marginal Likelihood from the Gibbs Output."
+ \emph{Journal of the American Statistical Association}. 90:
+ 1313-1321.
}
\examples{
@@ -122,7 +128,7 @@
true.beta2 <- c(1, -2)
true.beta3 <- c(1, 2)
- ## true two breaks at (50, 100)
+ ## set true two breaks at (50, 100)
true.s <- rep(1:3, each=n/3)
mu1 <- exp(1 + x1[true.s==1]*1)
mu2 <- exp(1 + x1[true.s==2]*-2)
@@ -131,6 +137,7 @@
y <- as.ts(c(rpois(n/3, mu1), rpois(n/3, mu2), rpois(n/3, mu3)))
formula = y ~ x1
+ ## fit multiple models with a varying number of breaks
model1 <- MCMCpoissonChange(formula, m=1,
mcmc = 1000, burnin = 1000, verbose = 500,
b0 = rep(0, 2), B0 = 5*diag(2), marginal.likelihood = "Chib95")
@@ -147,14 +154,15 @@
mcmc = 1000, burnin = 1000, verbose = 500,
b0 = rep(0, 2), B0 = 5*diag(2), marginal.likelihood = "Chib95")
+ ## find the most reasonable one
print(BayesFactor(model1, model2, model3, model4, model5))
- ## Draw plots using the "right" model
+ ## draw plots using the "right" model
par(mfrow=c(attr(model2, "m") + 1, 1), mai=c(0.4, 0.6, 0.3, 0.05))
plotState(model2, legend.control = c(1, 0.6))
plotChangepoint(model2, verbose = TRUE, ylab="Density", start=1, overlay=TRUE)
- ## No covariate
+ ## No covariate case
model2.1 <- MCMCpoissonChange(y ~ 1, m = 2, c0 = 2, d0 = 1,
mcmc = 1000, burnin = 1000, verbose = 500,
marginal.likelihood = "Chib95")
diff --git a/man/MCMCprobit.Rd b/man/MCMCprobit.Rd
index fbdb03c..d8b3912 100644
--- a/man/MCMCprobit.Rd
+++ b/man/MCMCprobit.Rd
@@ -14,7 +14,7 @@
MCMCprobit(formula, data = NULL, burnin = 1000, mcmc = 10000,
thin = 1, verbose = 0, seed = NA, beta.start = NA,
b0 = 0, B0 = 0, bayes.resid = FALSE,
- marginal.likelihood=c("none", "Laplace"), ...) }
+ marginal.likelihood=c("none", "Laplace", "Chib95"), ...) }
\arguments{
\item{formula}{Model formula.}
@@ -75,8 +75,9 @@ MCMCprobit(formula, data = NULL, burnin = 1000, mcmc = 10000,
\item{marginal.likelihood}{How should the marginal likelihood be
calculated? Options are: \code{none} in which case the marginal
- likelihood will not be calculated or \code{Laplace} in which case the
- Laplace approximation (see Kass and Raftery, 1995) is used.}
+ likelihood will not be calculated, \code{Laplace} in which case the
+ Laplace approximation (see Kass and Raftery, 1995) is used, or \code{Chib95}
+ in which case Chib (1995) method is used.}
\item{...}{further arguments to be passed}
}
@@ -104,12 +105,16 @@ MCMCprobit(formula, data = NULL, burnin = 1000, mcmc = 10000,
}
\references{
- Albert, J. H. and S. Chib. 1993. ``Bayesian Analysis of Binary and
- Polychotomous Response Data.'' \emph{J. Amer. Statist. Assoc.} 88, 669-679
+ Albert, J. H. and S. Chib. 1993. ``Bayesian Analysis of Binary and
+ Polychotomous Response Data.'' \emph{J. Amer. Statist. Assoc.} 88, 669-679
- Albert, J. H. and S. Chib. 1995. ``Bayesian Residual Analysis for
- Binary Response Regression Models.'' \emph{Biometrika.} 82, 747-759.
+ Albert, J. H. and S. Chib. 1995. ``Bayesian Residual Analysis for
+ Binary Response Regression Models.'' \emph{Biometrika.} 82, 747-759.
+ Siddhartha Chib. 1995. "Marginal Likelihood from the Gibbs Output."
+ \emph{Journal of the American Statistical Association}. 90:
+ 1313-1321.
+
Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007.
\emph{Scythe Statistical Library 1.0.} \url{http://scythe.wustl.edu}.
@@ -122,9 +127,15 @@ MCMCprobit(formula, data = NULL, burnin = 1000, mcmc = 10000,
\examples{
\dontrun{
data(birthwt)
- posterior <- MCMCprobit(low~age+as.factor(race)+smoke, data=birthwt)
- plot(posterior)
- summary(posterior)
+ out1 <- MCMCprobit(low~as.factor(race)+smoke, data=birthwt,
+ b0 = 0, B0 = 10, marginal.likelihood="Chib95")
+ out2 <- MCMCprobit(low~age+as.factor(race), data=birthwt,
+ b0 = 0, B0 = 10, marginal.likelihood="Chib95")
+ out3 <- MCMCprobit(low~age+as.factor(race)+smoke, data=birthwt,
+ b0 = 0, B0 = 10, marginal.likelihood="Chib95")
+ BayesFactor(out1, out2, out3)
+ plot(out3)
+ summary(out3)
}
}
diff --git a/man/MCMCprobitChange.Rd b/man/MCMCprobitChange.Rd
new file mode 100644
index 0000000..139e9ff
--- /dev/null
+++ b/man/MCMCprobitChange.Rd
@@ -0,0 +1,167 @@
+\name{MCMCprobitChange}
+\alias{MCMCprobitChange}
+
+\title{Markov Chain Monte Carlo for a linear Gaussian Multiple Changepoint Model}
+\description{
+ This function generates a sample from the posterior distribution
+ of a linear Gaussian 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{MCMCprobitChange(formula, data=parent.frame(), m = 1,
+ burnin = 10000, mcmc = 10000, thin = 1, verbose = 0,
+ seed = NA, beta.start = NA, P.start = NA,
+ b0 = NULL, B0 = NULL, a = NULL, b = NULL,
+ marginal.likelihood = c("none", "Chib95"), ...)}
+
+\arguments{
+ \item{formula}{Model formula.}
+
+ \item{data}{Data frame.}
+
+ \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 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 are printed to
+ the screen every \code{verbose}th iteration.}
+
+ \item{seed}{The seed for the random number generator. If NA, the Mersenne
+ Twister generator is used with default seed 12345; if an integer is
+ 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 MLE
+ 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{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{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{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, the
+ log-likelihood of the model (\code{loglike}), and
+ the log-marginal likelihood of the model (\code{logmarglike}).
+}
+
+\details{
+ \code{MCMCprobitChange} simulates from the posterior distribution of
+ a probit regression model with multiple parameter breaks. The simulation is based on Chib (1998) and Park (2011).
+
+ The model takes the following form:
+ \deqn{\Pr(y_t = 1) = \Phi(x_i'\beta_m) \;\; m = 1, \ldots, M}{
+ Pr(y_t = 1) = Phi(x_i'beta_m)}
+ Where \eqn{M}{M} is the number of states, and \eqn{\beta_m}{beta_m}
+ is a parameter when a state is \eqn{m}{m} at \eqn{t}{t}.
+}
+
+\author{Jong Hee Park, \email{jhp at uchicago.edu}, \url{http://home.uchicago.edu/~jhp/}.}
+
+\references{
+ Jong Hee Park. 2011. ``Changepoint Analysis of Binary and Ordinal Probit Models:
+ An Application to Bank Rate Policy Under the Interwar Gold Standard." \emph{Political Analysis}.
+
+ Siddhartha Chib. 1998. "Estimation and comparison of multiple change-point models."
+ \emph{Journal of Econometrics}. 86: 221-241.
+
+ Albert, J. H. and S. Chib. 1993. ``Bayesian Analysis of Binary and
+ Polychotomous Response Data.'' \emph{J. Amer. Statist. Assoc.} 88, 669-679
+}
+
+\examples{
+set.seed(1973)
+x1 <- rnorm(300, 0, 1)
+true.beta <- c(-.5, .2, 1)
+true.alpha <- c(.1, -1., .2)
+X <- cbind(1, x1)
+
+## set two true breaks at 100 and 200
+true.phi1 <- pnorm(true.alpha[1] + x1[1:100]*true.beta[1])
+true.phi2 <- pnorm(true.alpha[2] + x1[101:200]*true.beta[2])
+true.phi3 <- pnorm(true.alpha[3] + x1[201:300]*true.beta[3])
+
+## generate y
+y1 <- rbinom(100, 1, true.phi1)
+y2 <- rbinom(100, 1, true.phi2)
+y3 <- rbinom(100, 1, true.phi3)
+Y <- as.ts(c(y1, y2, y3))
+
+## fit multiple models with a varying number of breaks
+out0 <- MCMCprobit(formula=Y~X-1, data=parent.frame(),
+ mcmc=1000, burnin=1000, thin=1, verbose=1000,
+ b0 = 0, B0 = 10, marginal.likelihood = c("Laplace"))
+out1 <- MCMCprobitChange(formula=Y~X-1, data=parent.frame(), m=1,
+ mcmc=1000, burnin=1000, thin=1, verbose=1000,
+ b0 = 0, B0 = 10, a = 1, b = 1, marginal.likelihood = c("Chib95"))
+out2 <- MCMCprobitChange(formula=Y~X-1, data=parent.frame(), m=2,
+ mcmc=1000, burnin=1000, thin=1, verbose=1000,
+ b0 = 0, B0 = 10, a = 1, b = 1, marginal.likelihood = c("Chib95"))
+out3 <- MCMCprobitChange(formula=Y~X-1, data=parent.frame(), m=3,
+ mcmc=1000, burnin=1000, thin=1, verbose=1000,
+ b0 = 0, B0 = 10, a = 1, b = 1, marginal.likelihood = c("Chib95"))
+
+## find the most reasonable one
+BayesFactor(out0, out1, out2, out3)
+
+## draw plots using the "right" model
+plotState(out2)
+plotChangepoint(out2)
+}
+
+\keyword{models}
+
+\seealso{\code{\link{plotState}}, \code{\link{plotChangepoint}}}
diff --git a/man/MCMCtobit.Rd b/man/MCMCtobit.Rd
index 9494768..2b27b51 100644
--- a/man/MCMCtobit.Rd
+++ b/man/MCMCtobit.Rd
@@ -130,7 +130,10 @@ MCMCtobit(formula, data = NULL, below = 0, above = Inf,
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/}.
-
+
+ Siddhartha Chib. 1992. ``Bayes inference in the Tobit censored regression model."
+ \emph{Journal of Econometrics}. 51:79-99.
+
James Tobin. 1958. ``Estimation of relationships for limited dependent
variables." \emph{Econometrica.} 26:24-36.
}
diff --git a/src/MCMCoprobitChange.cc b/src/MCMCoprobitChange.cc
new file mode 100644
index 0000000..14ec4cb
--- /dev/null
+++ b/src/MCMCoprobitChange.cc
@@ -0,0 +1,747 @@
+// MCMCoprobitChange.cc is C++ code to estimate a oprobit changepoint model
+// with linear approximation
+//
+// Jong Hee Park
+// Dept. of Political Science
+// University of Chicago
+// jhp at uchicago.edu
+//
+// 07/06/2007
+// 11/02/2009 Modified
+
+#ifndef MCMCOPROBITCHANGE_CC
+#define MCMCOPROBITCHANGE_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>
+#include <R_ext/Utils.h>
+
+
+using namespace std;
+using namespace scythe;
+
+// density function for truncated normal
+static double dtnormLX(const double x,
+ const double mean,
+ const double sd,
+ const double lower,
+ const double upper){
+ double out = 0.0;
+ if (x>lower&x<upper){
+ double numer = dnorm(x, mean, sd);
+ double denom = pnorm(upper, mean, sd) - pnorm(lower, mean, sd);
+ out = numer/denom;
+ }
+ else{
+ Rprintf("\n x input for dtnormLX() %10.5f is out of bounds %10.5f %10.5f ", x, lower, upper, "\n");
+ out = 1;
+ }
+ return out;
+}
+
+// likelihood of oprobit
+static double oprobit_pdfLX(const int ncat,
+ const int Y,
+ double Xbeta,
+ const Matrix<>& gamma){
+
+ Matrix<> cat_prob(1, ncat-1);
+ Matrix<> prob(1, ncat);
+ for (int j=0; j< ncat-1; ++j){
+ cat_prob(0, j) = pnorm(gamma[j+1] - Xbeta, 0.0, 1.0);
+ }
+ prob(0, ncat-1) = 1 - cat_prob(0, ncat-2);
+ prob(0, 0) = cat_prob(0, 0);
+ for (int j=1; j<(ncat-1); ++j){
+ prob(0, j) = cat_prob(0,j) - cat_prob(0, j-1);
+ }
+ double like = prob(0,Y-1);
+
+ return like;
+}
+
+static double oprobit_log_postLX(unsigned int j,
+ const int ncat,
+ const Matrix<>& gamma_p,
+ const Matrix<>& gamma,
+ const Matrix<>& Y,
+ const Matrix<>& X,
+ const Matrix<>& beta,
+ const Matrix<>& tune,
+ const int gammafixed){
+ const int N = Y.rows();
+ double loglikerat = 0.0;
+ double loggendenrat = 0.0;
+ Matrix<> Xbeta = X*t(beta(j,_));
+
+ if (gammafixed==1){
+ for (unsigned int i=0; i<N; ++i){
+ int yi = Y[i];
+ if (yi == ncat){
+ loglikerat = loglikerat
+ + log(1.0 - pnorm(gamma_p(yi-1) - Xbeta[i], 0.0, 1.0) )
+ - log(1.0 - pnorm(gamma(yi-1) - Xbeta[i], 0.0, 1.0) );
+ }
+ else if (yi == 1){
+ loglikerat = loglikerat + log(pnorm(gamma_p(yi) - Xbeta[i], 0.0, 1.0) )
+ - log(pnorm(gamma(yi) - Xbeta[i], 0.0, 1.0) );
+ }
+ else{
+ loglikerat = loglikerat
+ + log(pnorm(gamma_p(yi) - Xbeta[i], 0.0, 1.0) -
+ pnorm(gamma_p(yi-1) - Xbeta[i], 0.0, 1.0) )
+ - log(pnorm(gamma(yi) - Xbeta[i], 0.0, 1.0) -
+ pnorm(gamma(yi - 1) - Xbeta[i], 0.0, 1.0) );
+ }
+ }
+ for (unsigned int k =2; k <ncat; ++k){
+ loggendenrat = loggendenrat +
+ log(pnorm(gamma(k+1), gamma(k), tune[j]) -
+ pnorm(gamma_p(k-1), gamma(k), tune[j])) -
+ log(pnorm(gamma_p(k+1), gamma_p(k), tune[j]) +
+ pnorm(gamma(k-1), gamma_p(k), tune[j]));
+ }
+ }
+ else{
+ for (unsigned int i=0; i<N; ++i){
+ int yi = Y[i];
+ if (yi == ncat){
+ loglikerat = loglikerat
+ + log(1.0 - pnorm(gamma_p(j, yi-1) - Xbeta[i], 0.0, 1.0) )
+ - log(1.0 - pnorm(gamma(j, yi-1) - Xbeta[i], 0.0, 1.0) );
+ }
+ else if (yi == 1){
+ loglikerat = loglikerat + log(pnorm(gamma_p(j, yi) - Xbeta[i], 0.0, 1.0) )
+ - log(pnorm(gamma(j, yi) - Xbeta[i], 0.0, 1.0) );
+ }
+ else{
+ loglikerat = loglikerat
+ + log(pnorm(gamma_p(j, yi) - Xbeta[i], 0.0, 1.0) -
+ pnorm(gamma_p(j, yi-1) - Xbeta[i], 0.0, 1.0) )
+ - log(pnorm(gamma(j, yi) - Xbeta[i], 0.0, 1.0) -
+ pnorm(gamma(j, yi - 1) - Xbeta[i], 0.0, 1.0) );
+ }
+ }
+ for (unsigned int k =2; k <ncat; ++k){
+ loggendenrat = loggendenrat +
+ log(pnorm(gamma(j ,k+1), gamma(j, k), tune[j]) -
+ pnorm(gamma_p(j, k-1), gamma(j, k), tune[j])) -
+ log(pnorm(gamma_p(j, k+1), gamma_p(j, k), tune[j]) +
+ pnorm(gamma(j, k-1), gamma_p(j, k), tune[j]));
+ }
+ }
+ return loglikerat + loggendenrat;
+}
+
+template <typename RNGTYPE>
+Matrix<> gaussian_ordinal_state_sampler_fixedsigma(rng<RNGTYPE>& stream,
+ const int m,
+ const Matrix<>& Y,
+ const Matrix<>& X,
+ const Matrix<>& beta,
+ const Matrix<>& Sigma,
+ const Matrix<>& 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);
+
+ for (int t=0; t<n ; ++t){
+ Matrix<> mu = X(t,_)*::t(beta);
+ for (int j = 0; j< ns; ++j){
+ py[j] = dnorm(Y[t], mu[j], sqrt(Sigma[0]));
+ }
+ if (t==0) pstyt1 = pr1;
+ else {
+ pstyt1 = ::t(F(t-1,_)*P);
+ }
+ Matrix<> unnorm_pstyt = pstyt1%py;
+ const Matrix<> pstyt = unnorm_pstyt/sum(unnorm_pstyt);
+ for (int j=0; j<ns ; ++j) F(t,j) = pstyt(j);
+
+ }
+
+ Matrix<int> s(n, 1);
+ Matrix<> ps = Matrix<>(n, ns);
+ ps(n-1,_) = F(n-1,_);
+ s(n-1) = ns;
+
+ Matrix<> pstyn = Matrix<>(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;
+ }
+ Matrix<> Sout(n, ns+1);
+ Sout(_, 0) = s(_,0);
+ for (int j = 0; j<ns; ++j){
+ Sout(_,j+1) = ps(_, j);
+ }
+ return Sout;
+}
+
+
+////////////////////////////////////////////
+// MCMCoprobitChangeLinearAppxpoint implementation.
+////////////////////////////////////////////
+template <typename RNGTYPE>
+void MCMCoprobitChange_impl(rng<RNGTYPE>& stream,
+ const int m, const int ncat,
+ const Matrix<>& Y, const Matrix<>& X,
+ Matrix<>& beta, Matrix<>& beta_linear, Matrix<>& gamma, Matrix<>& P,
+ Matrix<>& Sigma,
+ Matrix<>& b0, Matrix<>& B0, const Matrix<>& A0,
+ unsigned int burnin, unsigned int mcmc, unsigned int thin,
+ unsigned int verbose, const Matrix<>& tune, // const Matrix<int>& tdf,
+ bool chib, bool gammafixed,
+ Matrix<>& beta_store, Matrix<>& beta_linear_store,
+ Matrix<>& gamma_store, Matrix<>& Z_store,
+ Matrix<>& P_store, Matrix<>& ps_store, Matrix<int>& s_store,
+ double& logmarglike, double& loglike)
+{
+ const unsigned int tot_iter = burnin + mcmc;
+ const unsigned int nstore = mcmc / thin;
+ const unsigned int N = Y.rows();
+ const unsigned int ns = m + 1;
+ const unsigned int k = X.cols();
+ const unsigned int gk = ncat + 1;
+ const Matrix<> B0inv = invpd(B0);
+ Matrix<> Z(N, 1);
+ Matrix<> accepts(ns, 1);
+ Matrix<> gamma_p = gamma;
+
+ //MCMC loop
+ unsigned int count = 0;
+ for (int iter = 0; iter < tot_iter; ++iter){
+ // 1. Sample s
+ Matrix<> Sout = gaussian_ordinal_state_sampler_fixedsigma(stream, m, Y, X, beta_linear, Sigma, P);
+ Matrix<int> s = Sout(_, 0);
+ Matrix <double> ps(N, ns);
+ for (int j = 0; j<ns; ++j){
+ ps(_,j) = Sout(_,j+1);
+ }
+ // 2. Sample Z
+ for (unsigned int i = 0; i<N; ++i) {
+ Matrix<> mu = X(i,_)*t(beta(s[i]-1,_));
+ int yi = Y[i];
+ Z[i] = stream.rtnorm_combo(mu[0], 1.0, gamma(s[i]-1, yi-1), gamma(s[i]-1, yi));
+ }
+
+ // 3. Sample beta
+ int beta_count = 0;
+ Matrix<int> beta_count_storage(ns, 1);
+ Matrix<int> nstate(ns, 1);
+ for (int j = 0; j<ns; ++j){
+ for (int i = 0; i<N; ++i){
+ if (s[i] == j + 1) {
+ nstate[j] = nstate[j] + 1;
+ }
+ }
+ beta_count = beta_count + nstate[j];
+ Matrix<> yj = Y((beta_count - nstate[j]), 0, (beta_count - 1), 0);
+ Matrix<> Xj = X((beta_count - nstate[j]), 0, (beta_count - 1), k-1);
+ Matrix<> Zj = Z((beta_count - nstate[j]), 0, (beta_count - 1), 0);
+ Matrix<> XpX = t(Xj)*Xj;
+ Matrix<> XpZ = t(Xj)*Zj;
+ Matrix<> XpY = t(Xj)*yj;
+ Matrix<> Bn = invpd(B0inv + XpX/Sigma[0]);
+ Matrix<> bn = Bn*(B0inv*b0 + XpY/Sigma[0]);
+ beta_linear(j,_) = stream.rmvnorm(bn, Bn);
+ Matrix<> Bn2 = invpd(B0inv + XpX);
+ Matrix<> bn2 = Bn2*(B0inv*b0 + XpZ);
+ beta(j,_) = stream.rmvnorm(bn2, Bn2);
+ beta_count_storage[j] = beta_count;
+ }
+
+ // 4. Sample gamma
+ for (int j = 0; j < ns ; ++j){
+ for (int i = 2; i< ncat; ++i){
+ if (i==(ncat-1)){
+ gamma_p(j, i) = stream.rtbnorm_combo(gamma(j, i), ::pow(tune[j], 2.0),
+ gamma_p(j, i-1));
+ }
+ else {
+ gamma_p(j, i) = stream.rtnorm_combo(gamma(j, i), ::pow(tune[j], 2.0),
+ gamma_p(j, i-1), gamma(j, i+1));
+ }
+ }
+ Matrix<> Yj = Y((beta_count_storage[j] - nstate[j]), 0, (beta_count_storage[j] - 1), 0);
+ Matrix<> Xj = X((beta_count_storage[j] - nstate[j]), 0, (beta_count_storage[j] - 1), k-1);
+ double alpha = oprobit_log_postLX(j, ncat, gamma_p, gamma, Yj, Xj, beta, tune, gammafixed);
+
+ if (stream.runif() <= exp(alpha)){
+ gamma(j,_) = gamma_p(j,_);
+ accepts[j] = accepts[j] + 1;
+ }
+ }
+
+ // 5. Sample P
+ double shape1 = 0;
+ double shape2 = 0;
+ P(ns-1, ns-1) = 1;
+
+ for (int j =0; j< (ns-1); ++j){
+ shape1 = A0(j,j) + (double)nstate[j] - 1;
+ shape2 = A0(j,j+1) + 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)){
+ Matrix<> tbeta = ::t(beta);
+ Matrix<> tbetaLX = ::t(beta_linear);
+ for (int i=0; i<(ns*k); ++i){
+ beta_store(count,i) = tbeta[i];
+ beta_linear_store(count,i) = tbetaLX[i];
+ }
+ Matrix<> tgamma = ::t(gamma);
+ for (int i=0; i<(ns*gk); ++i)
+ gamma_store(count,i) = tgamma[i];
+ for (int j=0; j<ns*ns; ++j)
+ P_store(count,j)= P[j];
+ for (int l=0; l<N ; ++l)
+ ps_store(l,_) = ps_store(l,_) + ps(l,_);
+ s_store(count,_) = s;
+ Z_store(count,_) = Z;
+
+ ++count;
+
+ } // end of if(iter...)
+
+
+ // print output to stdout
+ if(iter > 1 && verbose > 0 && iter % verbose == 0){
+ Rprintf("\n\nMCMCoprobitChange iteration %i of %i \n", (iter+1), tot_iter);
+ for (int j = 0;j<ns; ++j){
+ double acceptancerate = accepts[j]/iter;
+ Rprintf("\n\n Acceptance rate for state %i is %10.5f \n", j+1, acceptancerate);
+ }
+ for (int j = 0;j<ns; ++j){
+ Rprintf("\n The number of observations in state %i is %10.5d", j+1, nstate[j]);
+ }
+ for (int j = 0; j<ns; ++j){
+ Rprintf("\n beta %i = ", j);
+ for (int i = 0; i<k; ++i){
+ Rprintf("%10.5f", beta(j, i));
+ }
+ }
+ // for (int j = 0; j<ns; ++j){
+ // Rprintf("\n beta_linear %i = ", j);
+ // for (int i = 0; i<k; ++i){
+ // Rprintf("%10.5f", beta_linear(j, i));
+ // }
+ // }
+ for (int j = 0; j<ns; ++j){
+ Rprintf("\n gamma %i = ", j);
+ for (int i = 0; i<(ncat-2); ++i){
+ Rprintf("%10.5f", gamma(j, i+2));
+ }
+ }
+ }
+
+ R_CheckUserInterrupt();
+
+ }// end MCMC loop
+
+ if(chib==1){
+ Matrix<> betast = meanc(beta_store);
+ Matrix<> betastLX = meanc(beta_linear_store);
+ Matrix<double, Row> beta_st(ns, k);
+ Matrix<double, Row> beta_linear_st(ns, k);
+ for (int j = 0; j<ns*k; ++j){
+ beta_st[j] = betast[j];
+ beta_linear_st[j] = betastLX[j];
+ }
+ Matrix<> gammast = meanc(gamma_store);
+ Matrix<double, Row> gamma_st(ns, gk);
+ for (int j = 0; j<ns*gk; ++j){
+ gamma_st[j] = gammast[j];
+ }
+
+ // for (int j = 0; j<ns; ++j){
+ // Rprintf("\n gamma_st %i = ", j);
+ // for (int i = 0; i<gk; ++i){
+ // Rprintf("%10.5f", gamma_st(j, i));
+ // }
+ // }
+
+ 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];
+ }
+ // storage
+ Matrix<> pdf_numer_store(nstore, 1);
+ Matrix<> pdf_alpha_store(nstore, 1);
+ Matrix<> pdf_P_store(nstore, ns);
+
+ // 1. gamma
+ Matrix<> densityq(nstore, ns);
+ Matrix<> alpha(nstore, ns);
+ for (int iter = 0; iter < nstore; ++iter){
+ int beta_count = 0;
+ Matrix<int> nstate(ns, 1);
+
+ Matrix<double, Row> gamma_g(ns, gk);
+ for (int h = 0; h<(ns*gk); ++h){
+ gamma_g[h] = gamma_store(iter, h);
+ }
+
+ Matrix<double, Row> beta_g(ns, k);
+ for (int h = 0; h<(ns*k); ++h){
+ beta_g[h] = beta_store(iter, h);
+ }
+ Matrix<> pdf_numer(ns, 1);
+ for (int j = 0; j <ns ; ++j){
+ for (int i = 0; i<N; ++i){
+ if (s_store(iter, i) == (j+1)) {
+ nstate[j] = nstate[j] + 1;
+ }
+ }
+ beta_count = beta_count + nstate[j];
+
+ Matrix<int> Yj = Y((beta_count - nstate[j]), 0, (beta_count - 1), 0);
+ Matrix<> Xj = X((beta_count - nstate[j]), 0, (beta_count - 1), k-1);
+ pdf_numer(j) = oprobit_log_postLX(j, ncat, gamma_st, gamma_g, Yj, Xj, beta_g,
+ tune, gammafixed);
+ for (int h = 2; h<ncat; ++h){
+ if (h == (ncat-1)){
+ densityq(iter, j) = densityq(iter, j) + log(dtnormLX(gamma_st(j, h), gamma_g(j, h), tune[j],
+ gamma_st(j, h-1), 300));
+ }
+ else {
+ densityq(iter, j) = densityq(iter, j) + log(dtnormLX(gamma_st(j, h), gamma_g(j, h), tune[j],
+ gamma_st(j, h-1), gamma_g(j, h+1)));
+ }
+ }
+ }
+
+ if (sum(pdf_numer) > 0){
+ pdf_numer_store(iter) = 0;
+ }
+ else{
+ pdf_numer_store(iter) = sum(pdf_numer);
+ }
+
+ }
+ double numerator = sum(meanc(pdf_numer_store)) + sum(meanc(densityq));
+
+ for (int iter = 0; iter < nstore; ++iter){
+ Matrix<> Sout = gaussian_ordinal_state_sampler_fixedsigma(stream, m, Y, X, beta_linear, Sigma, P);
+ Matrix<int> s = Sout(_, 0);
+ for (unsigned int i = 0; i<N; ++i) {
+ const Matrix<> mu = X(i,_)*t(beta(s[i]-1,_));
+ Z[i] = stream.rtnorm_combo(mu[0], 1.0, gamma_st(s[i]-1, Y[i]-1), gamma_st(s[i]-1, Y[i]));
+ }
+ int beta_count = 0;
+ Matrix<int> beta_count_storage(ns, 1);
+ Matrix<int> nstate(ns, 1);
+ for (int j = 0; j <ns ; ++j){
+ for (int i = 0; i<N; ++i){
+ if (s[i] == (j+1)) {
+ nstate[j] = nstate[j] + 1;
+ }
+ }
+ beta_count = beta_count + nstate[j];
+ Matrix<> yj = Y((beta_count - nstate[j]), 0, (beta_count - 1), 0);
+ Matrix<> Xj = X((beta_count - nstate[j]), 0, (beta_count - 1), k-1);
+ Matrix<> Zj = Z((beta_count - nstate[j]), 0, (beta_count - 1), 0);
+ Matrix<> XpX = t(Xj)*Xj;
+ Matrix<> XpZ = t(Xj)*Zj;
+ Matrix<> XpY = t(Xj)*yj;
+ Matrix<> Bn = invpd(B0inv + XpX/Sigma[0]);
+ Matrix<> bn = Bn*(B0inv*b0 + XpY/Sigma[0]);
+ beta_linear(j,_) = stream.rmvnorm(bn, Bn);
+ Matrix<> Bn2 = invpd(B0inv + XpX);
+ Matrix<> bn2 = Bn2*(B0inv*b0 + XpZ);
+ beta(j,_) = stream.rmvnorm(bn2, Bn2);
+ beta_count_storage[j] = beta_count;
+ }
+ // Sample P
+ double shape1 = 0;
+ double shape2 = 0;
+ P(ns-1, ns-1) = 1;
+
+ for (int j =0; j< (ns-1); ++j){
+ shape1 = A0(j,j) + nstate[j] - 1;
+ shape2 = A0(j,j+1) + 1; //
+ P(j,j) = stream.rbeta(shape1, shape2);
+ P(j,j+1) = 1 - P(j,j);
+ }
+ Matrix<> alpha(ns, 1);
+ for (int j = 0; j < ns ; ++j){
+ for (int i = 2; i< ncat; ++i){
+ if (i==(ncat-1)){
+ gamma_p(j, i) = stream.rtbnorm_combo(gamma_st(j, i), ::pow(tune[j], 2.0),
+ gamma_p(j, i-1));
+ }
+ else {
+ gamma_p(j, i) = stream.rtnorm_combo(gamma_st(j, i), ::pow(tune[j], 2.0),
+ gamma_p(j, i-1), gamma_st(j, i+1));
+ }
+ }
+ Matrix<int> Yj = Y((beta_count_storage[j] - nstate[j]), 0, (beta_count_storage[j] - 1), 0);
+ Matrix<> Xj = X((beta_count_storage[j] - nstate[j]), 0, (beta_count_storage[j] - 1), k-1);
+ alpha[j] = oprobit_log_postLX(j, ncat, gamma_p, gamma_st, Yj, Xj, beta, tune, gammafixed);
+ }
+ if (sum(alpha) > 0){
+ pdf_alpha_store(iter) = 0;
+ }
+ else{
+ pdf_alpha_store(iter) = sum(alpha);
+ }
+
+ }
+ double denominator = mean(pdf_alpha_store);
+ double pdf_gamma = numerator - denominator;
+
+ // 2. beta
+ Matrix<> density_beta(nstore, ns);
+ for (int iter = 0; iter < nstore; ++iter){
+ Matrix<> Sout = gaussian_ordinal_state_sampler_fixedsigma(stream, m, Y, X, beta_linear, Sigma, P);
+ Matrix<int> s = Sout(_, 0);
+ for (unsigned int i = 0; i<N; ++i) {
+ const Matrix<> mu = X(i,_)*t(beta(s[i]-1,_));
+ Z[i] = stream.rtnorm_combo(mu[0], 1.0, gamma_st(s[i]-1, Y[i]-1), gamma_st(s[i]-1, Y[i]));
+ }
+ int beta_count = 0;
+ Matrix<int> beta_count_storage(ns, 1);
+ Matrix<int> nstate(ns, 1);
+ for (int j = 0; j <ns ; ++j){
+ for (int i = 0; i<N; ++i){
+ if (s[i] == (j+1)) {
+ nstate[j] = nstate[j] + 1;
+ }
+ }
+ beta_count = beta_count + nstate[j];
+ Matrix<> yj = Y((beta_count - nstate[j]), 0, (beta_count - 1), 0);
+ Matrix<> Xj = X((beta_count - nstate[j]), 0, (beta_count - 1), k-1);
+ Matrix<> Zj = Z((beta_count - nstate[j]), 0, (beta_count - 1), 0);
+ Matrix<> XpX = t(Xj)*Xj;
+ Matrix<> XpZ = t(Xj)*Zj;
+ Matrix<> XpY = t(Xj)*yj;
+ Matrix<> Bn = invpd(B0inv + XpX/Sigma[0]);
+ Matrix<> bn = Bn*(B0inv*b0 + XpY/Sigma[0]);
+ beta_linear(j,_) = stream.rmvnorm(bn, Bn);
+ Matrix<> Bn2 = invpd(B0inv + XpX);
+ Matrix<> bn2 = Bn2*(B0inv*b0 + XpZ);
+ beta(j,_) = stream.rmvnorm(bn2, Bn2);
+ beta_count_storage[j] = beta_count;
+ density_beta(iter, j) = exp(lndmvn(t(beta_st(j,_)), bn2, Bn2));
+ }
+
+ // Sample P
+ double shape1 = 0;
+ double shape2 = 0;
+ P(ns-1, ns-1) = 1;
+
+ for (int j =0; j< (ns-1); ++j){
+ shape1 = A0(j,j) + nstate[j] - 1;
+ shape2 = A0(j,j+1) + 1; //
+ P(j,j) = stream.rbeta(shape1, shape2);
+ P(j,j+1) = 1 - P(j,j);
+ }
+ }
+
+ double pdf_beta = log(prod(meanc(density_beta)));
+
+ // 3. P
+ Matrix<> density_P(nstore, ns);
+ for (int iter = 0; iter < nstore; ++iter){
+ Matrix<> Sout = gaussian_ordinal_state_sampler_fixedsigma(stream, m, Y, X, beta_linear_st, Sigma, P);
+ Matrix <double> s = Sout(_, 0);
+ double shape1 = 0;
+ double shape2 = 0;
+ P(ns-1, ns-1) = 1;
+ Matrix <double> P_addN(ns, 1);
+ for (int j = 0; j<ns; ++j){
+ for (int i = 0; i<N; ++i){
+ if (s[i] == (j+1)) {
+ P_addN[j] = P_addN[j] + 1;
+ }
+ }
+ }
+
+ 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);
+ }
+ density_P(iter, ns-1) = 1;
+
+ }
+ double pdf_P = log(prod(meanc(density_P)));
+
+ // likelihood
+ Matrix<> F = Matrix<>(N, ns);
+ Matrix<> like(N, 1);
+ Matrix<> pr1 = Matrix<>(ns, 1);
+ pr1[0] = 1;
+ Matrix<> py(ns, 1);
+ Matrix<> pstyt1(ns, 1);
+
+ for (int t=0; t< N ; ++t){
+ Matrix<> mu = X(t,_)*::t(beta_st);
+ for (int j=0; j<ns ; ++j){
+ int yt = Y[t];
+ py[j] = oprobit_pdfLX(ncat, yt, mu[j], gamma_st(j,_));
+ }
+ if (t==0) pstyt1 = pr1;
+ else {
+ pstyt1 = ::t(F(t-1,_)*P_st);
+ }
+ Matrix<> unnorm_pstyt = pstyt1%py;
+ Matrix<> pstyt = unnorm_pstyt/sum(unnorm_pstyt);
+ for (int j=0; j<ns ; ++j){
+ F(t,j) = pstyt(j);
+ }
+ like[t] = sum(unnorm_pstyt);
+ }
+
+ loglike = sum(log(like));
+
+ // log prior ordinate
+ Matrix<> density_beta_prior(ns, 1);
+ Matrix<> density_P_prior(ns, 1);
+ density_P[ns-1] = 1; //
+
+ for (int j=0; j<ns ; ++j){
+ density_beta_prior[j] = lndmvn(::t(beta_st(j,_)), b0, B0);
+ }
+
+ 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)));
+ }
+ double density_gamma_prior = ns*(ncat-2)*log(dunif(1, 0, 10));
+
+ double logprior = sum(density_beta_prior) + sum(density_P_prior) + density_gamma_prior;
+ logmarglike = (loglike + logprior) - (pdf_beta + pdf_P + pdf_gamma);
+
+ /*
+ Rprintf("\n ----------------- marginal likelihood outputs ----------------- \n");
+ Rprintf("\n logmarglike %10.5f", logmarglike, "\n");
+ Rprintf("\n loglike %10.5f", loglike, "\n");
+ Rprintf("\n pdf_beta %10.5f", pdf_beta, "\n");
+ Rprintf("\n pdf_gamma %10.5f", pdf_gamma, "\n");
+ Rprintf("\n pdf_P %10.5f", pdf_P, "\n");
+ Rprintf("\n logprior %10.5f", logprior, "\n");
+ Rprintf("\n density_beta_prior %10.5f", sum(density_beta_prior), "\n");
+ Rprintf("\n density_gamma_prior %10.5f", density_gamma_prior, "\n");
+ Rprintf("\n density_P_prior %10.5f", sum(density_P_prior), "\n");
+ */
+ } // end of marginal likelihood
+
+}//end
+
+extern "C"{
+ void MCMCoprobitChange(double *betaout, double *betalinearout,
+ double *gammaout, double *Pout, double *psout, double *sout,
+ const double *Ydata,
+ const double *Xdata, const int *Xrow, const int *Xcol,
+ const int *m, const int *ncat,
+ const int *burnin, const int *mcmc, const int *thin, const int *verbose,
+ const double *tunedata, // const int *tdfdata,
+ const int *uselecuyer, const int *seedarray, const int *lecuyerstream,
+ const double *betastart, const double *betalinearstart,
+ const double *gammastart, const double *Pstart,
+ const double *sigmastart, const double *a, const double *b,
+ const double *b0data, const double *B0data,
+ const double *A0data,
+ double *logmarglikeholder, double *loglikeholder,
+ const int *chib, const int *gammafixed)
+ {
+
+ // pull together Matrix objects
+ const Matrix<> Y(*Xrow, 1, Ydata);
+ const Matrix<> X(*Xrow, *Xcol, Xdata);
+ const unsigned int nstore = *mcmc / *thin;
+ const unsigned int N = *Xrow;
+ const unsigned int k = *Xcol;
+ const unsigned int gk = *ncat + 1;
+ const unsigned int ns = *m + 1;
+
+ // generate starting values
+ Matrix<> beta(ns, k, betastart);
+ Matrix<> beta_linear(ns, k, betalinearstart);
+ Matrix<> Sigma(1, 1, sigmastart);
+ Matrix<> P(ns, ns, Pstart);
+ Matrix<> b0(k, 1, b0data);
+ Matrix<> B0(k, k, B0data);
+ Matrix<> tune(ns, 1, tunedata);
+ Matrix<> A0(ns, ns, A0data);
+ double logmarglike;
+ double loglike;
+
+ // storage matrices
+ Matrix<> beta_store(nstore, ns*k);
+ Matrix<> beta_linear_store(nstore, ns*k);
+ Matrix<> Z_store(nstore, N);
+ Matrix<> P_store(nstore, ns*ns);
+ Matrix<> ps_store(N, ns);
+ Matrix<int> s_store(nstore, N);
+
+ Matrix<> gamma(ns, gk, gammastart);
+ Matrix<> gamma_store(nstore, ns*gk);
+ MCMCPACK_PASSRNG2MODEL(MCMCoprobitChange_impl, *m, *ncat,
+ Y, X, beta, beta_linear, gamma, P, Sigma,
+ b0, B0, A0,
+ *burnin, *mcmc, *thin, *verbose, tune,
+ *chib, *gammafixed,
+ beta_store, beta_linear_store, gamma_store, Z_store,
+ P_store, ps_store, s_store,
+ logmarglike, loglike);
+
+ logmarglikeholder[0] = logmarglike;
+ loglikeholder[0] = loglike;
+
+ // return output
+ for (int i = 0; i<(nstore*ns*k); ++i){
+ betaout[i] = beta_store[i];
+ betalinearout[i] = beta_linear_store[i];
+ }
+ for (int i = 0; i<(nstore*ns*gk); ++i){
+ gammaout[i] = gamma_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];
+ }
+ }
+}
+#endif
+
+
diff --git a/src/MCMCpoissonChange.cc b/src/MCMCpoissonChange.cc
index 7c071b5..e9af4c5 100644
--- a/src/MCMCpoissonChange.cc
+++ b/src/MCMCpoissonChange.cc
@@ -298,6 +298,7 @@ void MCMCpoissonChangepoint_impl(rng<RNGTYPE>& stream,
const double *b,
const double *A0data,
double *logmarglikeholder,
+ double *loglikeholder,
const int *chib)
{
@@ -506,6 +507,7 @@ void MCMCpoissonChangepoint_impl(rng<RNGTYPE>& stream,
const double logmarglike = (loglike + logprior) - (pdf_lambda + pdf_P);
logmarglikeholder[0] = logmarglike;
+ loglikeholder[0] = loglike;
}
R_CheckUserInterrupt();
@@ -556,6 +558,7 @@ void MCMCpoissonRegChangepoint_impl(rng<RNGTYPE>& stream,
const double *B0data,
const double *A0data,
double *logmarglikeholder,
+ double *loglikeholder,
const double* wrin,
const double* mrin,
const double* srin,
@@ -850,7 +853,8 @@ void MCMCpoissonRegChangepoint_impl(rng<RNGTYPE>& stream,
const double logmarglike = (loglike + logprior) - (pdf_beta + pdf_P);
logmarglikeholder[0] = logmarglike;
-
+ loglikeholder[0] = loglike;
+
}
R_CheckUserInterrupt();
@@ -901,6 +905,7 @@ extern "C" {
const double *B0data,
const double *A0data,
double *logmarglikeholder,
+ double *loglikeholder,
const double *wrin,
const double *mrin,
const double *srin,
@@ -913,7 +918,7 @@ extern "C" {
burnin, mcmc, thin, verbose,
betastart, Pstart,
a, b, A0data,
- logmarglikeholder, chib)
+ logmarglikeholder, loglikeholder, chib)
}
else{
MCMCPACK_PASSRNG2MODEL(MCMCpoissonRegChangepoint_impl,
@@ -924,7 +929,7 @@ extern "C" {
betastart, Pstart,
taustart, componentstart,
a, b, b0data, B0data, A0data,
- logmarglikeholder,
+ logmarglikeholder, loglikeholder,
wrin, mrin, srin, chib);
}
} // end MCMC
diff --git a/src/MCMCprobit.cc b/src/MCMCprobit.cc
index d4c954c..fb9b2cf 100644
--- a/src/MCMCprobit.cc
+++ b/src/MCMCprobit.cc
@@ -50,59 +50,97 @@ using namespace scythe;
*/
template <typename RNGTYPE>
void MCMCprobit_impl (rng<RNGTYPE>& stream, const Matrix<>& Y,
- const Matrix<>& X, Matrix<>& beta, const Matrix<>& b0,
- const Matrix<>& B0, unsigned int burnin, unsigned int mcmc,
- unsigned int thin, unsigned int verbose, Matrix<>& result) {
-
- // define constants and from cross-product matrices
- const unsigned int tot_iter = burnin + mcmc; // total iterations
- const unsigned int nstore = mcmc / thin; // number of draws to store
- const unsigned int k = X.cols();
- const unsigned int N = X.rows();
- const Matrix<> XpX = crossprod(X);
+ const Matrix<>& X, Matrix<>& beta, const Matrix<>& b0,
+ const Matrix<>& B0, unsigned int burnin, unsigned int mcmc,
+ unsigned int thin, unsigned int verbose, bool chib,
+ Matrix<>& result,
+ double& logmarglike) {
+
+ // define constants and from cross-product matrices
+ const unsigned int tot_iter = burnin + mcmc; // total iterations
+ const unsigned int nstore = mcmc / thin; // number of draws to store
+ const unsigned int k = X.cols();
+ const unsigned int N = X.rows();
+ const Matrix<> XpX = crossprod(X);
+ const Matrix<> B0inv = invpd(B0);
+
+ // storage matrix or matrices
+ Matrix<> beta_store(nstore, k);
+ Matrix<> Z_store(nstore, N);
- // storage matrix or matrices
- Matrix<> storemat(nstore, k);
+ // initialize Z
+ Matrix<> Z(N,1);
+
+ // MCMC sampling starts here
+ unsigned int count = 0;
+ for (unsigned int iter = 0; iter < tot_iter; ++iter) {
- // initialize Z
- Matrix<> Z(N,1);
-
- // MCMC sampling starts here
- unsigned int count = 0;
- for (unsigned int iter = 0; iter < tot_iter; ++iter) {
-
- // [Z| beta, y]
- const Matrix<> Z_mean = X * beta;
- for (unsigned int i=0; i<N; ++i){
- if (Y[i] == 1.0)
- Z[i] = stream.rtbnorm_combo(Z_mean[i], 1.0, 0);
- if (Y[i] == 0.0)
- Z[i] = stream.rtanorm_combo(Z_mean[i], 1.0, 0);
- }
-
- // [beta|Z]
- 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)){
- storemat(count,_) = beta;
- ++count;
- }
-
- // print output to stdout
- if(verbose > 0 && iter % verbose == 0){
- Rprintf("\n\nMCMCprobit 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]);
- }
-
- R_CheckUserInterrupt(); // allow user interrupts
+ // [Z| beta, y]
+ const Matrix<> Z_mean = X * beta;
+ for (unsigned int i=0; i<N; ++i){
+ if (Y[i] == 1.0)
+ Z[i] = stream.rtbnorm_combo(Z_mean[i], 1.0, 0);
+ if (Y[i] == 0.0)
+ Z[i] = stream.rtanorm_combo(Z_mean[i], 1.0, 0);
+ }
+
+ // [beta|Z]
+ 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)){
+ beta_store(count,_) = beta;
+ Z_store(count,_) = Z;
+ ++count;
+ }
+
+ // print output to stdout
+ if(verbose > 0 && iter % verbose == 0){
+ Rprintf("\n\nMCMCprobit 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]);
+ }
+
+ R_CheckUserInterrupt(); // allow user interrupts
+ } // end of MCMC loop
+
+
+ if(chib==1){
+ // Rprintf("\n Marginal Likelihood Computation Starts!\n");
+ Matrix<double> beta_star = meanc(beta_store);
+ const Matrix<double> Z_reduced(N, 1);
+ double logbeta_sum = 0;
+ for (int iter = 0; iter<nstore; ++iter){
+ Z_reduced(_,0) = Z_store(iter,_);
+ const Matrix<double> XpZ = (::t(X)*Z_reduced);
+ const Matrix<double> Bn = invpd(B0inv + XpX);
+ const Matrix<double> bn = Bn*gaxpy(B0inv, b0, XpZ);
+ logbeta_sum += lndmvn(beta_star, bn, Bn);
}
- result = storemat;
+ double logbeta = logbeta_sum/nstore;
+
+ double loglike = 0.0;
+ Matrix<> eta = X * beta_star;
+ for (unsigned int i = 0; i < N; ++i) {
+ double phi = pnorm(eta(i), 0, 1);
+ loglike += log(dbinom(Y(i), 1, phi));
+ }
+
+ // calculate log prior ordinate
+ double logprior = lndmvn(beta_star, b0, B0inv);
+
+ logmarglike = loglike + logprior - logbeta;
+
+ // Rprintf("\n logmarglike %10.5f", logmarglike, "\n");
+ // Rprintf("\n loglike %10.5f", loglike, "\n");
+
+ }// end of marginal likelihood computation
+
+ result = beta_store;
}
-
+
extern "C"{
void MCMCprobit(double *sampledata, const int *samplerow,
@@ -114,7 +152,9 @@ extern "C"{
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 *B0data, const int *B0row, const int *B0col,
+ double *logmarglikeholder, // double *loglikeholder,
+ const int *chib) {
// pull together Matrix objects
const Matrix <> Y(*Yrow, *Ycol, Ydata);
@@ -123,11 +163,17 @@ extern "C"{
betastartdata);
const Matrix <> b0(*b0row, *b0col, b0data);
const Matrix <> B0(*B0row, *B0col, B0data);
+ double logmarglike;
+ // double loglike;
Matrix<> storagematrix;
MCMCPACK_PASSRNG2MODEL(MCMCprobit_impl, Y, X, beta, b0, B0, *burnin,
- *mcmc, *thin, *verbose, storagematrix);
-
+ *mcmc, *thin, *verbose, *chib,
+ storagematrix,
+ logmarglike);
+ logmarglikeholder[0] = logmarglike;
+ // loglikeholder[0] = loglike;
+
const unsigned int size = *samplerow * *samplecol;
for (unsigned int i=0; i<size; ++i)
sampledata[i] = storagematrix(i);
diff --git a/src/MCMCprobitChange.cc b/src/MCMCprobitChange.cc
new file mode 100644
index 0000000..ced6161
--- /dev/null
+++ b/src/MCMCprobitChange.cc
@@ -0,0 +1,415 @@
+// MCMCprobitChange.cc is C++ code to estimate a probit changepoint model
+// with a beta prior
+//
+// Jong Hee Park
+// Dept. of Political Science
+// University of Chicago
+// jhp at uchicago.edu
+//
+// Written by JHP 07/06/2007
+// Modifieed by JHP 11/02/2009
+// Included by JHP 1/21/2011
+// Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+// and Jong Hee Park
+//////////////////////////////////////////////////////////////////////////
+
+#ifndef MCMCPROBITCHANGE_CC
+#define MCMCPROBITCHANGE_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>
+#include <R_ext/Utils.h>
+
+
+using namespace std;
+using namespace scythe;
+
+// probit state sampler
+template <typename RNGTYPE>
+
+Matrix<> probit_state_sampler(rng<RNGTYPE>& stream,
+ const int m,
+ const Matrix<double>& Y,
+ const Matrix<double>& X,
+ const Matrix<double>& beta,
+ const Matrix<double>& P){
+
+ const int ns = m + 1;
+ const int n = Y.rows();
+ Matrix<double> F = Matrix<double>(n, ns);
+ Matrix<double> pr1 = Matrix<double>(ns, 1);
+ pr1[0] = 1;
+ Matrix<double> py(ns, 1);
+ Matrix<double> pstyt1(ns, 1);
+
+ for (int t=0; t<n ; ++t){
+ int yt = (int) Y[t];
+ Matrix<double> mu = X(t,_)*::t(beta);
+ for (int j=0; j<ns ; ++j){
+ py[j] = dbinom(yt, 1, pnorm(mu[j], 0, 1));
+ }
+ if (t==0) pstyt1 = pr1;
+ else {
+ pstyt1 = ::t(F(t-1,_)*P);
+ }
+ Matrix<double> unnorm_pstyt = pstyt1%py;
+ const Matrix<double> pstyt = unnorm_pstyt/sum(unnorm_pstyt);
+ for (int j=0; j<ns ; ++j){
+ F(t,j) = pstyt(j);
+ }
+ }
+
+ Matrix<int> s(n, 1);
+ Matrix<double> ps = Matrix<double>(n, ns);
+ ps(n-1,_) = F(n-1,_);
+ s(n-1) = ns;
+
+ Matrix<double> pstyn = Matrix<double>(ns, 1);
+ double pone = 0.0;
+ int t = n-2;
+ while (t >= 0){
+ int st = s(t+1);
+ Matrix<double> Pst_1 = ::t(P(_,st-1));
+ Matrix<double> 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;
+ }
+
+ Matrix<double> Sout(n, ns+1);
+ Sout(_, 0) = s(_,0);
+ for (int j = 0; j<ns; ++j){
+ Sout(_,j+1) = ps(_, j);
+ }
+
+ return Sout;
+
+}
+
+////////////////////////////////////////////
+// MCMCprobitChangepoint implementation.
+////////////////////////////////////////////
+template <typename RNGTYPE>
+void MCMCprobitChange_impl(rng<RNGTYPE>& stream,
+ const int m,
+ const Matrix<>& Y, const Matrix<>& X,
+ Matrix<>& beta, Matrix<>& P,
+ Matrix<>& b0, Matrix<>& B0,
+ const Matrix<>& A0,
+ unsigned int burnin, unsigned int mcmc, unsigned int thin,
+ unsigned int verbose, bool chib,
+ Matrix<>& beta_store, Matrix<>& Z_store,
+ Matrix<>& P_store, Matrix<>& ps_store, Matrix<int>& s_store,
+ double& logmarglike, double& loglike)
+{
+ const int tot_iter = burnin + mcmc;
+ const int nstore = mcmc / thin;
+ const int n = Y.rows();
+ const int ns = m + 1;
+ const int k = X.cols();
+ const Matrix<> B0inv = invpd(B0);
+ Matrix<> Z(n, 1);
+
+ unsigned int count = 0;
+ for (int iter = 0; iter < tot_iter; ++iter){
+
+ // 1. Sample s
+ Matrix<> Sout = probit_state_sampler(stream, m, Y, X, beta, P);
+ Matrix<int> s = Sout(_, 0);
+ Matrix <double> ps(n, ns);
+ for (int j = 0; j<ns; ++j){
+ ps(_,j) = Sout(_,j+1);
+ }
+
+ // 2. Sample Z
+ for (unsigned int i = 0; i < n; ++i) {
+ const Matrix<> mu = X(i,_)*t(beta(s[i]-1,_));
+ double muj = mu[0];
+ if(muj>200){
+ muj = 200;
+ }
+ if(muj<-200){
+ muj = -200;
+ }
+ if (Y[i] == 1.0)
+ Z[i] = stream.rtbnorm_combo(muj, 1.0, 0);
+ if (Y[i] == 0.0)
+ Z[i] = stream.rtanorm_combo(muj, 1.0, 0);
+ }
+
+ // 3. Sample beta
+ int beta_count = 0;
+ Matrix<int> nstate(ns, 1);
+ for (int j = 0; j <ns ; ++j){
+ for (int i = 0; i<n; ++i){
+ if (s[i] == (j+1)) {
+ nstate[j] = nstate[j] + 1;
+ }
+ }
+ beta_count = beta_count + nstate[j];
+ const Matrix<double> Zj = Z((beta_count - nstate[j]), 0, (beta_count - 1), 0);
+ const Matrix<double> Xj = X((beta_count - nstate[j]), 0, (beta_count - 1), k-1);
+ const Matrix<double> XpX = t(Xj)*Xj;
+ const Matrix<double> XpZ = t(Xj)*Zj;
+ beta(j,_) = NormNormregress_beta_draw(XpX, XpZ, b0, B0, 1.0, stream);
+ }
+
+ // 4. Sample P
+ double shape1 = 0;
+ double shape2 = 0;
+ P(ns-1, ns-1) = 1;
+
+ for (int j =0; j< (ns-1); ++j){
+ shape1 = A0(j,j) + nstate[j] - 1;
+ shape2 = A0(j,j+1) + 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)){
+ Matrix<double> tbeta = ::t(beta);
+ for (int i=0; i<(ns*k); ++i)
+ beta_store(count,i) = tbeta[i];
+ for (int j=0; j<ns*ns; ++j)
+ P_store(count,j)= P[j];
+ for (int l=0; l<n ; ++l)
+ ps_store(l,_) = ps_store(l,_) + ps(l,_);
+ s_store(count,_) = s;
+ Z_store(count,_) = Z;
+
+ ++count;
+
+ } // end of if(iter...)
+
+
+ // print output to stdout
+ if(verbose > 0 && iter % verbose == 0){
+ Rprintf("\n\nMCMCprobitChange iteration %i of %i \n", (iter+1), tot_iter);
+ for (int j = 0;j<ns; ++j){
+ Rprintf("\n The number of observations in state %i is %10.5d", j+1, nstate[j]);
+ }
+ for (int j = 0; j<ns; ++j){
+ Rprintf("\n beta %i is ", j);
+ for (int i = 0; i<k; ++i){
+ Rprintf("%10.5f", beta(j, i));
+ }
+ }
+ }
+
+ R_CheckUserInterrupt();
+
+ }// end MCMC loop
+
+ if(chib==1){
+ Matrix<double> betast = meanc(beta_store);
+ Matrix<double, Row> beta_st(ns, k);
+ for (int j = 0; j<ns*k; ++j){
+ beta_st[j] = betast[j];
+ }
+ Matrix<double> P_vec_st = meanc(P_store);
+ const Matrix<double> P_st(ns, ns);
+ for (int j = 0; j< ns*ns; ++j){
+ P_st[j] = P_vec_st[j];
+ }
+
+ // 1. beta
+ Matrix<double> density_beta(nstore, ns);
+ for (int iter = 0; iter<nstore; ++iter){
+ Matrix<int> nstate(ns, 1);
+ int beta_count = 0;
+ const Matrix<double> Z(n, 1);
+ Z(_,0) = Z_store(iter,_);
+ for (int j = 0; j<ns ; ++j){
+ for (int i = 0; i<n; ++i){
+ if (s_store(iter, i) == (j+1)) {
+ nstate[j] = nstate[j] + 1;
+ }
+ }
+ beta_count = beta_count + nstate[j];
+ const Matrix<double> Zj = Z((beta_count - nstate[j]), 0, (beta_count - 1), 0);
+ const Matrix<double> Xj = X((beta_count - nstate[j]), 0, (beta_count - 1), k-1);
+ const Matrix<double> XpX = (::t(Xj)*Xj);
+ const Matrix<double> XpZ = (::t(Xj)*Zj);
+ const Matrix<double> Bn = invpd(B0inv + XpX);
+ const Matrix<double> bn = Bn*gaxpy(B0inv, b0, XpZ);
+ density_beta(iter, j) = exp(lndmvn(::t(beta_st(j,_)), bn, Bn));
+ }
+ }
+
+ double pdf_beta = log(prod(meanc(density_beta)));
+
+ // 2. P
+ Matrix<double> density_P(nstore, ns);
+ for (int iter = 0; iter < nstore; ++iter){
+ Matrix <double> Sout = probit_state_sampler(stream, m, Y, X, beta_st, P);
+ Matrix <double> s = Sout(_, 0);
+ Matrix <double> 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 <double> P_addN(ns, 1);
+ for (int j = 0; j<ns; ++j){
+ for (int i = 0; i<n; ++i){
+ if (s[i] == (j+1)) {
+ P_addN[j] = P_addN[j] + 1;
+ }
+ }
+ }
+ 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);
+ }
+ density_P(iter, ns-1) = 1;
+
+ }
+ double pdf_P = log(prod(meanc(density_P)));
+
+ // likelihood
+ Matrix<double> F = Matrix<double>(n, ns);
+ Matrix<double> like(n, 1);
+ Matrix<double> pr1 = Matrix<double>(ns, 1);
+ pr1[0] = 1;
+ Matrix<double> py(ns, 1);
+ Matrix<double> pstyt1(ns, 1);
+
+ for (int t=0; t<n ; ++t){
+ int yt = (int) Y[t];
+ Matrix<double> mu = X(t,_)*::t(beta_st);
+ for (int j=0; j<ns ; ++j){
+ py[j] = dbinom(yt, 1, pnorm(mu[j], 0, 1));
+ }
+ if (t==0) pstyt1 = pr1;
+ else {
+ pstyt1 = ::t(F(t-1,_)*P_st);
+ }
+ Matrix<double> unnorm_pstyt = pstyt1%py;
+ Matrix<double> pstyt = unnorm_pstyt/sum(unnorm_pstyt);
+ for (int j=0; j<ns ; ++j){
+ F(t,j) = pstyt(j);
+ }
+ like[t] = sum(unnorm_pstyt);
+ }
+
+ loglike = sum(log(like));
+
+ //////////////////////
+ // log prior ordinate
+ //////////////////////
+ Matrix<double> density_beta_prior(ns, 1);
+ Matrix<double> density_P_prior(ns, 1);
+ density_P[ns-1] = 1; //
+
+ for (int j=0; j<ns ; ++j){
+ density_beta_prior[j] = lndmvn(::t(beta_st(j,_)), b0, B0);
+ }
+
+ 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_beta_prior) + sum(density_P_prior);
+ logmarglike = (loglike + logprior) - (pdf_beta + pdf_P);
+
+ // Rprintf("\n density_beta_prior %10.5f", density_beta_prior[0], "\n");
+ // Rprintf("\n density_P_prior %10.5f", density_P_prior[0], "\n");
+ Rprintf("\n logmarglike %10.5f", logmarglike, "\n");
+ Rprintf("\n loglike %10.5f", loglike, "\n");
+ // Rprintf("\n logprior %10.5f", logprior, "\n");
+ // Rprintf("\n pdf_beta %10.5f", pdf_beta, "\n");
+ // Rprintf("\n pdf_P %10.5f", pdf_P, "\n");
+ } // end of marginal likelihood
+}//end extern "C"
+
+extern "C"{
+ void MCMCprobitChange(double *betaout, double *Pout, double *psout, double *sout,
+ const double *Ydata, const int *Yrow, const int *Ycol,
+ const double *Xdata, const int *Xrow, const int *Xcol,
+ const int *m,
+ const int *burnin, const int *mcmc, const int *thin, const int *verbose,
+ const int *uselecuyer, const int *seedarray, const int *lecuyerstream,
+ const double *betastart, const double *Pstart,
+ const double *a, const double *b,
+ const double *b0data, const double *B0data,
+ const double *A0data,
+ double *logmarglikeholder, double *loglikeholder,
+ const int *chib){
+
+ // pull together Matrix objects
+ const Matrix <> Y(*Yrow, *Ycol, Ydata);
+ const Matrix <> X(*Xrow, *Xcol, Xdata);
+ const unsigned int tot_iter = *burnin + *mcmc;
+ const unsigned int nstore = *mcmc / *thin;
+ const int n = Y.rows();
+ const int k = X.cols();
+ const int ns = *m + 1;
+
+ // generate starting values
+ Matrix <> beta(ns, k, betastart);
+ Matrix <> P(ns, ns, Pstart);
+ Matrix <> b0(k, 1, b0data);
+ Matrix <> B0(k, k, B0data);
+ const Matrix <> A0(ns, ns, A0data);
+ double logmarglike;
+ double loglike;
+
+ // storage matrices
+ Matrix<> beta_store(nstore, ns*k);
+ Matrix<> Z_store(nstore, n);
+ Matrix<> P_store(nstore, ns*ns);
+ Matrix<> ps_store(n, ns);
+ Matrix<int> s_store(nstore, n);
+
+ MCMCPACK_PASSRNG2MODEL(MCMCprobitChange_impl, *m,
+ Y, X, beta, P, b0, B0, A0,
+ *burnin, *mcmc, *thin, *verbose, *chib,
+ beta_store, Z_store,
+ P_store, ps_store, s_store,
+ logmarglike, loglike);
+ logmarglikeholder[0] = logmarglike;
+ loglikeholder[0] = loglike;
+
+ // return output
+ for (int i = 0; i<(nstore*ns*k); ++i){
+ betaout[i] = beta_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];
+ }
+
+ }
+}
+#endif
+
+
diff --git a/src/MCMCprobitres.cc b/src/MCMCprobitres.cc
index e893c83..ea13701 100644
--- a/src/MCMCprobitres.cc
+++ b/src/MCMCprobitres.cc
@@ -56,8 +56,9 @@ void MCMCprobitres_impl (rng<RNGTYPE>& stream, const Matrix<>& Y,
const Matrix<>& b0,
const Matrix<>& B0, unsigned int burnin,
unsigned int mcmc,
- unsigned int thin, unsigned int verbose,
- Matrix<>& result) {
+ unsigned int thin, unsigned int verbose, bool chib,
+ Matrix<>& result,
+ double& logmarglike) {
// define constants and from cross-product matrices
const unsigned int tot_iter = burnin + mcmc; // total number of mcmc iterations
@@ -65,10 +66,12 @@ void MCMCprobitres_impl (rng<RNGTYPE>& stream, const Matrix<>& Y,
const unsigned int k = X.cols();
const unsigned int N = X.rows();
const Matrix<> XpX = crossprod(X);
+ const Matrix<> B0inv = invpd(B0);
- // holding matrices
- Matrix<> storemat(nstore, k+resvec.rows());
-
+ // storage matrix or matrices
+ Matrix<> beta_store(nstore, k);
+ Matrix<> Z_store(nstore, N);
+
// initialize Z
Matrix<> Z(N,1);
@@ -91,11 +94,13 @@ void MCMCprobitres_impl (rng<RNGTYPE>& stream, const Matrix<>& Y,
// store values in matrices
if (iter >= burnin && ((iter % thin)==0)){
- for (unsigned int j = 0; j < k; j++)
- storemat(count, j) = beta[j];
+ for (unsigned int j = 0; j < k; j++){
+ beta_store(count, j) = beta[j];
+ }
+ Z_store(count,_) = Z;
for (unsigned int j=0; j<(resvec.rows()); ++j){
const int i = static_cast<int>(resvec[j]) - 1;
- storemat(count, j+k) = Z[i] - Z_mean[i];
+ beta_store(count, j+k) = Z[i] - Z_mean[i];
}
++count;
}
@@ -112,8 +117,39 @@ void MCMCprobitres_impl (rng<RNGTYPE>& stream, const Matrix<>& Y,
} // end MCMC loop
- result = storemat;
+ if(chib==1){
+ Rprintf("\n Marginal Likelihood Computation Starts!\n");
+
+ Matrix<double> beta_star = meanc(beta_store);
+ Matrix<double> density_beta(nstore, 1);
+ for (int iter = 0; iter<nstore; ++iter){
+ const Matrix<> Z_reduced = Z_store(iter,_);
+ const Matrix<double> XpZ = (::t(X)*Z_reduced);
+ const Matrix<double> Bn = invpd(B0inv + XpX);
+ const Matrix<double> bn = Bn*gaxpy(B0inv, b0, XpZ);
+ density_beta(iter) = exp(lndmvn(beta_star, bn, Bn));
+ }
+ double logbeta = log(prod(meanc(density_beta)));
+
+ double loglike = 0.0;
+ Matrix<> eta = X * beta_star;
+ for (unsigned int i = 0; i < X.rows(); ++i) {
+ double phi = pnorm(eta(i), 0, 1);
+ loglike += log(dbinom(Y(i), 1, phi));
+ }
+
+ // calculate log prior ordinate
+ double logprior = lndmvn(beta_star, b0, B0inv);
+
+ logmarglike = loglike + logprior - logbeta;
+
+ Rprintf("\n logmarglike %10.5f", logmarglike, "\n");
+ Rprintf("\n loglike %10.5f", loglike, "\n");
+
+ }// end of marginal likelihood computation
+
+ result = beta_store;
}
@@ -133,7 +169,10 @@ extern "C"{
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 int *B0row, const int *B0col,
+ double *logmarglikeholder, // double *loglikeholder,
+ const int *chib) {
+
// pull together Matrix objects
const Matrix <> Y(*Yrow, *Ycol, Ydata);
@@ -144,12 +183,14 @@ extern "C"{
betastartdata);
const Matrix <> b0(*b0row, *b0col, b0data);
const Matrix <> B0(*B0row, *B0col, B0data);
-
+ double logmarglike;
+
Matrix<> storagematrix;
MCMCPACK_PASSRNG2MODEL(MCMCprobitres_impl, Y, X, beta, resvec,
b0, B0, *burnin,
- *mcmc, *thin, *verbose,
- storagematrix);
+ *mcmc, *thin, *verbose, *chib,
+ storagematrix,
+ logmarglike);
// return output
const unsigned int size = *samplerow * *samplecol;
--
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