[r-cran-mcmcpack] 20/90: Imported Upstream version 0.8-1
Andreas Tille
tille at debian.org
Fri Dec 16 09:07:35 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 f89446bd0d8bd4a25207f83887c166d43c2c5351
Author: Andreas Tille <tille at debian.org>
Date: Fri Dec 16 08:07:08 2016 +0100
Imported Upstream version 0.8-1
---
DESCRIPTION | 10 +-
HISTORY | 16 +++
INDEX | 64 ----------
NAMESPACE | 4 +
R/MCMCpoissonChangepoint.R | 272 ++++++++++++++++++++++++++++++++++++++++++
R/automate.R | 4 +-
R/btsutil.R | 262 ++++++++++++++++++++++++++++++++++++++++
R/zzz.R | 1 +
README | 64 ++++------
man/MCMCSVDreg.Rd | 23 ++--
man/MCMCdynamicEI.Rd | 14 ++-
man/MCMCfactanal.Rd | 2 +-
man/MCMChierEI.Rd | 9 +-
man/MCMCirt1d.Rd | 2 +-
man/MCMCirtKd.Rd | 10 +-
man/MCMCirtKdRob.Rd | 12 +-
man/MCMClogit.Rd | 6 +-
man/MCMCmetrop1R.Rd | 2 +-
man/MCMCmixfactanal.Rd | 8 +-
man/MCMCmnl.Rd | 2 +-
man/MCMCoprobit.Rd | 4 +-
man/MCMCordfactanal.Rd | 2 +-
man/MCMCpanel.Rd | 4 +-
man/MCMCpoisson.Rd | 2 +-
man/MCMCpoissonChangepoint.Rd | 137 +++++++++++++++++++++
man/MCMCprobit.Rd | 2 +-
man/MCMCregress.Rd | 2 +-
man/MCMCtobit.Rd | 7 +-
man/Nethvote.Rd | 2 +-
man/dtomog.Rd | 13 +-
man/plotPostChangepoint.Rd | 28 +++++
man/plotPostState.Rd | 30 +++++
man/procrust.Rd | 6 +-
man/readscythe.Rd | 8 +-
man/tomog.Rd | 5 +-
man/writescythe.Rd | 8 +-
36 files changed, 858 insertions(+), 189 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index 98d46c2..5a1e4b2 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,14 +1,14 @@
Package: MCMCpack
-Version: 0.7-4
-Date: 2006-9-29
+Version: 0.8-1
+Date: 2007-1-11
Title: Markov chain Monte Carlo (MCMC) Package
Author: Andrew D. Martin <admartin at wustl.edu>, and
Kevin M. Quinn <kevin_quinn at harvard.edu>
Maintainer: Andrew D. Martin <admartin at wustl.edu>
-Depends: R (>= 2.2.0), coda (>= 0.10-3), MASS
+Depends: R (>= 2.2.0), coda (>= 0.10-5), MASS, stats
Description: This package contains functions to perform Bayesian
inference using posterior simulation for a number of statistical
- models. All simulation is done in compiled C++ written in the Scythe
+ models. Most simulation is done in compiled C++ written in the Scythe
Statistical Library Version 1.0. All models return coda mcmc objects
that can then be summarized using the coda package. MCMCpack
also contains some useful utility functions, including some
@@ -17,4 +17,4 @@ Description: This package contains functions to perform Bayesian
sampling algorithm, and tools for visualization.
License: GPL version 2
URL: http://mcmcpack.wustl.edu
-Packaged: Thu Sep 28 21:16:37 2006; adm
+Packaged: Thu Jan 11 17:48:15 2007; adm
diff --git a/HISTORY b/HISTORY
index 5ab22f6..53637bd 100644
--- a/HISTORY
+++ b/HISTORY
@@ -2,9 +2,25 @@
// Changes and Bug Fixes
//
+0.7-4 to 0.8-1
+ * added MCMCPoissonChangepoint() model (authored by Jong Hee Park)
+ * added two plot methods for changepoint models:
+ plotPostChangepoint() and plotPostState() (authored by Jong Hee Park)
+ * package cleaned up prior to submission of piece to JStatSoft, including
+ major edit of documentation
+
0.7-3 to 0.7-4
* fixed minor bug in MCMCpoisson() that was causing the function to
not work on Windows machines.
+ * changed
+ test <- grep("\.nonconst$", c.names)
+ to
+ test <- grep("\\.nonconst$", c.names)
+ and
+ c.names <- sub("\.nonconst$", "", c.names)
+ to
+ c.names <- sub("\\.nonconst$", "", c.names)
+ in automate.R. [Thanks to Kurt Hornik for noticing this.]
0.7-2 to 0.7-3
* following posting by Radford Neal at:
diff --git a/INDEX b/INDEX
deleted file mode 100644
index 6f3e753..0000000
--- a/INDEX
+++ /dev/null
@@ -1,64 +0,0 @@
-BayesFactor Create an object of class BayesFactor from
- MCMCpack output
-Dirichlet The Dirichlet Distribution
-InvGamma The Inverse Gamma Distribution
-InvWishart The Inverse Wishart Distribution
-MCMCSVDreg Markov Chain Monte Carlo for SVD Regression
-MCMCdynamicEI Markov Chain Monte Carlo for Quinn's Dynamic
- Ecological Inference Model
-MCMCfactanal Markov Chain Monte Carlo for Normal Theory
- Factor Analysis Model
-MCMChierEI Markov Chain Monte Carlo for Wakefield's
- Hierarchial Ecological Inference Model
-MCMCirt1d Markov Chain Monte Carlo for One Dimensional
- Item Response Theory Model
-MCMCirtKd Markov Chain Monte Carlo for K-Dimensional Item
- Response Theory Model
-MCMCirtKdRob Markov Chain Monte Carlo for Robust
- K-Dimensional Item Response Theory Model
-MCMClogit Markov Chain Monte Carlo for Logistic
- Regression
-MCMCmetrop1R Metropolis Sampling from User-Written R
- function
-MCMCmixfactanal Markov Chain Monte Carlo for Mixed Data Factor
- Analysis Model
-MCMCmnl Markov Chain Monte Carlo for Multinomial
- Logistic Regression
-MCMCoprobit Markov Chain Monte Carlo for Ordered Probit
- Regression
-MCMCordfactanal Markov Chain Monte Carlo for Ordinal Data
- Factor Analysis Model
-MCMCpanel Markov Chain Monte Carlo for the General Linear
- Panel Model
-MCMCpoisson Markov Chain Monte Carlo for Poisson Regression
-MCMCprobit Markov Chain Monte Carlo for Probit Regression
-MCMCregress Markov Chain Monte Carlo for Gaussian Linear
- Regression
-MCMCtobit Markov Chain Monte Carlo for Gaussian Linear
- Regression with a Censored Dependent Variable
-MCbinomialbeta Monte Carlo Simulation from a Binomial
- Likelihood with a Beta Prior
-MCmultinomdirichlet Monte Carlo Simulation from a Multinomial
- Likelihood with a Dirichlet Prior
-MCnormalnormal Monte Carlo Simulation from a Normal Likelihood
- (with known variance) with a Normal Prior
-MCpoissongamma Monte Carlo Simulation from a Poisson
- Likelihood with a Gamma Prior
-Nethvote Dutch Voting Behavior in 1989
-NoncenHypergeom The Noncentral Hypergeometric Distribution
-PErisk Political Economic Risk Data from 62 Countries
- in 1987
-PostProbMod Calculate Posterior Probability of Model
-Senate 106th U.S. Senate Roll Call Vote Matrix
-SupremeCourt U.S. Supreme Court Vote Matrix
-Wishart The Wishart Distribution
-choicevar Handle Choice-Specific Covariates in
- Multinomial Choice Models
-dtomogplot Dynamic Tomography Plot
-procrustes Procrustes Transformation
-read.Scythe Read a Matrix from a File written by Scythe
-tomogplot Tomography Plot
-vech Extract Lower Triangular Elements from a
- Symmetric Matrix
-write.Scythe Write a Matrix to a File to be Read by Scythe
-xpnd Expand a Vector into a Symmetric Matrix
diff --git a/NAMESPACE b/NAMESPACE
index 8d7c36d..73a789b 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,6 +1,7 @@
useDynLib(MCMCpack)
import(coda)
import(MASS)
+import(stats)
export(
BayesFactor,
@@ -29,10 +30,13 @@ export(
MCMCordfactanal,
MCMCpanel,
MCMCpoisson,
+ MCMCpoissonChangepoint,
MCMCprobit,
MCMCregress,
MCMCSVDreg,
MCMCtobit,
+ plotPostState,
+ plotPostChangepoint,
PostProbMod,
procrustes,
rdirichlet,
diff --git a/R/MCMCpoissonChangepoint.R b/R/MCMCpoissonChangepoint.R
new file mode 100644
index 0000000..d4e9579
--- /dev/null
+++ b/R/MCMCpoissonChangepoint.R
@@ -0,0 +1,272 @@
+################################
+## Poisson Changepoint Model ##
+################################
+
+## This version does not allow covariates
+
+## Noninformative prior is not acceptable
+## in case of marginal likelihood computation
+
+## Originally written 12/23/2005 Jong Hee Park
+## Modified 02/13/2006 Jong Hee Park
+## Modified 10/21/2006 Jong Hee Park
+## Modified 12/13/2006 Jong Hee Park for JStatSoft article
+
+## some utility functions added
+## Jong Hee Park 07/14/2006
+
+## two plot functions added
+## Jong Hee Park 12/13/2006
+
+##############################################################
+## Andrew and Kevin, followings are a list of things to do
+##############################################################
+
+## 1. When you write a C code for MCMCpoissonChangepoint,
+## you need to export "prob.state" in addition to mcmc output.
+## This "prob.state" contains posterior probabilities of each state
+## and is essential to draw plots and posterior inference.
+
+## 2. I provide two special plot functions for changepoint models:
+## "plot.post.state" and "plot.post.changepoint."
+## These plot functions need to be changed based on the changes in
+## MCMCpoissonChangepoint() output.
+
+## 3. We need to change outputs of MCMCpoissonChangepoint make it accessible
+## by BayesFactor().
+
+## 4. All helper functions should be straightforward. "rdirichlet.cp" and
+## "ddirichlet.cp" might be redundant because they're simple a beta
+## distribution. In the future, I will replace these into beta functions
+## but for the time being, use these functions. Assigning zeros in the right
+## place is trickier than I thought. (I should have written these using
+## beta density in the first place...)
+
+
+
+#########################################################
+"MCMCpoissonChangepoint"<-
+#########################################################
+ function(data, m = 1, burnin = 1000, mcmc = 1000, thin = 1,
+ verbose = 0, seed = NA, c0, d0, a = NULL, b = NULL,
+ marginal.likelihood = c("none", "Chib95"), ...) {
+
+ ####################################################
+ ## Arguments
+ ## m : the number of changepoint
+ ## c0 and d0: gamma prior for lambda. NO DEFAULT values
+ ## a and b: 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.
+ ## ex) 300 years and 2 changepoints (3 states)
+ ## expected.duration <- 300/(2+1)
+ ## b <- 0.1
+ ## a <- b*expected.duration
+ ####################################################
+
+ # 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 <- data
+ n <- length(y)
+
+ ## prior
+ A0 <- trans.mat.prior(m=m, n=n, a=a, b=b)
+
+ ## get marginal likelihood argument
+ marginal.likelihood <- match.arg(marginal.likelihood)
+
+ ## storage matrices
+ lambda.store <- matrix(NA, mcmc/thin, m+1)
+ P.store <- matrix(NA, mcmc/thin, (m+1)^2)
+ ps.store <- matrix(0, n, m+1)
+ s1.store <- matrix(NA, mcmc/thin, n)
+ py <- rep(0, m+1)
+ pdf.P.store <- matrix(NA, mcmc, m+1)
+ lambda1 <- rep(NA, m+1)
+ P1 <- matrix(NA, m+1, m+1)
+
+ ## starting values
+ lambda0 <- runif(m+1)
+ P0 <- trans.mat.prior(m=m, n=n, a=0.9, b=0.1)
+
+ ####################################################################
+ ## MCMC iteration starts!
+ ####################################################################
+
+ for (iter in 1:totiter){
+
+ ####################################################################
+ ## Step 1: Sampling S
+ ####################################################################
+ state.out <- Poisson.state.sampler(m=m, y=y, lambda=lambda0, P=P0)
+ s1 <- state.out$s1
+ ps1<- state.out$ps1
+
+ ####################################################################
+ ## Step 2: Sampling lambda
+ ####################################################################
+ for (j in 1:(m+1)){
+ ej <- as.numeric(s1==j)
+ yj <- y[ej==1]
+ nj <- length(yj)
+ c1 <- sum(yj) + c0
+ d1 <- nj + d0
+ lambda1[j] <- rgamma(1, c1, d1)
+ }
+
+ ####################################################################
+ ## Step 3: Sampling P
+ ####################################################################
+ switch <- switchg(s1)
+ for (j in 1:(m+1)){
+ switch1 <- A0[j,] + switch[j,]
+ pj <- rdirichlet.cp(1, switch1)
+ P1[j,] <- pj
+ }
+
+ ## update
+ lambda0 <- lambda1
+ P0 <- P1
+
+ ####################################################################
+ ## end of one iteration
+ ####################################################################
+
+ ## store
+ if (iter > burnin && (iter %% thin == 0)) {
+ lambda.store[iter-burnin,] <- lambda1
+ P.store[iter-burnin,] <- as.vector(t(P1))
+ s1.store[iter-burnin,] <- s1
+ ps.store <- ps.store + ps1
+ }
+
+ ## report
+ if(verbose > 0 && iter%%verbose == 0){
+ cat("----------------------------------------------",'\n')
+ cat("iteration = ", iter, '\n')
+ cat("lambda = ", lambda1, '\n')
+ cat("Transition Matrix", '\n')
+ for(i in 1:(m+1))
+ cat(paste("", P1[i,]), fill=TRUE, labels=paste("{",i,"}:", sep=""), sep=",")
+ }
+
+ } ## end of MCMC sampling
+
+ ## marginal likelihood calculation if Chib == TRUE
+ if (marginal.likelihood == "Chib95"){
+
+ ############################################
+ ## Bayes Factor Calculation Starts ##
+ ############################################
+ lambda.st <- apply(lambda.store, 2, mean)
+ P.vec.st <- apply(P.store, 2, mean)
+ P.st <- t(matrix(P.vec.st, m+1, m+1))
+
+ #################################################################
+ ## 1. pdf.lambda
+ #################################################################
+ density.lambda <- matrix(NA, mcmc, m+1)
+ for (i in 1:mcmc){
+ for (j in 1:(m+1)){
+ # compute proposal density
+ ej <- as.numeric(s1.store[i,]==j)
+ yj <- y[ej==1]
+ nj <- length(yj)
+ c1 <- sum(yj) + c0
+ d1 <- nj + d0
+ density.lambda[i, j] <- dgamma(lambda.st[j], c1, d1)
+ }
+ }
+ pdf.lambda <- log(prod(apply(density.lambda, 2, mean)))
+
+ ######################################################################
+ ## 2. pdf.P
+ ######################################################################
+ for (g in 1:mcmc){
+ state.out <- Poisson.state.sampler(m=m, y=y, lambda=lambda.st, P=P0)
+ s1 <- state.out$s1
+ ps1 <- state.out$ps1
+ switch <- switchg(s1)
+
+ for (j in 1:(m+1)){
+ switch1 <- A0[j,] + switch[j,]
+ pj <- rdirichlet.cp(1, switch1)
+ P1[j,] <- pj
+ pdf.P.store[g,j]<- ddirichlet.cp(P.st[j,], switch1)
+ }
+ P0 <- P1
+ }
+ pdf.P <- log(prod(apply(pdf.P.store, 2, mean)))
+
+ ####################################################################
+ ## likelihood
+ ####################################################################
+ F <- matrix(NA, n, m+1)
+ like<- rep(NA, n)
+ pr1 <- c(1,rep(0,m))
+ for (t in 1:n){
+ py <- sapply(c(1:(m+1)), function(i){poisson.pdf(y[t], lambda.st[i])})
+ if(t==1) {pstyt1 = pr1}
+ else {pstyt1 <- F[t-1,]%*%P.st}
+ unnorm.pstyt <- pstyt1*py
+ pstyt <- unnorm.pstyt/sum(unnorm.pstyt)
+ F[t,] <- pstyt
+ like[t] <- sum(unnorm.pstyt)
+ }
+ loglik <- sum(log(like))
+
+ ####################################################################
+ ## prior ordinates
+ ####################################################################
+ nprior.lambda<- nprior.P <- rep(NA, m+1)
+ nprior.lambda<- sapply(c(1:(m+1)), function(i){dgamma(lambda.st[i], c0, d0, log=TRUE)})
+ nprior.P <- sapply(c(1:(m+1)), function(i){log(ddirichlet.cp(P.st[i,], A0[i,]))})
+
+ prior.lambda <- sum(nprior.lambda)
+ prior.P <- sum(nprior.P)
+
+ ############################################
+ ## Marginal Likelihood ##
+ ############################################
+ numerator <- loglik + prior.lambda + prior.P
+ denominator <- pdf.lambda + pdf.P
+ logmarglike <- numerator - denominator
+
+ ## print
+ if(verbose > 0){
+ cat("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n")
+ cat("Log Marginal Likelihood\n")
+ cat("-------------------------------------------------",'\n')
+ cat("log(marglike)= ", logmarglike, '\n')
+ cat("log(likelihood)= ", loglik, '\n')
+ cat("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n")
+ }
+ } ## end of marginal likelihood calculation
+
+ else {marginal <- loglik <- NULL}
+
+ ## return output
+ output <- as.mcmc(lambda.store)
+ varnames(output) <- paste("lambda.", 1:(m+1), sep = "")
+ attr(output, "title") <- "MCMCpoissonChangepoint Posterior Sample"
+ attr(output, "y") <- data
+ attr(output, "call") <- cl
+ attr(output, "logmarglike") <- logmarglike
+ attr(output, "loglik") <- loglik
+ attr(output, "prob.state") <- ps.store/mcmc
+ return(output)
+
+ }## end of MCMC function
+
diff --git a/R/automate.R b/R/automate.R
index 0d9ae9a..d2f397f 100644
--- a/R/automate.R
+++ b/R/automate.R
@@ -53,11 +53,11 @@
K <- length(objects)
c.names <- names(objects)
nonconst.indic <- rep(FALSE, K)
- test <- grep("\.nonconst$", c.names)
+ test <- grep("\\.nonconst$", c.names)
if(length(test)==0)
stop("something must be declared non-constant in auto.Scythe.call()\n")
nonconst.indic[test] <- TRUE
- c.names <- sub("\.nonconst$", "", c.names)
+ c.names <- sub("\\.nonconst$", "", c.names)
if(length(unique(c.names)) != K)
stop("non-unique nonconst names passed in auto.Scythe.call()\n")
diff --git a/R/btsutil.R b/R/btsutil.R
new file mode 100644
index 0000000..6babf17
--- /dev/null
+++ b/R/btsutil.R
@@ -0,0 +1,262 @@
+########## Utility Functions for Bayesian Times Series Models ##########
+
+## written and maintained by:
+## Jong Hee Park
+## Department of Political Science
+## University of Chicago
+## jhp at uchicago.edu
+
+## NOTE: only the plot functions are documented and exported in the
+## NAMESPACE.
+
+##############################################################
+## Helper functions for MCMCPoissonChangepoint()
+##############################################################
+
+ "switchg" <- function(s1){
+
+ ## "switchg" computes the frequency of elements in state vector
+ ## for one step ahead Markov process only
+ ## example: s1 <- c(1,1,1,2,2,2,3,4,4,4,4,5,5)
+ ## switchg(s1)
+ ## [,1] [,2] [,3] [,4] [,5]
+ ## [1,] 2 1 0 0 0
+ ## [2,] 0 2 1 0 0
+ ## [3,] 0 0 0 1 0
+ ## [4,] 0 0 0 3 1
+ ## [5,] 0 0 0 0 1
+
+ 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)
+ }
+
+ ## "trans.mat.prior" makes a transition matrix
+
+ ## a and b are beta prior for transition matrix
+ ## By default, the expected duration is computed and
+ ## corresponding a and b values are assigned. The expected
+ ## duration is the sample period divided by the number of states.
+ ## ex) 300 years and 2 changepoints (3 states)
+ ## expected.duration <- 300/(2+1)
+ ## b <- 0.1
+ ## a <- b*expected.duration
+
+ "trans.mat.prior" <- function(m, n, a=NULL, b=NULL){
+ if (!is.null(a)|!is.null(b)){
+ a <- a
+ 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)
+ }
+
+ "rdirichlet.cp" <- function(n, alpha){
+ ## "rdirichlet.cp" picks n random deviates from the Dirichlet function
+ ## with shape parameters alpha
+ ## Note that alpha can contain zero to deal with transition matrix rowwise.
+ ## It returns a matrix.
+
+ col <- length(alpha)
+ out <- matrix(NA, n, col)
+ out[,which(alpha==0)]<- 0
+ a <- alpha[alpha!=0]
+ l <- length(a);
+ x <- matrix(rgamma(l*n,a), ncol=l, byrow=TRUE);
+ sm <- x%*%rep(1,l);
+ dir.out <- x/as.vector(sm);
+ ## combine zero and nonzero parts prob
+ out[,which(alpha!=0)] <- dir.out
+ return(out)
+ }
+
+ "ddirichlet.cp" <- function (x, alpha){
+ ## "ddirichlet.cp" is a density function for the Dirichlet distribution
+ ## alpha and x can contain zero
+ ## but zero should be in the same place
+ ## x=c(0.9, 0.1, 0) and alpha=c(6,1,0)
+ ## It returns a matrix
+ col <- length(alpha)
+ out <- rep(NA, col)
+ out[which(alpha==0)] <- 0
+ a <- alpha[alpha!=0]
+ x1 <- x[x!=0]
+ dirichlet1 <- function(x1, a) {
+ logD <- sum(lgamma(a)) - lgamma(sum(a))
+ s <- sum((a - 1) * log(x1))
+ exp(sum(s) - logD)
+ }
+ if (!is.matrix(x1))
+ if (is.data.frame(x1))
+ x1 <- as.matrix(x1)
+ else x1 <- t(x1)
+ if (!is.matrix(a))
+ a <- matrix(a, ncol = length(a), nrow = nrow(x1),
+ byrow = TRUE)
+ if (any(dim(x1) != dim(a)))
+ stop("Mismatch between dimensions of x and alpha in
+ ddirichlet().\n")
+ pd <- vector(length = nrow(x1))
+ for (i in 1:nrow(x1)) pd[i] <- dirichlet1(x1[i, ], a[i,])
+ pd[apply(x1, 1, function(z) any(z < 0 | z > 1))] <- 0
+ pd[apply(x1, 1, function(z) all.equal(sum(z), 1) != TRUE)] <- 0
+ return(pd)
+ }
+
+ "poisson.pdf" <- function(y, lambda){
+ ## What this function does is same with "prod(dpois(y, lambda))"
+ n <- length(y)
+ log.like <- sum(sapply(c(1:n), function(i){-lambda + y[i]*log(lambda)- log(factorial(y[i]))}))
+ return(exp(log.like))
+ }
+
+ "Poisson.state.sampler" <- function(m, y, lambda, P){
+ ## "Poisson.state.sampler" samples state vector (s1) by forward and backward sampling
+ ## P is a (m+1 by m+1) transition matrix
+ ## F matrix that contains all the information of Pr(st|Yt)
+
+ ## Forward sampling: update F matrix
+ n <- length(y)
+ F <- matrix(NA, n, m+1) # storage for the Filtered probabilities
+ pr1 <- c(1,rep(0, m)) # initial probability Pr(s1=k|Y0, lambda)
+ for (t in 1:n){
+ py <- sapply(c(1:(m+1)), function(i){poisson.pdf(y[t], lambda[i])})
+ if(t==1) pstyt1 = pr1
+ else pstyt1 <- F[t-1,]%*%P
+ unnorm.pstyt <- pstyt1*py
+ pstyt <- unnorm.pstyt/sum(unnorm.pstyt) # Pr(st|Yt)
+ F[t,] <- pstyt
+ }
+
+ ## Backward recursions: sample s1 using F matrix and a transition matrix (P)
+ s1 <- matrix(NA, n, 1) ## holder for state variables
+ ps1 <- matrix(NA, n, m+1) ## holder for state probabilities
+ ps1[n,] <- F[n,] ## we know last elements of ps1 and s1
+ s1[n,1] <- m+1
+ t <- n-1
+
+ while (t>=1){
+ st1 <- s1[t+1]
+ unnorm.pstyn <- F[t,]*P[,st1]
+ ## normalize into a prob. density
+ pstyn <- unnorm.pstyn/sum(unnorm.pstyn)
+ if (st1==1) {s1[t]<-1}
+ ## If this is not the first period,
+ ## draw a state variable from a discrete prob. distribution("pstyn")
+ ## using the inverse CDF method.
+
+ ## What inverse CDF does is as follows
+ ## In backward recursion, the current state (t) is always at st1, which is one of (1...m+1)
+ ## The probability that transition happens at t-1, that is Pr(t-1=st1-1|t=st1) can be found
+ ## by using (dicrete case) inverse CDF method for the probability vector we have (pstyn).
+ ## draw u ~ U[0,1] then, if u<pstyn[st1-1], t-1=st1-1, otherwise t-1=st1
+ ## This is because Pr(x<t)=Pr(F^-1(u)<t)=Pr(u<F(t))=F(t)
+ ## Thus, when u ~ U[0,1], x=F^-1(u) ~ F
+
+ else {
+ pone <- pstyn[st1-1]
+ s1[t] <- ifelse (runif(1) < pone, st1-1, st1)}
+ ps1[t,] <- pstyn # probabilities pertaining to a certain state
+ t <- t-1
+ }
+ ## name and report outputs
+ out <- list(s1, ps1)
+ names(out)<-c("s1","ps1")
+ return(out)
+ }
+
+
+##############################################################
+## Plot functions for MCMCPoissonChangepoint()
+##############################################################
+
+ ## "plotPostState" draws a plot of posterior distribution of states
+
+ "plotPostState" <-
+ function(mcmcout, y, m, legend.control=NULL, cex=0.8, lwd=1.2){
+
+ ## To use legend.control, a user has to specify xy coordinate for legend.control.
+ ## ex) legend.control = c(x, y)
+
+ ## extract posterior state probability vector
+ out <- attr(mcmcout, "prob.state")
+
+ ## y must be a ts object. Otherwise, it is coerced as a ts object.
+ if (!is.ts(y)) y <- ts(y)
+ time.frame <- as.vector(time(y))
+
+ ## draw a plot for each state variable
+ plot(1, 0, xlim=range(time.frame), ylim=c(0,1), type="n", main="", xlab="Time", cex=cex, lwd=lwd,
+ ylab=expression(paste("Pr(", S[t], "= k |", Y[t], ")")))
+ for (i in 1:(m+1))
+ points(time.frame, out[,i], type="o", lty=i, lwd=lwd, col=i, cex=cex)
+
+ ## put a legend to indicate each line for the state
+ if (!is.null(legend.control)){
+ if (length(legend.control)!=2)
+ 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")
+ }
+ }
+
+ ## "plotPostChangepoint" draws a plot of posterior changepoint probability
+ "plotPostChangepoint" <-
+ function(mcmcout, y, m, xlab="Time", ylab=""){
+
+ ## extract posterior state probability vector
+ out <- attr(mcmcout, "prob.state")
+
+ ## y must be a ts object. Otherwise, it is coerced as a ts object.
+ if (!is.ts(y)) y <- ts(y)
+ time.frame <- as.vector(time(y))
+
+ ## divide space
+ if (m == 1){
+ pr.st <- c(0, diff(out[,(m+1)]))
+ plot(time.frame, pr.st, type="h", main="", xlab=xlab, ylab=ylab)
+ cp <- which(cumsum(pr.st) > 0.5)[1]-1
+ abline(v=cp+ time.frame[1], lty=3, col="red")
+ }
+ else {
+ par(mfrow=c(m,1));par(mar=c(2,4,1,1))
+ cp <- rep(NA, m) # holder for expected changepoints
+ for (i in 2:m){
+ pr.st <- c(0, diff(out[,i]))
+ pr.st <- ifelse(pr.st <0, 0, pr.st)
+ plot(time.frame, pr.st, type="h", main="", xlab=xlab, ylab=ylab)
+ cp[i-1] <- which(cumsum(pr.st)>0.5)[1]-1
+ abline(v=cp[i-1]+ time.frame[1], lty=3, col="red")
+
+ }
+ pr.st <- c(0, diff(out[,(m+1)]))
+ plot(time.frame, pr.st, type="h", main="", xlab=xlab, ylab=ylab)
+ cp[m] <- which(cumsum(pr.st)>0.5)[1]-1
+ abline(v=cp[m]+ time.frame[1], lty=3, col="red")
+ }
+
+ cp.means <- rep(NA, m+1)
+ cp.start <- c(1, cp + 1)
+ cp.end <- c(cp, length(y))
+ cat("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n")
+ cat("Expected changepoint(s) ", cp + time.frame[1], '\n')
+ for (i in 1:(m+1)) cp.means[i] <- mean(y[cp.start[i]:cp.end[i]])
+ cat("Local means for each regime are ", cp.means, '\n')
+ cat("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n")
+ }
diff --git a/R/zzz.R b/R/zzz.R
index 7666a6c..de34ba7 100644
--- a/R/zzz.R
+++ b/R/zzz.R
@@ -13,6 +13,7 @@
cat("## (Grants SES-0350646 and SES-0350613)\n##\n")
require(coda, quietly=TRUE)
require(MASS, quietly=TRUE)
+ require(stats, quietly=TRUE)
}
.onUnload <- function(libpath) {
diff --git a/README b/README
index 2589dc0..abcdb17 100644
--- a/README
+++ b/README
@@ -1,57 +1,39 @@
-//
-// MCMCpack, Version 0.6-6
-//
+/////////////////////
+// MCMCpack README //
+/////////////////////
-Release Date: November 14, 2005
-
-This package contains functions for posterior simulation for a number of
-statistical models. All simulation is done in compiled C++ written in
-the Scythe Statistical Library Version 1.0. All models return coda mcmc
-objects that can then be summarized using coda functions or the coda
-menu interface. The package also contains some useful utility
-functions, including some additional density functions and pseudo-random
-number generators for statistical distributions, a general purpose Metropolis sampling algorithm, and tools for visualization.
-
-//
// Authors
-//
Andrew D. Martin <admartin at wustl.edu>
Kevin M. Quinn <kevin_quinn at harvard.edu>
-//
// Compilation
-//
This package (along with Scythe) uses C++ and the Standard Template
-Library (STL). It requires use of the GCC compiler 3.0 or greater. It
-has been tested on GCC 3.1, 3.2, 3.3, 3.4, and 4.0 on Linux and MacOS X.
+Library (STL). We suggest using of the GCC compiler 4.0 or greater. The
+current package has been tested on GCC 3.1, 3.2, 3.3, 3.4, and 4.0 on Linux
+and MacOS X.
-We have also successfully compiled this package for Windows using
-a MinGW cross-compiler using binutils 2.13, gcc 3.2, mingw
-runtime 2.2, and w32api 2.0. This is essentially the same setup
-described by Yan and Rossini (2002), except for more recent
-versions of the compiler to successfully handle the C++ issues
-noted above. Many thanks to Dan Pemstein for helping with the
-cross-compilation, and to Kurt Hornik and Fritz Leisch for their
-help with debugging and service to the R community.
+We have also successfully compiled this package for Windows using a
+MinGW cross-compiler using the setup described by Yan and Rossini
+(2006). Many thanks to Dan Pemstein for helping with the
+cross-compilation, and to Kurt Hornik and Fritz Leisch for their help
+with debugging as well as their service to the R community.
-//
// Acknowledgments
-//
We gratefully acknowledge support from:
* National Science Foundation, Program in Methodology, Measurement, and
Statistics, Grants SES-0350646 and SES-0350613
-* Washington University, Department of Political Science and the
+* Washington University, Department of Political Science, the
Weidenbaum Center on the Economy, Government, and Public Policy
-http://wc.wustl.edu
+(http://wc.wustl.edu), and the Center for Empirical Research in the Law
+(http://cerl.wustl.edu)
* Harvard University, Department of Government and the
-Center for Basic Research in the Social Sciences
-http://cbrss.harvard.edu}
+Institute for Quantitative Social Sciences (http://iq.harvard.edu)
Neither the National Science Foundation, Washington University, or
Harvard University bear any responsibility for the content of this
@@ -62,15 +44,13 @@ you have any problems or questions.
--
Andrew D. Martin, Ph.D.
-Department of Political Science
-Washington University
-Campus Box 1063
-One Brookings Drive
-St. Louis, MO 63130
+Professor and CERL Director, School of Law
+Professor, Political Science, Arts & Sciences
+Washington University in St. Louis
+
(314) 935-5863 (Office)
-(314) 753-8377 (Cell)
-(314) 935-5856 (Fax)
+(314) 935-5150 (Fax)
-Office: Eliot Hall 326
+Office: Anheuser-Busch 470
Email: admartin at wustl.edu
-WWW: http://adm.wustl.edu
+WWW: http://adm.wustl.edu
diff --git a/man/MCMCSVDreg.Rd b/man/MCMCSVDreg.Rd
index 5f9a045..b181bc1 100644
--- a/man/MCMCSVDreg.Rd
+++ b/man/MCMCSVDreg.Rd
@@ -100,7 +100,7 @@ MCMCSVDreg(formula, data=parent.frame(), burnin = 1000, mcmc = 10000,
\deqn{y = X \beta + \varepsilon}{y = X beta + epsilon}
Where the errors are assumed to be iid Gaussian:
\deqn{\varepsilon_{i} \sim \mathcal{N}(0, \sigma^2)}{epsilon_i ~ N(0,
- sigma^2)}.
+ sigma^2)}
Let \eqn{N}{N} denote the number of rows of \eqn{X}{X} and \eqn{P}{P}
the number of columns of \eqn{X}{X}. Unlike the standard regression
@@ -133,15 +133,15 @@ MCMCSVDreg(formula, data=parent.frame(), burnin = 1000, mcmc = 10000,
}
\references{
- West, Mike; Josheph Nevins; Jeffrey Marks; Rainer Spang, and Harry
- Zuzan. 2000. "DNA MICROARRAY DATA ANALYSIS AND REGRESSION MODELING
- FOR GENETIC EXPRESSION PROFILING" Duke ISDS working paper.
+ Mike West, Josheph Nevins, Jeffrey Marks, Rainer Spang, and Harry
+ Zuzan. 2000. ``DNA Microarray Data Analysis and Regression Modeling
+ for Genetic Expression Profiling." Duke ISDS working paper.
- Gottardo, Raphael, and Adrian Raftery. 2004. "Markov chain Monte
- Carlo with mixtures of singular distributions". Statistics
- department, University of Washington, Technical report 470.
-
- Andrew D. Martin, Kevin M. Quinn, and Daniel Pemstein. 2004.
+ Gottardo, Raphael, and Adrian Raftery. 2004. ``Markov chain Monte
+ Carlo with mixtures of singular distributions.'' Statistics
+ Department, University of Washington, Technical Report 470.
+
+ Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007.
\emph{Scythe Statistical Library 1.0.} \url{http://scythe.wustl.edu}.
Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2002.
@@ -150,11 +150,6 @@ MCMCSVDreg(formula, data=parent.frame(), burnin = 1000, mcmc = 10000,
}
-\examples{
-\dontrun{
-}
-}
-
\keyword{models}
\seealso{\code{\link[coda]{plot.mcmc}},
diff --git a/man/MCMCdynamicEI.Rd b/man/MCMCdynamicEI.Rd
index d549a41..02c02d9 100644
--- a/man/MCMCdynamicEI.Rd
+++ b/man/MCMCdynamicEI.Rd
@@ -139,19 +139,21 @@ Dependence." In \emph{Ecological Inference: New Methodological
Strategies}. Gary King, Ori Rosen, and Martin A. Tanner (eds.). New
York: Cambridge University Press.
-Jonathan C. Wakefield. 2003. ``Ecological inference for 2x2 tables." Read
-before the Royal Statistical Society, on November 12th, 2003.
+
Radford Neal. 2003. ``Slice Sampling" (with discussion). \emph{Annals of
Statistics}, 31: 705-767.
- Andrew D. Martin, Kevin M. Quinn, and Daniel Pemstein. 2004.
- \emph{Scythe Statistical
- Library 1.0.} \url{http://scythe.wustl.edu}.
-
+ Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007.
+ \emph{Scythe Statistical Library 1.0.} \url{http://scythe.wustl.edu}.
+
Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2002.
\emph{Output Analysis and Diagnostics for MCMC (CODA)}.
\url{http://www-fis.iarc.fr/coda/}.
+
+Jonathan C. Wakefield. 2004. ``Ecological Inference for 2 x 2 Tables.'' \emph{Journal
+ of the Royal Statistical Society, Series A}. 167(3): 385�445.
+
}
\examples{
diff --git a/man/MCMCfactanal.Rd b/man/MCMCfactanal.Rd
index fa866f7..5a20c1d 100644
--- a/man/MCMCfactanal.Rd
+++ b/man/MCMCfactanal.Rd
@@ -157,7 +157,7 @@ MCMCfactanal(x, factors, lambda.constraints=list(),
}
\references{
- Andrew D. Martin, Kevin M. Quinn, and Daniel Pemstein. 2004.
+ Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007.
\emph{Scythe Statistical Library 1.0.} \url{http://scythe.wustl.edu}.
Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2002.
diff --git a/man/MCMChierEI.Rd b/man/MCMChierEI.Rd
index 34a8102..8728270 100644
--- a/man/MCMChierEI.Rd
+++ b/man/MCMChierEI.Rd
@@ -139,15 +139,14 @@ MCMChierEI(r0, r1, c0, c1, burnin=5000, mcmc=50000, thin=1,
}
\references{
-Jonathan C. Wakefield. 2003. ``Ecological inference for 2x2 tables." Read
-before the Royal Statistical Society, on November 12th, 2003.
+ Jonathan C. Wakefield. 2004. ``Ecological Inference for 2 x 2 Tables.'' \emph{Journal
+ of the Royal Statistical Society, Series A}. 167(3): 385�445.
Radford Neal. 2003. ``Slice Sampling" (with discussion). \emph{Annals of
Statistics}, 31: 705-767.
- Andrew D. Martin, Kevin M. Quinn, and Daniel Pemstein. 2004.
- \emph{Scythe Statistical
- Library 1.0.} \url{http://scythe.wustl.edu}.
+ Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007.
+ \emph{Scythe Statistical Library 1.0.} \url{http://scythe.wustl.edu}.
Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2002.
\emph{Output Analysis and Diagnostics for MCMC (CODA)}.
diff --git a/man/MCMCirt1d.Rd b/man/MCMCirt1d.Rd
index 8fbc827..96d94ae 100644
--- a/man/MCMCirt1d.Rd
+++ b/man/MCMCirt1d.Rd
@@ -184,7 +184,7 @@ MCMCirt1d(datamatrix, theta.constraints=list(), burnin = 1000,
Valen E. Johnson and James H. Albert. 1999. ``Ordinal Data Modeling."
Springer: New York.
- Andrew D. Martin, Kevin M. Quinn, and Daniel Pemstein. 2004.
+ Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007.
\emph{Scythe Statistical Library 1.0.} \url{http://scythe.wustl.edu}.
Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2002.
diff --git a/man/MCMCirtKd.Rd b/man/MCMCirtKd.Rd
index 3ff5294..8b27852 100644
--- a/man/MCMCirtKd.Rd
+++ b/man/MCMCirtKd.Rd
@@ -185,17 +185,17 @@ MCMCirtKd(datamatrix, dimensions, item.constraints=list(),
Curves Using Gibbs Sampling." \emph{Journal of Educational Statistics}.
17: 251-269.
- Joshua Clinton, Simon Jackman, and Douglas Rivers. 2000. ``The Statistical
- Analysis of Legislative Behavior: A Unified Approach." Paper presented at
- the Annual Meeting of the Political Methodology Society.
+ Joshua Clinton, Simon Jackman, and Douglas Rivers. 2004. ``The Statistical
+ Analysis of Roll Call Data." \emph{American Political Science Review}.
+ 98: 355-370.
Simon Jackman. 2001. ``Multidimensional Analysis of Roll Call Data
via Bayesian Simulation.'' \emph{Political Analysis.} 9: 227-241.
- Valen E. Johnson and James H. Albert. 1999. ``Ordinal Data Modeling."
+ Valen E. Johnson and James H. Albert. 1999. \emph{Ordinal Data Modeling}.
Springer: New York.
- Andrew D. Martin, Kevin M. Quinn, and Daniel Pemstein. 2004.
+ Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007.
\emph{Scythe Statistical Library 1.0.} \url{http://scythe.wustl.edu}.
Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2002.
diff --git a/man/MCMCirtKdRob.Rd b/man/MCMCirtKdRob.Rd
index af56250..81c457f 100644
--- a/man/MCMCirtKdRob.Rd
+++ b/man/MCMCirtKdRob.Rd
@@ -230,7 +230,7 @@ MCMCirtKdRob(datamatrix, dimensions, item.constraints=list(),
}
\details{
- \code{MCMCirtKd} simulates from the posterior using
+ \code{MCMCirtKdRob} simulates from the posterior using
the slice sampling algorithm of Neal (2003).
The simulation proper is done in
compiled C++ code to maximize efficiency. Please consult the
@@ -287,17 +287,17 @@ MCMCirtKdRob(datamatrix, dimensions, item.constraints=list(),
Kaplan. 2005. ``Practical Issues in Implementing and Understanding
Bayesian Ideal Point Estimation.'' \emph{Political Analysis}.
- Joshua Clinton, Simon Jackman, and Douglas Rivers. 2000. ``The Statistical
- Analysis of Legislative Behavior: A Unified Approach." Paper presented at
- the Annual Meeting of the Political Methodology Society.
+ Joshua Clinton, Simon Jackman, and Douglas Rivers. 2004. ``The Statistical
+ Analysis of Roll Call Data." \emph{American Political Science Review}.
+ 98: 355-370.
Simon Jackman. 2001. ``Multidimensional Analysis of Roll Call Data
via Bayesian Simulation.'' \emph{Political Analysis.} 9: 227-241.
- Valen E. Johnson and James H. Albert. 1999. ``Ordinal Data Modeling."
+ Valen E. Johnson and James H. Albert. 1999. \emph{Ordinal Data Modeling}.
Springer: New York.
- Andrew D. Martin, Kevin M. Quinn, and Daniel Pemstein. 2004.
+ Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007.
\emph{Scythe Statistical Library 1.0.} \url{http://scythe.wustl.edu}.
Radford Neal. 2003. ``Slice Sampling'' (with discussion). \emph{Annals of
diff --git a/man/MCMClogit.Rd b/man/MCMClogit.Rd
index 1636a37..81a21b7 100644
--- a/man/MCMClogit.Rd
+++ b/man/MCMClogit.Rd
@@ -115,7 +115,7 @@ MCMClogit(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
Additionally, arbitrary user-defined priors can be specified with the
\code{user.prior.density} argument.
- If the default multivariate normal prior is used, the Metropollis
+ If the default multivariate normal prior is used, the Metropolis
proposal distribution is centered at the current value of
\eqn{\beta}{beta} and has variance-covariance \eqn{V = T (B_0 +
C^{-1})^{-1} T }{V = T (B0 + C^{-1})^{-1} T}, where \eqn{T}{T} is a
@@ -124,7 +124,7 @@ MCMClogit(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
sample variance-covariance matrix of the MLEs. This last calculation
is done via an initial call to \code{glm}.
- If a user-defined prior is used, the Metropollis
+ If a user-defined prior is used, the Metropolis
proposal distribution is centered at the current value of
\eqn{\beta}{beta} and has variance-covariance \eqn{V = T
C T }{V = T C T}, where \eqn{T}{T} is a
@@ -136,7 +136,7 @@ MCMClogit(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
\references{
- Andrew D. Martin, Kevin M. Quinn, and Daniel Pemstein. 2004.
+ Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007.
\emph{Scythe Statistical Library 1.0.} \url{http://scythe.wustl.edu}.
Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2002.
diff --git a/man/MCMCmetrop1R.Rd b/man/MCMCmetrop1R.Rd
index b7f0d99..a933b3e 100644
--- a/man/MCMCmetrop1R.Rd
+++ b/man/MCMCmetrop1R.Rd
@@ -206,7 +206,7 @@ mode. This last calculation is done via an initial call to
B. Rubin. 2003. \emph{Bayesian Data Analysis}. 2nd Edition. Boca
Raton: Chapman & Hall/CRC.
- Andrew D. Martin, Kevin M. Quinn, and Daniel Pemstein. 2004.
+ Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007.
\emph{Scythe Statistical Library 1.0.} \url{http://scythe.wustl.edu}.
Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2002.
diff --git a/man/MCMCmixfactanal.Rd b/man/MCMCmixfactanal.Rd
index ed8aab3..a800392 100644
--- a/man/MCMCmixfactanal.Rd
+++ b/man/MCMCmixfactanal.Rd
@@ -216,8 +216,6 @@ MCMCmixfactanal(x, factors, lambda.constraints=list(),
}
\references{
- Kevin M. Quinn. 2004. ``Bayesian Factor Analysis for Mixed Ordinal and
- Continuous Responses.'' \emph{Political Analysis}. 12: 338-353.
M. K. Cowles. 1996. ``Accelerating Monte Carlo Markov Chain Convergence for
Cumulative-link Generalized Linear Models." \emph{Statistics and Computing.}
@@ -226,12 +224,16 @@ MCMCmixfactanal(x, factors, lambda.constraints=list(),
Valen E. Johnson and James H. Albert. 1999. ``Ordinal Data Modeling."
Springer: New York.
- Andrew D. Martin, Kevin M. Quinn, and Daniel Pemstein. 2004.
+ Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007.
\emph{Scythe Statistical Library 1.0.} \url{http://scythe.wustl.edu}.
Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2002.
\emph{Output Analysis and Diagnostics for MCMC (CODA)}.
\url{http://www-fis.iarc.fr/coda/}.
+
+ Kevin M. Quinn. 2004. ``Bayesian Factor Analysis for Mixed Ordinal and
+ Continuous Responses.'' \emph{Political Analysis}. 12: 338-353.
+
}
diff --git a/man/MCMCmnl.Rd b/man/MCMCmnl.Rd
index 7562f8e..accb635 100644
--- a/man/MCMCmnl.Rd
+++ b/man/MCMCmnl.Rd
@@ -147,7 +147,7 @@ that can be used to analyze the posterior sample.
\references{
- Andrew D. Martin, Kevin M. Quinn, and Daniel Pemstein. 2004.
+ Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007.
\emph{Scythe Statistical Library 1.0.} \url{http://scythe.wustl.edu}.
Radford Neal. 2003. ``Slice Sampling'' (with discussion). \emph{Annals of
diff --git a/man/MCMCoprobit.Rd b/man/MCMCoprobit.Rd
index 813c575..02ad1b6 100644
--- a/man/MCMCoprobit.Rd
+++ b/man/MCMCoprobit.Rd
@@ -109,10 +109,10 @@ that can be used to analyze the posterior sample.
Cumulative-link Generalized Linear Models." \emph{Statistics and Computing.}
6: 101-110.
- Valen E. Johnson and James H. Albert. 1999. ``Ordinal Data Modeling."
+ Valen E. Johnson and James H. Albert. 1999. \emph{Ordinal Data Modeling}.
Springer: New York.
- Andrew D. Martin, Kevin M. Quinn, and Daniel Pemstein. 2004.
+ Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007.
\emph{Scythe Statistical Library 1.0.} \url{http://scythe.wustl.edu}.
Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2002.
diff --git a/man/MCMCordfactanal.Rd b/man/MCMCordfactanal.Rd
index f8f7272..3f8f931 100644
--- a/man/MCMCordfactanal.Rd
+++ b/man/MCMCordfactanal.Rd
@@ -192,7 +192,7 @@ MCMCordfactanal(x, factors, lambda.constraints=list(),
Valen E. Johnson and James H. Albert. 1999. ``Ordinal Data Modeling."
Springer: New York.
- Andrew D. Martin, Kevin M. Quinn, and Daniel Pemstein. 2004.
+ Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007.
\emph{Scythe Statistical Library 1.0.} \url{http://scythe.wustl.edu}.
Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2002.
diff --git a/man/MCMCpanel.Rd b/man/MCMCpanel.Rd
index ae8331f..61dcecc 100644
--- a/man/MCMCpanel.Rd
+++ b/man/MCMCpanel.Rd
@@ -140,8 +140,8 @@ MCMCpanel(obs, Y, X, W, burnin = 1000, mcmc = 10000, thin = 5,
Hierarchical Longitudinal Models." \emph{Statistics and Computing.} 9:
17-26.
- Andrew D. Martin, Kevin M. Quinn, and Daniel Pemstein. 2004. \emph{Scythe
- Statistical Library 1.0.} \url{http://scythe.wustl.edu}.
+ Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007.
+ \emph{Scythe Statistical Library 1.0.} \url{http://scythe.wustl.edu}.
Andrew D. Martin and Kyle L. Saunders. 2002. ``Bayesian Inference for
Political Science Panel Data.'' Paper presented at the 2002 Annual Meeting
diff --git a/man/MCMCpoisson.Rd b/man/MCMCpoisson.Rd
index 40154c9..d41c4ed 100644
--- a/man/MCMCpoisson.Rd
+++ b/man/MCMCpoisson.Rd
@@ -110,7 +110,7 @@ MCMCpoisson(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
\references{
- Andrew D. Martin, Kevin M. Quinn, and Daniel Pemstein. 2004.
+ Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007.
\emph{Scythe Statistical Library 1.0.} \url{http://scythe.wustl.edu}.
Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2002.
diff --git a/man/MCMCpoissonChangepoint.Rd b/man/MCMCpoissonChangepoint.Rd
new file mode 100644
index 0000000..7a63728
--- /dev/null
+++ b/man/MCMCpoissonChangepoint.Rd
@@ -0,0 +1,137 @@
+\name{MCMCpoissonChangepoint}
+\alias{MCMCpoissonChangepoint}
+
+\title{Markov Chain Monte Carlo for a Poisson Multiple Changepoint Model}
+\description{
+ This function generates a sample from the posterior distribution
+ of a Poisson model with multiple changepoints. The function uses
+ the Markov chain Monte Carlo method of Chib (1998).
+ The user supplies data and priors, and
+ a sample from the posterior distribution is returned as an mcmc
+ object, which can be subsequently analyzed with functions
+ provided in the coda package.
+}
+
+\usage{MCMCpoissonChangepoint(data, m = 1, burnin = 1000, mcmc = 1000,
+ thin = 1, verbose = 0, seed = NA, c0, d0, a = NULL, b = NULL,
+ marginal.likelihood = c("none", "Chib95"), ...) }
+
+\arguments{
+ \item{data}{The data.}
+
+ \item{m}{The number of changepoints.}
+
+ \item{burnin}{The number of burn-in iterations for the sampler.}
+
+ \item{mcmc}{The number of MCMC iterations after burn-in.}
+
+ \item{thin}{The thinning interval used in the simulation. The number of
+ MCMC iterations must be divisible by this value.}
+
+ \item{verbose}{A switch which determines whether or not the progress of
+ the sampler is printed to the screen. If \code{verbose} is greater
+ than 0, the iteration number and the posterior density samples are printed to the screen every \code{verbose}th iteration.}
+
+ \item{seed}{The seed for the random number generator. If NA, current R
+ system seed is used.}
+
+ \item{c0}{\eqn{c_0}{c0} is the shape parameter for Gamma prior on \eqn{\lambda}{lambda}
+ (the mean).}
+
+ \item{d0}{\eqn{d_0}{d0} is the scale parameter for Gamma prior on \eqn{\lambda}{lambda}
+ (the mean).}
+
+ \item{a}{\eqn{a}{a} is the shape1 beta prior for transition probabilities. By default,
+ the expected duration is computed and corresponding a and b values are assigned. The expected
+ duration is the sample period divided by the number of states.}
+
+ \item{b}{\eqn{b}{b} is the shape2 beta prior for transition probabilities. By default,
+ the expected duration is computed and corresponding a and b values are assigned. The expected
+ duration is the sample period divided by the number of states.}
+
+ \item{marginal.likelihood}{How should the marginal likelihood be
+ calculated? Options are: \code{none} in which case the marginal
+ likelihood will not be calculated, and
+ \code{Chib95} in which case the method of Chib (1995) is used.}
+
+ \item{...}{further arguments to be passed}
+}
+
+\value{
+ An mcmc object that contains the posterior sample. This
+ object can be summarized by functions provided by the coda package. The object contains an attribute \code{prob.state} storage matrix that contains the probability of \eqn{state_i}{state_i} for each period, and the
+ log-likelihood of the model (\code{log.like}).
+}
+
+\details{
+ \code{MCMCpoissonChangepoint} simulates from the posterior distribution of
+ a Poisson model with multiple changepoints.
+
+ The model takes the following form:
+ \deqn{Y_t \sim \mathcal{P}oisson(\lambda_i),\;\; i = 1, \ldots, k}{Y_t ~ Poisson(lambda_i), i = 1,...,k.}
+ Where \eqn{k}{k} is the number of states.
+
+ We assume Gamma priors for \eqn{\lambda_{i}}{lambda_i} and Beta priors for transition probabilities:
+ \deqn{\lambda_i \sim \mathcal{G}amma(c_0, d_0)}{lambda_i ~ Gamma(c0, d0)}
+ \deqn{p_{ii} \sim \mathcal{B}eta{a}{b},\;\; i = 1, \ldots, k}{p_ii ~ Beta(a, b), i = 1,...,k.}
+ Where \eqn{k}{k} is the number of states.
+
+ Note that no default value is provided for Gamma priors.
+
+ The simulation in this model-fitting is performed in R.
+ }
+
+\author{Jong Hee Park, \email{jhp at uchicago.edu}, \url{http://home.uchicago.edu/~jhp/}.}
+
+\references{
+ Siddhartha Chib. 1995. "Marginal Likelihood from the Gibbs Output."
+ \emph{Journal of the American Statistical Association}. 90:
+ 1313-1321.
+
+ Siddhartha Chib. 1998. "Estimation and comparison of multiple change-point models."
+ \emph{Journal of Econometrics}. 86: 221-241.
+
+
+ Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2002.
+ \emph{Output Analysis and Diagnostics for MCMC (CODA)}.
+ \url{http://www-fis.iarc.fr/coda/}.
+}
+
+\examples{
+ \dontrun{
+ ## generate event count time series data with two breakpoints at 100 and 200
+ set.seed(1973)
+ n <- 300
+ true.lambda <- c(3, 2, 4)
+ y1 <- rpois(100, true.lambda[1])
+ y2 <- rpois(100, true.lambda[2])
+ y3 <- rpois(100, true.lambda[3])
+ y <- c(y1, y2, y3)
+
+ ## run the example
+ model1 <- MCMCpoissonChangepoint(y, m=1, burnin = 100, mcmc = 100,
+ thin=1, verbose=5, c0=3, d0=1, marginal.likelihood =
+ c("Chib95"))
+ model2 <- MCMCpoissonChangepoint(y, m=2, burnin = 100, mcmc = 100,
+ thin=1, verbose=5, c0=3, d0=1, marginal.likelihood =
+ c("Chib95"))
+ model3 <- MCMCpoissonChangepoint(y, m=3, burnin = 100, mcmc = 100,
+ thin=1, verbose=5, c0=3, d0=1, marginal.likelihood =
+ c("Chib95"))
+ model4 <- MCMCpoissonChangepoint(y, m=4, burnin = 100, mcmc = 100,
+ thin=1, verbose=5, c0=3, d0=1, marginal.likelihood =
+ c("Chib95"))
+
+ ## perform model comparison
+ BayesFactor(model1, model2, model3, model4)
+
+ ## Draw plots using the "right" model
+ plotPostState(model2, y, m=2, legend.control=NULL)
+ plotPostChangepoint(model2, y, m=2)
+ }
+}
+
+\keyword{models}
+
+\seealso{\code{\link{plotPostState}}, \code{\link{plotPostChangepoint}}}
+
diff --git a/man/MCMCprobit.Rd b/man/MCMCprobit.Rd
index 6986d6a..04552ca 100644
--- a/man/MCMCprobit.Rd
+++ b/man/MCMCprobit.Rd
@@ -110,7 +110,7 @@ MCMCprobit(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
Albert, J. H. and S. Chib. 1995. ``Bayesian Residual Analysis for
Binary Response Regression Models.'' \emph{Biometrika.} 82, 747-759.
- Andrew D. Martin, Kevin M. Quinn, and Daniel Pemstein. 2004.
+ Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007.
\emph{Scythe Statistical Library 1.0.} \url{http://scythe.wustl.edu}.
Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2002.
diff --git a/man/MCMCregress.Rd b/man/MCMCregress.Rd
index a563f26..fb694f2 100644
--- a/man/MCMCregress.Rd
+++ b/man/MCMCregress.Rd
@@ -125,7 +125,7 @@ MCMCregress(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
Robert E. Kass and Adrian E. Raftery. 1995. "Bayes Factors."
\emph{Journal of the American Statistical Association}. 90: 773-795.
- Andrew D. Martin, Kevin M. Quinn, and Daniel Pemstein. 2004.
+ Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007.
\emph{Scythe Statistical Library 1.0.} \url{http://scythe.wustl.edu}.
Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2002.
diff --git a/man/MCMCtobit.Rd b/man/MCMCtobit.Rd
index 3d21a5e..b0ceea2 100644
--- a/man/MCMCtobit.Rd
+++ b/man/MCMCtobit.Rd
@@ -124,7 +124,7 @@ MCMCtobit(formula, data = parent.frame(), below = 0, above = Inf,
}
\references{
- Andrew D. Martin, Kevin M. Quinn, and Daniel Pemstein. 2004.
+ Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007.
\emph{Scythe Statistical Library 1.0.} \url{http://scythe.wustl.edu}.
Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2002.
@@ -132,10 +132,11 @@ MCMCtobit(formula, data = parent.frame(), below = 0, above = Inf,
\url{http://www-fis.iarc.fr/coda/}.
James Tobin. 1958. ``Estimation of relationships for limited dependent
- variables." \emph{Econometrica.}, 26:24-36.
+ variables." \emph{Econometrica.} 26:24-36.
}
-\author{Ben Goodrich}
+\author{Ben Goodrich, \email{goodrich.ben at gmail.com},
+ \url{http://www.people.fas.harvard.edu/~goodrich/}}
\seealso{
\code{\link[coda]{plot.mcmc}},
diff --git a/man/Nethvote.Rd b/man/Nethvote.Rd
index 9df10be..88c4aa4 100644
--- a/man/Nethvote.Rd
+++ b/man/Nethvote.Rd
@@ -49,7 +49,7 @@ Dutch Voting Behavior in 1989.
Netherlands Central Bureau of Statistics, Amsterdam.
}
\references{
- K.M. Quinn and A.M. Martin. 2002. ``An Integrated Computational Model
+ Kevin M. Quinn and Andrew D. Martin. 2002. ``An Integrated Computational Model
of Multiparty Electoral Competition.'' \emph{Statistical Science}. 17:
405-419.
}
diff --git a/man/dtomog.Rd b/man/dtomog.Rd
index e9774a1..5454a27 100644
--- a/man/dtomog.Rd
+++ b/man/dtomog.Rd
@@ -74,13 +74,14 @@ dtomogplot(r0, r1, c0, c1, time.vec=NA, delay=0,
Gary King, 1997. \emph{A Solution to the Ecological Inference Problem}.
Princeton: Princeton University Press.
- Jonathan Wakefield. 2001. ``Ecological Inference for 2 x 2 Tables,''
- Center for Statistics and the Social Sciences Working Paper
- no. 12. University of Washington.
+ Jonathan C. Wakefield. 2004. ``Ecological Inference for 2 x 2 Tables.''
+ \emph{Journal of the Royal Statistical Society, Series A}. 167(3): 385�445.
+
+Kevin Quinn. 2004. ``Ecological Inference in the Presence of Temporal
+Dependence." In \emph{Ecological Inference: New Methodological
+Strategies}. Gary King, Ori Rosen, and Martin A. Tanner (eds.). New
+York: Cambridge University Press.
- Kevin M. Quinn. 2002. ``Ecological Inference in the Presence of
- Temporal Dependence.'' Paper prepared for Ecological Inference
- Conference, Harvard University, June 17-18, 2002.
}
\examples{
diff --git a/man/plotPostChangepoint.Rd b/man/plotPostChangepoint.Rd
new file mode 100644
index 0000000..9ba2ffd
--- /dev/null
+++ b/man/plotPostChangepoint.Rd
@@ -0,0 +1,28 @@
+\name{plotPostChangepoint}
+\alias{plotPostChangepoint}
+\title{Changepoint Location Plots}
+\description{Plot the posterior distribution over the location of the changepoints.}
+
+\usage{
+ plotPostChangepoint(mcmcout, y, m, xlab="Time", ylab="")
+}
+
+\arguments{
+
+\item{mcmcout}{The \code{mcmc} object containing the posterior density sample from a changepoint model. Note that this must have a \code{prob.state} attribute.}
+
+\item{y}{The data.}
+
+\item{m}{The number of changepoints.}
+
+\item{xlab}{Label for the x-axis.}
+
+\item{ylab}{Label for the y-axis.}
+
+}
+
+\author{Jong Hee Park, \email{jhp at uchicago.edu}, \url{http://home.uchicago.edu/~jhp/}.}
+
+\keyword{hplot}
+
+\seealso{\code{\link{MCMCpoissonChangepoint}}}
diff --git a/man/plotPostState.Rd b/man/plotPostState.Rd
new file mode 100644
index 0000000..5dc078d
--- /dev/null
+++ b/man/plotPostState.Rd
@@ -0,0 +1,30 @@
+\name{plotPostState}
+\alias{plotPostState}
+\title{Changepoint State Plot}
+\description{Plot the posterior probability that each time point is in each state.}
+
+\usage{
+ plotPostState(mcmcout, y, m, legend.control=NULL, cex=0.8, lwd=1.2)
+}
+
+\arguments{
+
+\item{mcmcout}{The \code{mcmc} object containing the posterior density sample from a changepoint model. Note that this must have a \code{prob.state} attribute.}
+
+\item{y}{The data.}
+
+\item{m}{The number of changepoints.}
+
+\item{legend.control}{Control the location of the legend. It is necessary to pass both the x and y locations; i.e., \code{c(x,y)}.}
+
+\item{cex}{Control point size.}
+
+\item{lwd}{Line width parameter.}
+
+}
+
+\author{Jong Hee Park, \email{jhp at uchicago.edu}, \url{http://home.uchicago.edu/~jhp/}.}
+
+\keyword{hplot}
+
+\seealso{\code{\link{MCMCpoissonChangepoint}}}
diff --git a/man/procrust.Rd b/man/procrust.Rd
index f222a82..a9143ee 100644
--- a/man/procrust.Rd
+++ b/man/procrust.Rd
@@ -30,13 +30,13 @@ procrustes(X, Xstar, translation=FALSE, dilation=FALSE)
}
\details{
- \code{R}, \code{tt}, and \code{s} are chosen so that
+ \code{R}, \code{tt}, and \code{s} are chosen so that:
- \deqn{s X R + 1 tt' \approx X^*}{s X R + 1 tt' approximately Xstar}.
+ \deqn{s X R + 1 tt' \approx X^*}{s X R + 1 tt' approximately Xstar}
\code{X.new} is given by:
- \deqn{X_{new} = s X R + 1 tt'}{X.new = s X R + 1 tt'}.
+ \deqn{X_{new} = s X R + 1 tt'}{X.new = s X R + 1 tt'}
}
diff --git a/man/readscythe.Rd b/man/readscythe.Rd
index 7d55b7d..c1a3c55 100644
--- a/man/readscythe.Rd
+++ b/man/readscythe.Rd
@@ -21,12 +21,14 @@
}
\examples{
- \dontrun{mymatrix <- read.Scythe("myfile.txt")}
+ \dontrun{
+ mymatrix <- read.Scythe("myfile.txt")
+ }
}
\references{
- Andrew D. Martin, Kevin M. Quinn, and Daniel Pemstein. 2004.
- \emph{Scythe Statistical Library 1.0.} \url{http://scythe.wustl.edu}.
+ Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007.
+ \emph{Scythe Statistical Library 1.0.} \url{http://scythe.wustl.edu}.
}
\keyword{file}
diff --git a/man/tomog.Rd b/man/tomog.Rd
index 08b621f..d3b1725 100644
--- a/man/tomog.Rd
+++ b/man/tomog.Rd
@@ -61,9 +61,8 @@ tomogplot(r0, r1, c0, c1, xlab="fraction of r0 in c0 (p0)",
Gary King, 1997. \emph{A Solution to the Ecological Inference Problem}.
Princeton: Princeton University Press.
- Jonathan Wakefield. 2001. ``Ecological Inference for 2 x 2 Tables,''
- Center for Statistics and the Social Sciences Working Paper
- no. 12. University of Washington.
+ Jonathan C. Wakefield. 2004. ``Ecological Inference for 2 x 2 Tables.''
+ \emph{Journal of the Royal Statistical Society, Series A}. 167(3): 385�445.
}
\examples{
diff --git a/man/writescythe.Rd b/man/writescythe.Rd
index bfdd788..86f1d6e 100644
--- a/man/writescythe.Rd
+++ b/man/writescythe.Rd
@@ -25,12 +25,14 @@
}
\examples{
- \dontrun{write.Scythe(mymatrix, "myfile.txt")}
+ \dontrun{
+ write.Scythe(mymatrix, "myfile.txt")
+ }
}
\references{
- Andrew D. Martin, Kevin M. Quinn, and Daniel Pemstein. 2004.
- \emph{Scythe Statistical Library 1.0.} \url{http://scythe.wustl.edu}.
+ Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007.
+ \emph{Scythe Statistical Library 1.0.} \url{http://scythe.wustl.edu}.
}
\keyword{file}
--
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