[r-cran-mcmcpack] 36/90: Imported Upstream version 1.0-2
Andreas Tille
tille at debian.org
Fri Dec 16 09:07:38 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 ef60022125cafc6dfa5a99ef4547b7e2fde9478c
Author: Andreas Tille <tille at debian.org>
Date: Fri Dec 16 08:07:14 2016 +0100
Imported Upstream version 1.0-2
---
DESCRIPTION | 8 +-
HISTORY | 10 +
NAMESPACE | 3 +-
R/MCMCbinaryChange.R | 100 +++
R/MCMCpoissonChange.R | 139 +++
R/MCMCpoissonChangepoint.R | 112 ---
R/btsutil.R | 300 ++++---
...MCpoissonChangepoint.Rd => MCMCbinaryChange.Rd} | 255 +++---
man/MCMCpoissonChange.Rd | 167 ++++
man/plotChangepoint.Rd | 14 +-
man/plotState.Rd | 11 +-
src/MCMCbinaryChange.cc | 363 ++++++++
src/MCMCpoissonChange.cc | 937 +++++++++++++++++++++
src/MCMCpoissonChangepoint.cc | 386 ---------
14 files changed, 2065 insertions(+), 740 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index 4190fa1..63ef4be 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
Package: MCMCpack
-Version: 1.0-1
-Date: 2009-6-29
+Version: 1.0-2
+Date: 2009-7-17
Title: Markov chain Monte Carlo (MCMC) Package
Author: Andrew D. Martin <admartin at wustl.edu>, Kevin M. Quinn
<kevin_quinn at harvard.edu>, Jong Hee Park <jhp at uchicago.edu>
@@ -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: 2009-06-30 09:23:24 UTC; adm
+Packaged: 2009-07-17 18:52:45 UTC; adm
Repository: CRAN
-Date/Publication: 2009-07-01 12:02:06
+Date/Publication: 2009-07-18 15:22:20
diff --git a/HISTORY b/HISTORY
index 94a4ebf..698091d 100644
--- a/HISTORY
+++ b/HISTORY
@@ -2,6 +2,16 @@
// Changes and Bug Fixes
//
+1.0-1 to 1.0-2
+ * Added Poisson regression changepoint analysis MCMCpoissonChange() [JHP]
+ * Old MCMCpoissonChangepoint() is deprecated [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]
+
+
+
0.9-6 to 1.0-1
* added one-dimensional dynamic IRT model
* added Rehnquist Court data to illustrate dynamic IRT model
diff --git a/NAMESPACE b/NAMESPACE
index eade73c..cda65b0 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -16,6 +16,7 @@ export(
MCpoissongamma,
MCnormalnormal,
MCmultinomdirichlet,
+ MCMCbinaryChange,
MCMCdynamicEI,
MCMCdynamicIRT1d,
MCMCfactanal,
@@ -31,7 +32,7 @@ export(
MCMCoprobit,
MCMCordfactanal,
MCMCpoisson,
- MCMCpoissonChangepoint,
+ MCMCpoissonChange,
MCMCprobit,
MCMCregress,
MCMCSVDreg,
diff --git a/R/MCMCbinaryChange.R b/R/MCMCbinaryChange.R
new file mode 100644
index 0000000..b720031
--- /dev/null
+++ b/R/MCMCbinaryChange.R
@@ -0,0 +1,100 @@
+## sample from the posterior distribution
+## of a binary model with multiple changepoints
+## using linked C++ code in Scythe
+##
+## Written by JHP 07/01/2007
+## Revised by JHP 07/16/2009
+
+"MCMCbinaryChange"<-
+ function(data, m = 1, c0 = 1, d0 = 1, a = NULL, b = NULL,
+ burnin = 10000, mcmc = 10000, thin = 1, verbose = 0,
+ seed = NA, phi.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()
+
+ ## seeds
+ seeds <- form.seeds(seed)
+ lecuyer <- seeds[[1]]
+ seed.array <- seeds[[2]]
+ lecuyer.stream <- seeds[[3]]
+ if(!is.na(seed)) set.seed(seed)
+
+ ## sample size
+ y <- as.matrix(data)
+ n <- nrow(y)
+ ns <- m+1
+
+ ## prior for transition matrix
+ A0 <- trans.mat.prior(m=m, n=n, a=a, b=b)
+
+ ## 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
+ }
+
+ ## starting values
+ Pstart <- check.P(P.start, m=m, a=a, b=b)
+ phistart <- check.theta(phi.start, ns, y, min=0, max=1)
+
+ nstore <- mcmc/thin
+
+ ## call C++ code to draw sample
+ posterior <- .C("MCMCbinaryChange",
+ phiout = 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),
+ verbose = as.integer(verbose),
+ lecuyer=as.integer(lecuyer),
+ seedarray=as.integer(seed.array),
+ lecuyerstream=as.integer(lecuyer.stream),
+ phistart = as.double(phistart),
+ 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
+ }
+
+ ## pull together matrix and build MCMC object to return
+ phi.holder <- matrix(posterior$phiout, nstore, )
+ P.holder <- matrix(posterior$Pout, nstore, )
+ s.holder <- matrix(posterior$sout, nstore, )
+ ps.holder <- matrix(posterior$psout, n, )
+
+ output <- mcmc(data=phi.holder, start=burnin+1, end=burnin + mcmc, thin=thin)
+ varnames(output) <- paste("phi.", 1:ns, sep = "")
+ attr(output,"title") <- "MCMCbinaryChange Posterior Sample"
+ attr(output, "y") <- y
+ attr(output, "m") <- m
+ attr(output, "call") <- cl
+ attr(output, "logmarglike") <- logmarglike
+ 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
new file mode 100644
index 0000000..4dc90e5
--- /dev/null
+++ b/R/MCMCpoissonChange.R
@@ -0,0 +1,139 @@
+################################
+## Poisson Changepoint Model
+##
+## 07/14/2009 Jong Hee Park
+################################
+
+"MCMCpoissonChange"<-
+ function(formula, data = parent.frame(), m = 1,
+ b0 = 0, B0 = 1, a = NULL, b = NULL, c0 = NA, d0 = NA,
+ burnin = 1000, mcmc = 1000, thin = 1, verbose = 0,
+ seed = NA, beta.start = NA, P.start = NA,
+ 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)
+ n.arrival<- y + 1
+ NT <- max(n.arrival)
+ tot.comp <- n + sum(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
+ A0 <- trans.mat.prior(m=m, n=n, a=a, b=b)
+
+ if (k==1){
+ if (is.na(c0)||is.na(d0))
+ stop("You have to prior for lambda (c0 and d0) when there is no covariate.\n")
+ }
+ else{
+ c0 <- d0 <- 0
+ mvn.prior <- form.mvn.prior(b0, B0, k)
+ b0 <- mvn.prior[[1]]
+ B0 <- mvn.prior[[2]]
+ }
+
+ ## get marginal likelihood argument
+ marginal.likelihood <- match.arg(marginal.likelihood)
+ logmarglike <- NULL
+ chib <- 0
+ if (marginal.likelihood == "Chib95"){
+ chib <- 1
+ }
+
+ ## get initial values of tau from observed y
+ Pstart <- check.P(P.start, m, a=a, b=b)
+ betastart <- beta.change.start(beta.start, ns, k, formula, family=poisson, data)
+ if (k == 1){
+ betastart <- exp(betastart)
+ }
+ taustart <- tau.initial(y, tot.comp)
+ componentstart <- round(runif(tot.comp, 1, 5))
+
+ ## normal mixture weights
+ wr <- c(0.2294, 0.2590, 0.2480, 0.1525, 0.0472)
+ mr <- c(0.0982, -1.5320, -0.7433, 0.8303, -3.1428)
+ sr <- sqrt(c(0.2401, 1.1872, 0.3782, 0.1920, 3.2375))
+
+ ## call C++ code to draw sample
+ posterior <- .C("MCMCpoissonChange",
+ 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),
+ betastart = as.double(betastart),
+ Pstart = as.double(Pstart),
+ taustart = as.double(taustart),
+ componentstart = as.double(componentstart),
+ a = as.double(a),
+ b = as.double(b),
+ c0 = as.double(c0),
+ d0 = as.double(d0),
+ lecuyer = as.integer(lecuyer),
+ seedarray = as.integer(seed.array),
+ lecuyerstream = as.integer(lecuyer.stream),
+ b0data = as.double(b0),
+ B0data = as.double(B0),
+ A0data = as.double(A0),
+ logmarglikeholder = as.double(0.0),
+ wrin = as.double(wr),
+ mrin = as.double(mr),
+ srin = as.double(sr),
+ chib = as.integer(chib),
+
+ PACKAGE="MCMCpack")
+
+ ## get marginal likelihood if Chib95
+ if (marginal.likelihood == "Chib95"){
+ logmarglike <- posterior$logmarglikeholder
+ }
+
+ ## pull together matrix and build MCMC object to return
+ beta.holder <- matrix(posterior$betaout, nstore, )
+ 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(xnames, "_regime", i, sep = "")})
+ attr(output, "title") <- "MCMCpoissonChange 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, "prob.state") <- ps.holder/nstore
+ attr(output, "s.store") <- s.holder
+ attr(output, "P.store") <- P.holder
+ return(output)
+
+ }
diff --git a/R/MCMCpoissonChangepoint.R b/R/MCMCpoissonChangepoint.R
deleted file mode 100644
index 8c21797..0000000
--- a/R/MCMCpoissonChangepoint.R
+++ /dev/null
@@ -1,112 +0,0 @@
-##########################################################################
-## 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
-##
-## This software is distributed under the terms of the GNU GENERAL
-## PUBLIC LICENSE Version 2, June 1991. See the package LICENSE
-## file for more information.
-##
-## Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
-## Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
-## and Jong Hee Park
-##########################################################################
-
-"MCMCpoissonChangepoint"<-
- function(data, m = 1, c0 = NA, d0 = NA, a = NULL, b = NULL,
- 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()
-
- ## ## 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
- }
-
- 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
- }
-
- ## 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, "prob.state") <- ps.holder/(mcmc/thin)
- attr(output, "s.store") <- s.holder
- return(output)
-
- }## end of MCMC function
-
diff --git a/R/btsutil.R b/R/btsutil.R
index cd5f857..d04a58a 100644
--- a/R/btsutil.R
+++ b/R/btsutil.R
@@ -24,138 +24,224 @@
##############################################################
-## Helper functions for MCMCPoissonChangepoint()
+## Helper functions for MCMCpoissonChange and MCMCbinaryChange()
##############################################################
- ## switch a state vector into a matrix containing the number of states
- "switchg" <- function(s1){
- s <- max(s1)
- out <- matrix(0,s,s)
+## switch a state vector into a matrix containing the number of states
+"switchg" <- function(s1){
+ s <- max(s1)
+ out <- matrix(0,s,s)
+
+ ## all P(i,i+1) are 1
+ for (i in 1:(s-1)){
+ out[i,i+1] <- 1}
+
+ ## diagonal elements is (the number of occurrence - 1)
+ diag(out) <- table(s1)-1
+ return(out)
+}
- ## all P(i,i+1) are 1
- for (i in 1:(s-1)){
- out[i,i+1] <- 1}
+## "trans.mat.prior" makes a transition matrix
+"trans.mat.prior" <- function(m, n, a=NULL, b=NULL){
+ if (!is.null(a)|!is.null(b)){
+ a <- a
+ b <- b
+ }
+ else {
+ expected.duration <- round(n/(m+1))
+ b <- 0.1
+ a <- b*expected.duration
+ }
+ trans <- diag(1, m+1)
+ ## put a as diagonal elements except the last row
+ diag(trans)[1:m]<-rep(a, m)
+ ## put b in trans[i, i+1]
+ for (i in 1:m){trans[i, i+1]<-b}
+ return(trans)
+}
- ## diagonal elements is (the number of occurrence - 1)
- diag(out) <- table(s1)-1
- return(out)
- }
-
- ## "trans.mat.prior" makes a transition matrix
- "trans.mat.prior" <- function(m, n, a=NULL, b=NULL){
- if (!is.null(a)|!is.null(b)){
- a <- a
- b <- b
- }
- else {expected.duration <- round(n/(m+1))
- b <- 0.1
- a <- b*expected.duration
- }
- trans <- diag(1, m+1)
- # put a as diagonal elements except the last row
- diag(trans)[1:m]<-rep(a, m)
- # put b in trans[i, i+1]
- for (i in 1:m){trans[i, i+1]<-b}
- return(trans)
- }
-
- ## "plotState" draws a plot of posterior distribution of states
- "plotState" <-
- function (mcmcout, legend.control = NULL, cex = 0.8, lwd = 1.2)
- {
+## "plotState" draws a plot of posterior distribution of states
+"plotState" <-
+ function (mcmcout, main="Posterior Regime Probability", ylab=expression(paste("Pr(", S[t], "= k |", Y[t], ")")),
+ legend.control = NULL, cex = 0.8, lwd = 1.2, start=1)
+ {
out <- attr(mcmcout, "prob.state")
y <- attr(mcmcout, "y")
m <- attr(mcmcout, "m")
+
if (!is.ts(y))
- y <- ts(y)
+ y <- ts(y, start)
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], ")")))
+
+ plot(start, 0, xlim = range(time.frame), ylim = c(0, 1), type = "n",
+ main = main, xlab = "Time", cex = cex, lwd = lwd,
+ ylab = ylab, axes=F)
+ axis(1, tick = FALSE, col="darkgrey")
+ axis(2, tick = FALSE, col="darkgrey")
+ for (i in 1:length(axTicks(1))) lines(c(axTicks(1)[i], axTicks(1)[i]), c(0,1), col="darkgrey")
+ for (i in 1:length(axTicks(2))) lines(c(axTicks(2)[i], max(axTicks(1))), c(axTicks(2)[i], axTicks(2)[i]), col="darkgrey")
+
for (i in 1:(m + 1)) points(time.frame, out[, i], type = "o",
- lty = i, lwd = lwd, col = i, cex = cex)
+ 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")
+ 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")
}
-}
+ }
+
- ## "plotChangepoint" draws a plot of posterior changepoint probability
- "plotChangepoint" <-
- function (mcmcout, xlab = "Time", ylab = "", verbose = FALSE)
- {
+## "plotChangepoint" draws a plot of posterior changepoint probability
+## Thanks to Patrick Brandt for providing the idea of overlaying.
+"plotChangepoint" <-
+ function (mcmcout, main="Posterior Density of Regime Change Probabilities", xlab = "Time", ylab = "",
+ verbose = FALSE, start=1, overlay=FALSE)
+ {
out <- attr(mcmcout, "prob.state")
y <- attr(mcmcout, "y")
m <- attr(mcmcout, "m")
+ if(overlay==FALSE){
+ par(mfrow = c(m, 1), mar = c(2, 4, 1, 1))
+ }
if (!is.ts(y))
- y <- ts(y)
+ y <- ts(y, start)
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")
+ pr.st <- c(0, diff(out[, (m + 1)]))
+ pr.st[pr.st<0] <- 0
+ plot(time.frame, pr.st, type = "h", lwd=2, main = main, xlab = xlab, ylab = ylab, axes=F)
+ axis(1, tick = FALSE, col="darkgrey")
+ axis(2, tick = FALSE, col="darkgrey")
+ for (i in 1:length(axTicks(1))) lines(c(axTicks(1)[i], axTicks(1)[i]), c(0, max(axTicks(2))), col="darkgrey")
+ for (i in 1:length(axTicks(2))) lines(c(axTicks(2)[i], max(axTicks(1))), c(axTicks(2)[i], axTicks(2)[i]), col="darkgrey")
+ cp <- which(cumsum(pr.st) > 0.5)[1] - 1
+ lines(c(cp + time.frame[1], cp + time.frame[1]), c(0, max(axTicks(2))), 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")
- }
- 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 <- 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", lwd=2, main = "", xlab = xlab, ylab = ylab, col="black", axes=FALSE)
+ axis(1, tick = FALSE, col="darkgrey")
+ axis(2, tick = FALSE, col="darkgrey")
+ for (k in 1:length(axTicks(1))) {lines(c(axTicks(1)[k], axTicks(1)[k]), c(0, max(axTicks(2))), col="darkgrey")}
+ for (k in 1:length(axTicks(2))) {lines(c(axTicks(2)[k], max(axTicks(1))), c(axTicks(2)[k], axTicks(2)[k]),
+ col="darkgrey")}
+ cp[i - 1] <- which(cumsum(pr.st) > 0.5)[1] - 1
+ lines(c(cp[i - 1] + time.frame[1], cp[i - 1] + time.frame[1]), c(0, max(axTicks(2))), lty = 3, col = "red")
+ }
+ pr.st <- c(0, diff(out[, (m + 1)]))
+ pr.st[pr.st<0] <- 0
+ plot(time.frame, pr.st, type = "h", lwd=2, main = main, xlab = xlab, ylab = ylab, col="black", axes=FALSE)
+ axis(1, tick = FALSE, col="darkgrey")
+ axis(2, tick = FALSE, col="darkgrey")
+ for (k in 1:length(axTicks(1))) {lines(c(axTicks(1)[k], axTicks(1)[k]), c(0, max(axTicks(2))), col="darkgrey")}
+ for (k in 1:length(axTicks(2))) {lines(c(axTicks(2)[k], max(axTicks(1))), c(axTicks(2)[k], axTicks(2)[k]),
+ col="darkgrey")}
+ cp[m] <- which(cumsum(pr.st) > 0.5)[1] - 1
+ lines(c(cp[m] + time.frame[1], cp[m] + time.frame[1]), c(0, max(axTicks(2))), lty = 3, col = "red")
}
+
cp.means <- rep(NA, m + 1)
cp.start <- c(1, cp + 1)
cp.end <- c(cp, length(y))
- if (verbose == TRUE){
- 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("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\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)
+
+ if (verbose == TRUE){
+ 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("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\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)
+}
+
- ## 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)
+## initial values of tau in MCMCpoissonChange
+"tau.initial" <- function(y, tot.comp){
+ tau <- rep(NA, tot.comp)
+ lambda.t <- 0.1
+ count <- 0
+ for (t in 1:length(y)){
+ nt <- y[t]
+ if (nt==0) {
+ taut <- 1 + rexp(1, lambda.t)
+ count <- count + nt + 1
+ }
+ else{
+ ut <- runif(nt)
+ uorder <- c(0, sort(ut))
+ tau.tj <- diff(uorder)
+ sum.tau.tj <- sum(tau.tj)
+ tau.last<- 1 - sum.tau.tj + rexp(1, y[t])
+ count <- count + nt + 1
+ taut <- c(tau.tj, tau.last)
}
+ tau[(count-nt):count] <- taut
+ }
+ return(tau)
+}
+
+## beta starting values in MCMCpoissonChange()
+"beta.change.start"<- function (beta.start, ns, k, formula, family, data){
+ ## if a user does not specify beta.start, use a coefficient vector from mle
+ if (is.na(beta.start[1])) {
+ b0 <- coef(glm(formula, family = family, data = data))
+ beta.start <- matrix(rep(b0, ns), ns, k, byrow=TRUE)
+ }
+ ## if beta.start is scalar or k by 1 vector, repeat this
+ else if (is.null(dim(beta.start))&&length(beta.start)<=k) {
+ beta.start <- beta.start * matrix(1, ns, k)
+ ## this alternates beta.start if beta.start is not a scalar
+ }
+ ## if the length of beta.start is same to ns*k, make this as a matrix
+ else if (is.null(dim(beta.start))&&length(beta.start)==ns*k) {
+ beta.start <- matrix(beta.start, ns, k)
+ }
+ else if (is.null(dim(beta.start))&&length(beta.start)>=k) {
+ cat("Error: Starting value for beta not conformable.\n")
+ stop("Please respecify and call ", calling.function(),
+ " again.\n", call. = FALSE)
+ }
+ ## else, report an error message and stop
+ else if (!all(dim(beta.start) == c(ns, k))) {
+ cat("Error: Starting value for beta not conformable.\n")
+ stop("Please respecify and call ", calling.function(),
+ " again.\n", call. = FALSE)
+ }
+ return(beta.start)
+}
diff --git a/man/MCMCpoissonChangepoint.Rd b/man/MCMCbinaryChange.Rd
similarity index 52%
rename from man/MCMCpoissonChangepoint.Rd
rename to man/MCMCbinaryChange.Rd
index 8e1c4f3..ba2a180 100644
--- a/man/MCMCpoissonChangepoint.Rd
+++ b/man/MCMCbinaryChange.Rd
@@ -1,124 +1,131 @@
-\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}}}
+\name{MCMCbinaryChange}
+\alias{MCMCbinaryChange}
+
+\title{Markov Chain Monte Carlo for a Binary Multiple Changepoint Model}
+\description{
+ This function generates a sample from the posterior distribution
+ of a binary 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{MCMCbinaryChange(data, m = 1, c0 = 1, d0 = 1, a = NULL, b = NULL,
+ burnin = 10000, mcmc = 10000, thin = 1, verbose = 0,
+ seed = NA, phi.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 shape1 parameter for Beta prior on \eqn{\phi}{phi}
+ (the mean).}
+
+ \item{d0}{\eqn{d_0}{d0} is the shape2 parameter for Beta prior on \eqn{\phi}{phi}
+ (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{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{phi.start}{The starting values for the mean. The default value of NA will use draws from the Uniform distribution.}
+
+ \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{MCMCbinaryChange} simulates from the posterior distribution of
+ a binary model with multiple changepoints.
+
+ The model takes the following form:
+ \deqn{Y_t \sim \mathcal{B}ernoulli(\phi_i),\;\; i = 1, \ldots, k}{Y_t ~ Bernoulli(phi_i), i = 1,...,k.}
+ Where \eqn{k}{k} is the number of states.
+
+ We assume Beta priors for \eqn{\phi_{i}}{phi_i} and for transition probabilities:
+ \deqn{\phi_i \sim \mathcal{B}eta(c_0, d_0)}{phi_i ~ Beta(c0, d0)}
+ And:
+ \deqn{p_{mm} \sim \mathcal{B}eta{a}{b},\;\; m = 1, \ldots, k}{p_mm ~ Beta(a, b), m = 1,...,M.}
+ Where \eqn{M}{M} is the number of states.
+ }
+
+\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{
+ set.seed(19173)
+ true.phi<- c(0.5, 0.8, 0.4)
+
+ ## two breaks at c(80, 180)
+ y1 <- rbinom(80, 1, true.phi[1])
+ y2 <- rbinom(100, 1, true.phi[2])
+ y3 <- rbinom(120, 1, true.phi[3])
+ y <- as.ts(c(y1, y2, y3))
+
+ model1 <- MCMCbinaryChange(y, m=1, c0=2, d0=2, mcmc=1000, burnin=1000, verbose=500, marginal.likelihood = "Chib95")
+ model2 <- MCMCbinaryChange(y, m=2, c0=2, d0=2, mcmc=1000, burnin=1000, verbose=500, marginal.likelihood = "Chib95")
+ model3 <- MCMCbinaryChange(y, m=3, c0=2, d0=2, mcmc=1000, burnin=1000, verbose=500, marginal.likelihood = "Chib95")
+ model4 <- MCMCbinaryChange(y, m=4, c0=2, d0=2, mcmc=1000, burnin=1000, verbose=500, marginal.likelihood = "Chib95")
+ model5 <- MCMCbinaryChange(y, m=5, c0=2, d0=2, mcmc=1000, burnin=1000, verbose=500, marginal.likelihood = "Chib95")
+
+ print(BayesFactor(model1, model2, model3, model4, model5))
+
+ ## plot two plots in one screen
+ 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)
+
+ }
+}
+
+\keyword{models}
+
+\seealso{\code{\link{MCMCpoissonChange}},\code{\link{plotState}}, \code{\link{plotChangepoint}}}
diff --git a/man/MCMCpoissonChange.Rd b/man/MCMCpoissonChange.Rd
new file mode 100644
index 0000000..697fedd
--- /dev/null
+++ b/man/MCMCpoissonChange.Rd
@@ -0,0 +1,167 @@
+\name{MCMCpoissonChange}
+\alias{MCMCpoissonChange}
+
+\title{Markov Chain Monte Carlo for a Poisson Regression Changepoint Model}
+\description{
+ This function generates a sample from the posterior distribution
+ of a Poisson regression 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{MCMCpoissonChange(
+ formula, data = parent.frame(), m = 1,
+ b0 = 0, B0 = 1, a = NULL, b = NULL, c0 = NA, d0 = NA,
+ burnin = 1000, mcmc = 1000, thin = 1, verbose = 0,
+ seed = NA, beta.start = NA, P.start = NA,
+ marginal.likelihood = c("none", "Chib95"), ...)}
+
+\arguments{
+ \item{formula}{Model formula.}
+
+ \item{data}{Data frame.}
+
+ \item{m}{The number of changepoints.}
+
+ \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{c0}{\eqn{c_0}{c0} is the shape parameter for Gamma prior on \eqn{\lambda}{lambda}
+ (the mean). When there is no covariate, this should be provided by users. No default value is provided.}
+
+ \item{d0}{\eqn{d_0}{d0} is the scale parameter for Gamma prior on \eqn{\lambda}{lambda}
+ (the mean). When there is no covariate, this should be provided by users. No default value is provided.}
+
+ \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{beta.start}{The starting values for the beta vector. This can either be a scalar or a column vector with dimension equal to the number of betas. 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 betas. When there is no covariate, the log value of means should be used.}
+
+ \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{MCMCpoissonChange} simulates from the posterior distribution of
+ a Poisson regression model with multiple changepoints.
+
+ The model takes the following form:
+ \deqn{y_t \sim \mathcal{P}oisson(\mu_t)}{y_t ~ Poisson(mu_t)}
+ \deqn{\mu_t = x_t ' \beta_m,\;\; m = 1, \ldots, M}{mu_t = x_t'beta_m, m = 1,...,M.}
+ Where \eqn{M}{M} is the number of states and \eqn{\beta_m}{beta_m} is paramters when a state is \eqn{m}{m} at \eqn{t}{t}.
+
+ We assume Gaussian distribution for prior of \eqn{\beta}{beta}:
+ \deqn{\beta_m \sim \mathcal{N}(b_0,B_0^{-1}),\;\; m = 1, \ldots, M}{beta_m ~ N(b0,B0^(-1)), m = 1,...,M.}
+ And:
+ \deqn{p_{mm} \sim \mathcal{B}eta{a}{b},\;\; m = 1, \ldots, k}{p_mm ~ Beta(a, b), m = 1,...,M.}
+ Where \eqn{M}{M} is the number of states.
+ }
+
+\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(11119)
+ n <- 150
+ x1 <- runif(n, 0, 0.5)
+ true.beta1 <- c(1, 1)
+ true.beta2 <- c(1, -2)
+ true.beta3 <- c(1, 2)
+
+ ## 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)
+ mu3 <- exp(1 + x1[true.s==3]*2)
+
+ y <- as.ts(c(rpois(n/3, mu1), rpois(n/3, mu2), rpois(n/3, mu3)))
+ formula = y ~ x1
+
+ model1 <- MCMCpoissonChange(formula, m=1,
+ mcmc = 1000, burnin = 1000, verbose = 500,
+ b0 = rep(0, 2), B0 = 5*diag(2), marginal.likelihood = "Chib95")
+ model2 <- MCMCpoissonChange(formula, m=2,
+ mcmc = 1000, burnin = 1000, verbose = 500,
+ b0 = rep(0, 2), B0 = 5*diag(2), marginal.likelihood = "Chib95")
+ model3 <- MCMCpoissonChange(formula, m=3,
+ mcmc = 1000, burnin = 1000, verbose = 500,
+ b0 = rep(0, 2), B0 = 5*diag(2), marginal.likelihood = "Chib95")
+ model4 <- MCMCpoissonChange(formula, m=4,
+ mcmc = 1000, burnin = 1000, verbose = 500,
+ b0 = rep(0, 2), B0 = 5*diag(2), marginal.likelihood = "Chib95")
+ model5 <- MCMCpoissonChange(formula, m=5,
+ mcmc = 1000, burnin = 1000, verbose = 500,
+ b0 = rep(0, 2), B0 = 5*diag(2), marginal.likelihood = "Chib95")
+
+ print(BayesFactor(model1, model2, model3, model4, model5))
+
+ ## 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
+ model2.1 <- MCMCpoissonChange(y ~ 1, m = 2, c0 = 2, d0 = 1,
+ mcmc = 1000, burnin = 1000, verbose = 500,
+ marginal.likelihood = "Chib95")
+ print(BayesFactor(model2, model2.1))
+ }
+}
+
+\keyword{models}
+
+\seealso{\code{\link{MCMCbinaryChange}}, \code{\link{plotState}}, \code{\link{plotChangepoint}}}
diff --git a/man/plotChangepoint.Rd b/man/plotChangepoint.Rd
index ecc8482..1f58d92 100644
--- a/man/plotChangepoint.Rd
+++ b/man/plotChangepoint.Rd
@@ -1,26 +1,32 @@
\name{plotChangepoint}
\alias{plotChangepoint}
-\title{Changepoint Location Plots}
-\description{Plot the posterior distribution over the location of the changepoints.}
+\title{Posterior Density of Regime Change Plot}
+\description{Plot the posterior density of regime change.}
\usage{
- plotChangepoint(mcmcout, xlab="Time", ylab="", verbose = FALSE)
+ plotChangepoint(mcmcout, main="Posterior Density of Regime Change Probabilities",
+ xlab = "Time", ylab = "", verbose = FALSE, start=1, overlay=FALSE)
}
\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{main}{Title of the plot}
+
\item{xlab}{Label for the x-axis.}
\item{ylab}{Label for the y-axis.}
\item{verbose}{If verbose=TRUE, expected changepoints are printed.}
+\item{start}{The time of the first observation to be shown in the time series plot.}
+
+\item{overlay}{If overlay=TRUE, the probability of each regime change is drawn separately, which will be useful to draw multiple plots in one screen. See the example in \code{MCMCpoissonChange}. Otherwise, multiple plots of regime change probabilities will be drawn.}
}
\author{Jong Hee Park, \email{jhp at uchicago.edu}, \url{http://home.uchicago.edu/~jhp/}.}
\keyword{hplot}
-\seealso{\code{\link{MCMCpoissonChangepoint}}}
+\seealso{\code{\link{MCMCpoissonChange}}, \code{\link{MCMCbinaryChange}}}
diff --git a/man/plotState.Rd b/man/plotState.Rd
index 1d21369..2c55e49 100644
--- a/man/plotState.Rd
+++ b/man/plotState.Rd
@@ -4,23 +4,30 @@
\description{Plot the posterior probability that each time point is in each state.}
\usage{
- plotState(mcmcout, legend.control=NULL, cex=0.8, lwd=1.2)
+ plotState(mcmcout, main="Posterior Regime Probability", ylab=expression(paste("Pr(", S[t], "= k |", Y[t], ")")),
+ legend.control = NULL, cex = 0.8, lwd = 1.2, start=1)
}
\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{main}{Title of the plot.}
+
+\item{ylab}{Label for the y-axis.}
+
\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.}
\item{lwd}{Line width parameter.}
+\item{start}{The time of the first observation to be shown in the time series plot.}
+
}
\author{Jong Hee Park, \email{jhp at uchicago.edu}, \url{http://home.uchicago.edu/~jhp/}.}
\keyword{hplot}
-\seealso{\code{\link{MCMCpoissonChangepoint}}}
+\seealso{\code{\link{MCMCpoissonChange}}, \code{\link{MCMCbinaryChange}}}
diff --git a/src/MCMCbinaryChange.cc b/src/MCMCbinaryChange.cc
new file mode 100644
index 0000000..87d8a55
--- /dev/null
+++ b/src/MCMCbinaryChange.cc
@@ -0,0 +1,363 @@
+// MCMCbinaryChange.cc is C++ code to estimate a binary changepoint model
+// with a beta prior
+//
+// Jong Hee Park
+// Dept. of Political Science
+// University of Chicago
+// jhp at uchicago.edu
+//
+// 03/03/2009
+
+#ifndef MCMCBINARYCHANGE_CC
+#define MCMCBINARYCHANGE_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>
+
+
+template <typename RNGTYPE>
+// bianry state sampler
+Matrix<> binary_state_sampler(rng<RNGTYPE>& stream,
+ const int m,
+ const Matrix<>& Y,
+ const Matrix<>& theta,
+ const Matrix<>& P){
+
+ const int ns = m + 1;
+ const int n = Y.rows();
+ Matrix<> F = Matrix<>(n, ns);
+ Matrix<> pr1 = Matrix<>(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] = dbinom(yt, 1, theta[j]);
+ }
+ 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;
+
+}
+
+////////////////////////////////////////////
+// MCMCbinaryChangepoint implementation.
+////////////////////////////////////////////
+template <typename RNGTYPE>
+void MCMCbinaryChange_impl(rng<RNGTYPE>& stream, const Matrix<>& Y,
+ Matrix<>& phi, 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<>& phi_store, Matrix<>& P_store,
+ Matrix<>& ps_store, Matrix<>& s_store,
+ double& logmarglike)
+{
+ const unsigned int tot_iter = burnin + mcmc;
+ const unsigned int nstore = mcmc / thin;
+ const int n = Y.rows();
+ const int ns = m + 1;
+
+ //MCMC loop
+ unsigned int count = 0;
+ for (int iter = 0; iter < tot_iter; ++iter){
+
+ //////////////////////
+ // 1. Sample s
+ //////////////////////
+ Matrix<> Sout = binary_state_sampler(stream, m, Y, phi, P);
+ Matrix<> s = Sout(_, 0);
+ Matrix<> ps(n, ns);
+ for (int j = 0; j<ns; ++j){
+ ps(_,j) = Sout(_,j+1);
+ }
+
+ //////////////////////
+ // 2. Sample phi
+ //////////////////////
+ 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)) {
+ addN[j] = addN[j] + 1;
+ addY[j] = addY[j] + Y[i];
+ }
+ }
+ double c1 = addY[j] + c0;
+ double d1 = addN[j] - addY[j] + d0;
+ phi[j] = stream.rbeta(c1, d1);
+
+ }
+ //////////////////////
+ // 3. 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) + addN[j] - 1;
+ shape2 = A0(j,j+1) + 1;
+ P(j,j) = stream.rbeta(shape1, shape2);
+ P(j,j+1) = 1 - P(j,j);
+ }
+
+ if (iter >= burnin && ((iter % thin)==0)){
+ for (int i=0; i<ns; ++i)
+ phi_store(count,i) = phi[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,_);
+
+ ++count;
+
+ }
+
+
+ if(verbose > 0 && iter % verbose == 0){
+ Rprintf("\n\n MCMCbinaryChange iteration %i of %i", (iter+1), tot_iter);
+ for (int j = 0;j<ns; ++j){
+ Rprintf("\n The number of observations in state %i is %10.5f", j + 1, addN[j]);
+ }
+ for (int j=0; j<ns; ++j)
+ Rprintf("\n phi in state %i is %10.5f", j + 1, phi[j]);
+ }
+
+ R_CheckUserInterrupt();
+
+ }// end MCMC loop
+
+ if(chib ==1){
+
+ Matrix<> phi_st = meanc(phi_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];
+ }
+
+ //////////////////////
+ // phi
+ //////////////////////
+ Matrix<> density_phi(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];
+ }
+ }
+ double c1 = addY[j] + c0;
+ double d1 = addN[j] - addY[j] + d0;
+ density_phi(iter, j) = dbeta(phi_st[j], c1, d1);
+ }
+
+ }
+ double pdf_phi = log(prod(meanc(density_phi)));
+
+ //////////////////////
+ // P
+ //////////////////////
+ Matrix<> density_P(nstore, ns);
+ for (int iter = 0; iter < nstore; ++iter){
+ Matrix<> Sout = binary_state_sampler(stream, m, Y, phi_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)) {
+ 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(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] = dbinom(yt, 1, phi_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);
+ }
+
+ double loglike = sum(log(like));
+
+ //////////////////////
+ // prior
+ //////////////////////
+ Matrix<> density_phi_prior(ns, 1);
+ Matrix<> density_P_prior(ns, 1);
+
+ for (int j=0; j<ns ; ++j){
+ density_phi_prior[j] = log(dbeta(phi_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_phi_prior) + sum(density_P_prior);
+ logmarglike = (loglike + logprior) - (pdf_phi + pdf_P);
+ }
+}
+////////////////////////////////////////////
+// Start MCMCbinaryChangepoint function
+///////////////////////////////////////////
+extern "C"{
+ void MCMCbinaryChange(double *phiout, 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 *verbose,
+ const int *uselecuyer, const int *seedarray, const int *lecuyerstream,
+ const double *phistart, 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 <> Y(*Yrow, *Ycol, Ydata);
+ const unsigned int tot_iter = *burnin + *mcmc;
+ const unsigned int nstore = *mcmc / *thin;
+ const int n = Y.rows();
+ const int ns = *m + 1;
+
+ // generate starting values
+ Matrix <> phi(ns, 1, phistart);
+ const Matrix <> A0(ns, ns, A0data);
+ Matrix <> P(ns, ns, Pstart);
+
+ double logmarglike;
+
+ // storage matrices
+ Matrix<> phi_store(nstore, ns);
+ Matrix<> P_store(nstore, ns*ns);
+ Matrix<> ps_store(n, ns);
+ Matrix<> s_store(nstore, n);
+
+ MCMCPACK_PASSRNG2MODEL(MCMCbinaryChange_impl, Y,
+ phi, P, A0, *m, *c0, *d0,
+ *burnin, *mcmc, *thin, *verbose, *chib,
+ phi_store, P_store, ps_store, s_store,
+ logmarglike)
+
+ logmarglikeholder[0] = logmarglike;
+
+ // return output
+ for (int i = 0; i<(nstore*ns); ++i){
+ phiout[i] = phi_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
new file mode 100644
index 0000000..7c071b5
--- /dev/null
+++ b/src/MCMCpoissonChange.cc
@@ -0,0 +1,937 @@
+// MCMCpoissonRegChange.cc is C++ code to estimate a Poisson changepoint model
+//
+// Jong Hee Park
+// Dept. of Political Science
+// University of Chicago
+// jhp at uchicago.edu
+//
+// 07/06/2007
+// 12/20/2007 included in MCMCpack
+
+#ifndef MCMCPOISSONCHANGE_CC
+#define MCMCPOISSONCHANGE_CC
+
+#include<vector>
+#include<algorithm>
+
+#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;
+
+//tau and component sampler
+template <typename RNGTYPE>
+Matrix<> tau_comp_sampler(rng<RNGTYPE>& stream,
+ const int m,
+ const int totcomp,
+ const Matrix<>& Y,
+ const Matrix<>& X,
+ const Matrix<>& wr,
+ const Matrix<>& mr,
+ const Matrix<>& sr,
+ const Matrix<>& beta,
+ const Matrix<>& s){
+
+ // itialize
+ const int ns = m + 1;
+ const int n = Y.rows();
+ const int k = X.cols();
+ Matrix<int> component(totcomp, 1);
+ Matrix<> tau(totcomp, 1);
+ Matrix<> post_taut_mat(totcomp, 5);
+ Matrix<> post_tau_mat(totcomp, 5);
+ Matrix<> norm_post_tau(totcomp, 5);
+ Matrix<> cumsum_post_tau(totcomp, 5);
+ Matrix<> xt(1, k);
+ Matrix<> wr_mat = eye(5);
+ for (int i=0; i<5; ++i) {
+ wr_mat(i,i) = wr[i];
+ }
+
+ int tau_start = 0;
+ int tau_end = 0;
+
+ for(int t=0; t<n; ++t){
+ int yt = (int)Y[t];
+ xt = X(t,_);
+ tau_start = tau_end;
+ tau_end = tau_start + yt + 1;
+
+ int st = (int)s[t];
+ Matrix<> mu_t = exp(xt* ::t(beta(st-1,_)));
+ double mut = mu_t[0];
+
+ if (yt == 0){
+ double tau_double = 1 + stream.rexp(mut);
+ tau[tau_end - 1] = tau_double;
+ for (int h=0; h<5; ++h){
+ double first = 1/(sr[h]*tau_double);
+ double second = (log(tau_double) + log(mut) - mr[h])/sr[h];
+ post_taut_mat(tau_end-1, h) = first*exp(-0.5*second*second);
+ }
+ }
+
+ else {
+ Matrix<> ut = stream.runif(yt, 1);
+ Matrix<> sort_ut = ::sort(ut);
+ Matrix<> tau_tj(yt, 1);
+ for(int i=1; i<yt; ++i){
+ tau_tj[i] = sort_ut[i] - sort_ut[i-1];
+ }
+ tau_tj[0] = sort_ut[0];
+
+ double sum_tau_tj = sum(tau_tj);
+ double tau_last = 1 - sum_tau_tj + stream.rexp(mut);
+
+ Matrix<> tau_mat(yt+1, 1);
+ tau_mat(0, 0, yt-1, 0) = tau_tj(_,0);
+ tau_mat[yt] = tau_last;
+
+ for (int i = 0; i<(yt+1); ++i){
+ tau[i + tau_start] = tau_mat[i];
+ for (int h=0; h<5; ++h){
+ double first = 1/(sr[h]*tau_mat[i]);
+ double second = (log(tau_mat[i]) + log(mut) - mr[h])/sr[h];
+ post_taut_mat(i+tau_start, h) = first*exp(-0.5*second*second);
+ }
+ }
+ }
+
+ }
+ post_tau_mat = post_taut_mat*wr_mat;
+
+ for(int i = 0; i<totcomp ; ++i){
+ norm_post_tau(i,_) = post_tau_mat(i,_)/sum(post_tau_mat(i,_));
+ cumsum_post_tau(i,0) = norm_post_tau(i,0);
+ for (int j=1; j<5; ++j){
+ cumsum_post_tau(i,j) = cumsum_post_tau(i,j-1) + norm_post_tau(i,j);
+ }
+ double U = stream.runif();
+ if (U < cumsum_post_tau(i,0)){
+ component[i] = 1;
+ }
+ if (cumsum_post_tau(i,0)<=U&&U<cumsum_post_tau(i,1)){
+ component[i] = 2;
+ }
+ if (cumsum_post_tau(i,1)<=U&&U<cumsum_post_tau(i,2)){
+ component[i] = 3;
+ }
+ if (cumsum_post_tau(i,2)<=U&&U<cumsum_post_tau(i,3)){
+ component[i] = 4;
+ }
+ if (cumsum_post_tau(i,3)<=U&&U<=1){
+ component[i] = 5;
+ }
+ }
+
+ Matrix<> TAUout(totcomp, 2);
+ TAUout(_, 0) = tau(_, 0);
+ TAUout(_, 1) = component(_, 0);
+
+ return TAUout;
+
+}
+
+template <typename RNGTYPE>
+Matrix<> poisson_state_sampler(rng<RNGTYPE>& stream,
+ const int& m,
+ const Matrix<>& Y,
+ const Matrix<>& lambda,
+ 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){
+ 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);
+ for (int j=0; j<ns ; ++j){
+ F(t,j) = pstyt(j);
+ }
+ }
+ 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;
+
+}
+
+template <typename RNGTYPE>
+Matrix<> poisson_reg_state_sampler(rng<RNGTYPE>& stream,
+ const int m,
+ const Matrix<>& Y,
+ const Matrix<>& X,
+ const Matrix<>& beta,
+ 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){
+ int yt = (int) Y[t];
+ Matrix<> lambda = exp(X(t,_)*::t(beta));
+ 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;
+
+ 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;
+}
+
+
+////////////////////////////////////////////
+// MCMCpoissonChangepoint implementation.
+////////////////////////////////////////////
+template <typename RNGTYPE>
+void MCMCpoissonChangepoint_impl(rng<RNGTYPE>& stream,
+ double *betaout,
+ double *Pout,
+ double *psout,
+ double *sout,
+ const double *Ydata,
+ const int *Yrow,
+ const int *Ycol,
+ const int *m,
+ const double *c0,
+ const double *d0,
+ const int *burnin,
+ const int *mcmc,
+ const int *thin,
+ const int *verbose,
+ const double *betastart,
+ const double *Pstart,
+ const double *a,
+ const double *b,
+ const double *A0data,
+ double *logmarglikeholder,
+ const int *chib)
+
+{
+ const Matrix <> Y(*Yrow, *Ycol, Ydata);
+ const int tot_iter = *burnin + *mcmc;
+ const int nstore = *mcmc / *thin;
+ const int n = Y.rows();
+ const int ns = *m + 1;
+ const Matrix <> A0(ns, ns, A0data);
+
+ Matrix <> lambda(ns, 1, betastart);
+ Matrix <> P(ns, ns, Pstart);
+
+ Matrix<> lambda_store(nstore, ns);
+ Matrix<> P_store(nstore, ns*ns);
+ Matrix<> ps_store(n, ns);
+ Matrix<> s_store(nstore, n);
+
+ //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)) {
+ addN[j] = addN[j] + 1;
+ addY[j] = addY[j] + Y[i];
+ }
+ }
+ 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;
+
+ for (int j =0; j< (ns-1); ++j){
+ shape1 = A0(j,j) + addN[j] - 1;
+ shape2 = A0(j,j+1) + 1;
+ P(j,j) = stream.rbeta(shape1, shape2);
+ P(j,j+1) = 1 - P(j,j);
+ }
+
+ 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,_);
+
+ ++count;
+
+ }
+
+ if(*verbose > 0 && iter % *verbose == 0){
+ Rprintf("\n\n MCMCpoissonChange iteration %i of %i", (iter+1), tot_iter);
+ for (int j = 0;j<ns; ++j){
+ Rprintf("\n The number of observations in state %i is %10.5f", j + 1, addN[j]);
+ }
+ for (int j=0; j<ns; ++j)
+ Rprintf("\n lambda in state %i is %10.5f", j + 1, lambda[j]);
+ }
+
+ R_CheckUserInterrupt();
+
+ }
+
+ 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];
+ }
+
+ //////////////////////
+ // 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];
+ }
+ }
+ double c1 = addY[j] + *c0;
+ double d1 = addN[j] + *d0;
+ density_lambda(iter, j) = dgamma(lambda_st[j], c1, 1/d1);
+ }
+
+ }
+ double pdf_lambda = log(prod(meanc(density_lambda)));
+
+ //////////////////////
+ // 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)) {
+ 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(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);
+ }
+ 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);
+ }
+
+ double loglike = sum(log(like));
+
+ //////////////////////
+ // prior
+ //////////////////////
+ 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
+ const double logprior = sum(density_lambda_prior) + sum(density_P_prior);
+ const double logmarglike = (loglike + logprior) - (pdf_lambda + pdf_P);
+
+ logmarglikeholder[0] = logmarglike;
+ }
+
+ R_CheckUserInterrupt();
+
+ for (int i = 0; i<(nstore*ns); ++i){
+ betaout[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];
+ }
+
+
+}
+
+////////////////////////////////////////////
+// MCMCpoissonRegChange implementation.
+////////////////////////////////////////////
+template <typename RNGTYPE>
+void MCMCpoissonRegChangepoint_impl(rng<RNGTYPE>& stream,
+ 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 double *betastart,
+ const double *Pstart,
+ const double *taustart,
+ const double *componentstart,
+ const double *a,
+ const double *b,
+ const double *b0data,
+ const double *B0data,
+ const double *A0data,
+ double *logmarglikeholder,
+ const double* wrin,
+ const double* mrin,
+ const double* srin,
+ const int *chib)
+{
+ const Matrix <> Y(*Yrow, *Ycol, Ydata);
+ const Matrix <> X(*Xrow, *Xcol, Xdata);
+
+ const int tot_iter = *burnin + *mcmc;
+ const int nstore = *mcmc / *thin;
+ const int n = Y.rows();
+ const int k = X.cols();
+ const int ns = *m + 1;
+ const int totcomp = n + (int) sum(Y);
+ const Matrix <> b0(k, 1, b0data);
+ const Matrix <> B0(k, k, B0data);
+ const Matrix <> B0inv = invpd(B0);
+ Matrix <> wr(5, 1, wrin);
+ Matrix <> mr(5, 1, mrin);
+ Matrix <> sr(5, 1, srin);
+ const Matrix <> A0(ns, ns, A0data);
+
+ Matrix <> beta(ns, k, betastart);
+ Matrix <> tau(totcomp, 1, taustart);
+ Matrix <> component(totcomp, 1, componentstart);
+ Matrix <> P(ns, ns, Pstart);
+
+ Matrix<> beta_store(nstore, ns*k);
+ Matrix<> P_store(nstore, ns*ns);
+ Matrix<> ps_store(n, ns);
+ Matrix<> s_store(nstore, n);
+ Matrix<> component_store(nstore, totcomp);
+ Matrix<> tau_store(nstore, totcomp);
+
+ Matrix<> y_tilde(n ,1);
+ Matrix<> Sigma_inv_sum(n, 1);
+
+ //MCMC loop
+ int count = 0;
+ for (int iter = 0; iter < tot_iter; ++iter){
+
+ int y_count = 0;
+ for (int t = 0; t<n ; ++t){
+ int yt = (int) Y[t];
+ y_count = y_count + yt + 1;
+
+ Matrix<> Yt_over_Sigma(yt + 1, 1);
+ Matrix<> Sigma_inv(yt + 1, 1);
+
+ for(int j = (y_count - yt - 1); j< y_count; ++j){
+ int jone = (int) component[j] - 1 ; //zero base in C!
+ Sigma_inv[j-(y_count-yt-1)] = 1/(sr[jone]*sr[jone]);
+ Yt_over_Sigma[j-(y_count-yt-1)] = (log(tau[j])- mr[jone])*Sigma_inv[j-(y_count-yt-1)];
+ }
+ y_tilde[t] = sum(Yt_over_Sigma);
+ Sigma_inv_sum[t] = sum(Sigma_inv);
+ }
+ //////////////////////
+ // 1. Sample s
+ //////////////////////
+ Matrix <> Sout = poisson_reg_state_sampler(stream, *m, Y, X, beta, P);
+ Matrix <> s = Sout(_, 0);
+ Matrix <> ps(n, ns);
+ for (int j=0; j<ns; ++j){
+ ps(_,j) = Sout(_,j+1);
+ }
+
+ //////////////////////
+ // 2. 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];
+ Matrix<> yj = y_tilde((beta_count - nstate[j]), 0, (beta_count - 1), 0);
+ Matrix<> Xj = X((beta_count - nstate[j]), 0, (beta_count - 1), k-1);
+ Matrix<> wi = Sigma_inv_sum((beta_count - nstate[j]), 0, (beta_count - 1), 0);
+ Matrix<> Xwj(nstate[j], k);
+ for (int h = 0; h<nstate[j]; ++h){
+ Xwj(h, _) = Xj(h,_)*wi[h];
+ }
+ Matrix<> Bn = invpd(B0inv + ::t(Xj)*Xwj);
+ Matrix<> bn = Bn*gaxpy(B0inv, b0, -1*::t(Xj)*yj);
+ beta(j,_) = stream.rmvnorm(bn, Bn);
+ }
+
+ //////////////////////
+ // 3. 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);
+ }
+
+ //////////////////////
+ // 4. Sample tau
+ //////////////////////
+ Matrix <> TAUout = tau_comp_sampler(stream, *m, totcomp, Y, X, wr, mr, sr, beta, s);
+ tau = TAUout(_, 0);
+ component = TAUout(_, 1);
+
+ if (iter >= *burnin && ((iter % *thin)==0)){
+ Matrix<> 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(_, 0);
+ tau_store(count,_) = tau(_, 0);
+ component_store(count,_) = component(_, 0);
+
+ ++count;
+
+ }
+
+ if(*verbose > 0 && iter % *verbose == 0){
+ Rprintf("\n\n MCMCpoissonChange 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.5f", j+1, static_cast<double>(nstate[j]));
+ }
+ for (int i = 0; i<ns; ++i){
+ for (int j = 0; j<k; ++j){
+ Rprintf("\n beta(%i) in state %i is %10.5f", j+1, i+1, beta(i, j));
+ }
+ }
+ }
+
+ }// end MCMC loop
+
+
+ //////////////////////////////////////
+ if(*chib ==1){
+ //////////////////////////////////////
+
+ Matrix<> 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<> 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];
+ }
+
+ //////////////////////
+ // beta
+ //////////////////////
+ Matrix<> density_beta(nstore, ns);
+
+ for (int iter = 0; iter<nstore; ++iter){
+ int y_count = 0;
+ for (int t = 0; t<n ; ++t){
+ int yt = (int) Y[t];
+ y_count = y_count + yt + 1;
+ Matrix<> Yt_over_Sigma(yt + 1, 1);
+ Matrix<> Sigma_inv(yt + 1, 1);
+ for(int j = (y_count - yt - 1); j< y_count; ++j){
+ int jone = (int)component_store(iter, j) - 1 ;
+ Sigma_inv[j-(y_count-yt-1)] = 1/(sr[jone]*sr[jone]);
+ Yt_over_Sigma[j-(y_count-yt-1)] = (log(tau_store(iter, j))- mr[jone])*Sigma_inv[j-(y_count-yt-1)];
+ }
+ y_tilde[t] = sum(Yt_over_Sigma);
+ Sigma_inv_sum[t] = sum(Sigma_inv);
+ }
+
+ 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_store(iter, i) == (j+1)) {
+ nstate[j] = nstate[j] + 1;
+ }
+ }
+ beta_count = beta_count + nstate[j];
+ Matrix<> yj = y_tilde((beta_count - nstate[j]), 0, (beta_count - 1), 0);
+ Matrix<> Xj = X((beta_count - nstate[j]), 0, (beta_count - 1), k-1);
+ Matrix<> wi = Sigma_inv_sum((beta_count - nstate[j]), 0, (beta_count - 1), 0);
+ Matrix<> Xwj(nstate[j], k);
+ for (int h = 0; h<nstate[j]; ++h){
+ Xwj(h, _) = Xj(h,_)*wi[h];
+ }
+ Matrix<> Bn = invpd(B0inv + ::t(Xj)*Xwj);
+ Matrix<> bn = Bn*gaxpy(B0inv, b0, -1*::t(Xj)*yj);
+ density_beta(iter, j) = exp(lndmvn(::t(beta_st(j,_)), bn, Bn));
+ }
+ }
+
+ double pdf_beta = log(prod(meanc(density_beta)));
+
+ //////////////////////
+ // P
+ //////////////////////
+ Matrix<> density_P(nstore, ns);
+ for (int iter = 0; iter < nstore; ++iter){
+ Matrix <> Sout = poisson_reg_state_sampler(stream, *m, Y, X, beta_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)) {
+ 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){
+ int yt = (int) Y[t];
+ Matrix<> lambda = exp(X(t,_)*::t(beta_st));
+ 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_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);
+ }
+
+ const double loglike = sum(log(like));
+
+ //////////////////////
+ // prior
+ //////////////////////
+ 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)));
+ }
+
+ // compute marginal likelihood
+ const double logprior = sum(density_beta_prior) + sum(density_P_prior);
+ const double logmarglike = (loglike + logprior) - (pdf_beta + pdf_P);
+
+ logmarglikeholder[0] = logmarglike;
+
+ }
+
+ R_CheckUserInterrupt();
+
+ 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];
+ }
+
+}
+
+extern "C" {
+ void MCMCpoissonChange( 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 double *betastart,
+ const double *Pstart,
+ const double *taustart,
+ const double *componentstart,
+ const double *a,
+ const double *b,
+ const double *c0,
+ const double *d0,
+ const int* uselecuyer,
+ const int* seedarray,
+ const int* lecuyerstream,
+ const double *b0data,
+ const double *B0data,
+ const double *A0data,
+ double *logmarglikeholder,
+ const double *wrin,
+ const double *mrin,
+ const double *srin,
+ const int *chib){
+ if(*Xcol == 1){
+ MCMCPACK_PASSRNG2MODEL(MCMCpoissonChangepoint_impl,
+ betaout, Pout, psout, sout,
+ Ydata, Yrow, Ycol,
+ m, c0, d0,
+ burnin, mcmc, thin, verbose,
+ betastart, Pstart,
+ a, b, A0data,
+ logmarglikeholder, chib)
+ }
+ else{
+ MCMCPACK_PASSRNG2MODEL(MCMCpoissonRegChangepoint_impl,
+ betaout, Pout, psout, sout,
+ Ydata, Yrow, Ycol,
+ Xdata, Xrow, Xcol,
+ m, burnin, mcmc, thin, verbose,
+ betastart, Pstart,
+ taustart, componentstart,
+ a, b, b0data, B0data, A0data,
+ logmarglikeholder,
+ wrin, mrin, srin, chib);
+ }
+ } // end MCMC
+
+} // end extern "C"
+
+
+#endif
+
+
diff --git a/src/MCMCpoissonChangepoint.cc b/src/MCMCpoissonChangepoint.cc
deleted file mode 100644
index 17bc6dd..0000000
--- a/src/MCMCpoissonChangepoint.cc
+++ /dev/null
@@ -1,386 +0,0 @@
-//////////////////////////////////////////////////////////////////////////
-// 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) 2003-2007 Andrew D. Martin and Kevin M. Quinn
-// Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
-// and Jong Hee Park
-//////////////////////////////////////////////////////////////////////////
-
-#ifndef MCMCPOISSONCHANGEPOINT_CC
-#define MCMCPOISSONCHANGEPOINT_CC
-
-#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
--
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