[r-cran-mcmcpack] 75/90: Imported Upstream version 1.3-6
Andreas Tille
tille at debian.org
Fri Dec 16 09:07:50 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 3856757c29b57285d6dcfeb05e3f8c7b76a6bb4d
Author: Andreas Tille <tille at debian.org>
Date: Fri Dec 16 08:07:25 2016 +0100
Imported Upstream version 1.3-6
---
DESCRIPTION | 25 +-
HISTORY | 465 -----------
MD5 | 117 ++-
NAMESPACE | 21 +-
R/BayesFactors.R | 11 +-
R/MCMCdynamicIRT1d-b.R | 2 +-
R/MCMCdynamicIRT1d.R | 2 +-
R/MCMCintervention.R | 231 ------
R/MCMCirt1d.R | 2 +-
R/MCMCirtHier1d.R | 2 +-
R/MCMCirtKdRob.R | 2 +-
R/MCMCmetrop1R.R | 15 +-
R/MCMCmixfactanal.R | 2 +-
R/MCMCordfactanal.R | 2 +-
R/MCMCprobit.R | 2 +-
R/testpanelSubjectBreak.R | 2 +-
README | 20 +-
cleanup | 2 +-
man/HMMpanelFE.Rd | 18 +-
man/HMMpanelRE.Rd | 6 +-
man/MCMCSVDreg.Rd | 8 +-
man/MCMCbinaryChange.Rd | 24 +-
man/MCMCdynamicEI.Rd | 13 +-
man/MCMCfactanal.Rd | 7 +-
man/MCMChierEI.Rd | 7 +-
man/MCMChlogit.Rd | 7 +-
man/MCMChpoisson.Rd | 7 +-
man/MCMChregress.Rd | 7 +-
man/MCMCintervention.Rd | 209 -----
man/MCMCirt1d.Rd | 7 +-
man/MCMCirtHier1d.Rd | 16 +-
man/MCMCirtKd.Rd | 7 +-
man/MCMCirtKdHet.Rd | 7 +-
man/MCMCirtKdRob.Rd | 11 +-
man/MCMClogit.Rd | 7 +-
man/MCMCmetrop1R.Rd | 7 +-
man/MCMCmixfactanal.Rd | 7 +-
man/MCMCmnl.Rd | 7 +-
man/MCMCoprobit.Rd | 7 +-
man/MCMCoprobitChange.Rd | 15 +-
man/MCMCordfactanal.Rd | 7 +-
man/MCMCpoisson.Rd | 8 +-
man/MCMCpoissonChange.Rd | 4 +-
man/MCMCprobit.Rd | 8 +-
man/MCMCprobitChange.Rd | 2 -
man/MCMCquantreg.Rd | 8 +-
man/MCMCregress.Rd | 10 +-
man/MCMCregressChange.Rd | 23 +-
man/MCMCtobit.Rd | 9 +-
man/PErisk.Rd | 4 +-
man/Rehnquist.Rd | 2 +-
man/SSVSquantreg.Rd | 10 +-
man/Senate.Rd | 2 +-
man/SupremeCourt.Rd | 2 +-
man/plotChangepoint.Rd | 2 -
man/plotIntervention.Rd | 37 -
man/plotState.Rd | 5 +-
man/testpanelGroupBreak.Rd | 6 +-
man/testpanelSubjectBreak.Rd | 11 +-
src/MCMCintervention.cc | 1664 ----------------------------------------
src/MCMCregressChange.cc | 4 +-
src/{Makevars => Makevars.win} | 0
62 files changed, 294 insertions(+), 2870 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index 1d94975..f2ba3fa 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,24 +1,25 @@
Package: MCMCpack
-Version: 1.3-3
-Date: 2013-5-1
-Title: Markov chain Monte Carlo (MCMC) Package
+Version: 1.3-6
+Date: 2016-4-16
+Title: Markov Chain Monte Carlo (MCMC) Package
Author: Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park
Maintainer: Jong Hee Park <jongheepark at snu.ac.kr>
Depends: R (>= 2.10.0), coda (>= 0.11-3), MASS, stats
-Description: This package contains functions to perform Bayesian
+Imports: graphics, grDevices, lattice, methods, utils, mcmc, quantreg
+Description: Contains functions to perform Bayesian
inference using posterior simulation for a number of
statistical models. Most simulation is done in compiled C++
written in the Scythe Statistical Library Version 1.0.3. All
models return coda mcmc objects that can then be summarized
- using the coda package. MCMCpack also contains some useful
- utility functions, including some additional density functions
- and pseudo-random number generators for statistical
- distributions, a general purpose Metropolis sampling algorithm,
- and tools for visualization.
+ using the coda package. Some useful
+ utility functions such as density functions,
+ pseudo-random number generators for statistical
+ distributions, a general purpose Metropolis sampling algorithm,
+ and tools for visualization are provided.
License: GPL-3
SystemRequirements: gcc (>= 4.0)
-URL: http://mcmcpack.wustl.edu
-Packaged: 2013-05-01 00:37:02 UTC; parkjonghee
+URL: http://mcmcpack.berkeley.edu
+Packaged: 2016-04-15 00:19:09 UTC; parkjonghee
NeedsCompilation: yes
Repository: CRAN
-Date/Publication: 2013-05-01 07:02:08
+Date/Publication: 2016-04-15 08:32:53
diff --git a/HISTORY b/HISTORY
deleted file mode 100644
index 056bb9c..0000000
--- a/HISTORY
+++ /dev/null
@@ -1,465 +0,0 @@
-//
-// Changes and Bug Fixes
-//
-1.3-2 to 1.3-3
- * fixed an error in MCMCpoissonChange: is.na(lambda.mu) and dropped offset argument
- (thanks to Brian Ripley)
-
-1.3-1 to 1.3-2
- * fixed two mistakes in marginal likelihood computation in MCMCregress
- (thanks to Sid Chib). In computing log_sigma2,
- - exp(digamma) -> digamma
- - averaging over tot_iter -> averaging over nstore
- * added MCMCintervention(), MCMCregressChange(), and plotIntervention()
- * sigma.mu and sigma.var are added in MCMCregress()
- * lambda.mu and lambda.sigma are added in MCMCpoissonChange()
- * fixed minor mistakes and change expressions for a recent version of g++
- (thanks to Jeffrey Oldham)
-
-1.2-4 to 1.3-1
- * fixed many C++ issues related to clang and Solaris (many thanks to
- Brian Ripley for his significant contributions)
-
-1.2-3 to 1.2-4
- * fixed a bug in HMMmultivariateGaussian.cc (comparison of unsigned and
- singed integers)
- * cleaned up g++ 4.7 compatibility issues (thanks to Dan Pemstein and
- Chris Lawrence)
-
-1.2-2 to 1.2-3
- * fixed extraneous warnings in hidden.R (thanks to Jeffrey Arnold)
- * added make.breaklist.R
- * fixed example codes in HMMpanelFE.Rd and testpanelSubjectBreak.Rd
- * fixed testpanelSubjectBreak() to handle missing values in subject.id
-
-1.2-1 to 1.2-2
- * fixed bug in dwish() (thanks to Taráz E. Buck)
- * fixed .onAttach() to no longer require packages and use packageStartupMessage()
- (thanks to Brian Ripley)
- * updated to most recent version of Scythe library
- * cleaned up issues with cout, cerr, abort(), exit(), and pthread.h
- (thanks to Brian Ripley)
-
-1.1-5 to 1.2-1
- * added HMMpanelFE()
- * added HMMpanelRE()
- * added testpanelGroupBreak()
- * added testpanelSubjectBreak()
- * added MCMCresidualBreakAnalysis()
- * deleted #undef DO, #undef DS, #undef SO, #undef SS
- due to a symbol clash with Solaris' headers (thanks to Brian Ripley)
-
-1.1-4 to 1.1-5
- * fixed an error in MCMCregress: loglike should be dropped
-
-1.1-3 to 1.1-4
- * fixed an error in MCMChregress: NQ has to be replaced by NP (line 223)
- * in MCMChregress, MCMChpoisson, MCMChlogit, R_CheckUserInterrupt(void) is replaced by R_CheckUserInterrupt()
- * fixed an error in the marginal likelihood computation
- in MCMCproit.cc and MCMCregress.cc (thanks to Sid Chib)
- * changed B0 as prior precision in MCMCprobitChange, MCMCoprobitChange and MCMCpoissonChange
- * the acceptance rate is reported if verbose > 0
- * added the reporting of marginal likelihood components if verbose > 0
-
-1.0-11 to 1.1-2
- * fixed an error in MCMCoprobit.R
- * changed the documentation of MCMCoprobit.
- * Added m=0 option to all changepoint models. Now all changepoint models allow users
- to set m=0 for marginal likelihood computation.
- * Fixed a bug in MCMCbinaryChange. I don't know why but log(dbeta) in Scythe often
- generates NaN results while lndbeta1() does not.
-
-1.0-10 to 1.0-11
- * updated CITATION file in inst/ to coincide with JStatSoft publication
- * added reference to JStatSoft in references
-
-1.0-9 to 1.0-10
- * added CITATION in inst/ to coincide with JStatSoft publication
- * compressed two large data files (Senate and Nethvote)
-
-1.0-8 to 1.0-9
- * added MCMCprobitChange()
- * added MCMCoprobitChange()
- * modified MCMCprobit and MCMCprobitres for marginal likelihood estimation
-
-1.0-7 to 1.0-8
- * fixed some NAMESPACE issues [thanks to Shawn Treier]
-
-1.0-6 to 1.0-7
- * slight modifications to internals of rinvgamma to make parameterization
- more clear [thanks to Daniel Runcie and Richard Morey]
- * fixed some NAMESPACE issues with SSVSquantreg and its methods [written by Craig
- Reed]
-
-1.0-5 to 1.0-6
- * added SSVSquantreg [written by Craig Reed]
- * added functions to handle SSVSquantreg output [written by Craig Reed]
- * revisions to the MCMCquantreg function [written by Craig Reed]
- * revisions to the MCMCquantreg function [written by Craig Reed]
- * fixed parameterization issue with rinvgamma [thanks to Julian Stander]
-
-1.0-4 to 1.0-5
- * added heteroskedastic IRT model [written by Ben Lauderdale]
- * minor changes to error status in hierarchical IRT [written by Mike Malecki]
-
-1.0-3 to 1.0-4
- * little change to parameterization of BQR [contributed by Craig Reed]
-
-1.0-2 to 1.0-3
- * fix to hierarchical IRT documentation [written by Mike Malecki]
- * added Bayesian quantile regression [contributed by Craig Reed]
-
-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
-
-0.9-5 to 0.9-6
- * fixed bug in MCMCmetrop1R.R [thanks to many for spotting memory issue]
- * fixed by in Lecuyer seeds [thanks to Eduardo Leoni]
-
-0.9-4 to 0.9-5
- * changed how missing values are handled when calculating agreement scores
- for starting values in factor analysis and IRT models
- * added warning about incomparable models to BayesFactor
- * fixed bug with MCMCmnl example [thanks to Duncan Murdoch]
- * fixed formula interface bug for all models [thanks to Duncan Murdoch]
- * added hierarchical IRT code [written by Mike Malecki]
-
-0.9-3 to 0.9-4
- * fixed issue with MCMCmetrop1R.R [thanks to Jim Albert]
- * fixed Linux compilation issue [thanks to Chris Lawrence]
-
-0.9-2 to 0.9-3
- * added gcc 4.0 check in configure.ac and SystemRequirements
- * added verbose option to plotChangepoint()
-
-0.9-1 to 0.9-2
- * fixed Makevars per Ripley's email [PKG_CXXFLAGS to PKG_CPPFLAGS]
- * fixed encoding on some documentation files [thanks to Kurt Hornik]
- * added other estimation algorithms for MCMCoprobit
- * added other estimation algorithms for MCMCmnl
- * added MCMCpoissonChangepoint with estimation in compiled C++
- * a variety of other minor fixes
-
-0.8-2 to 0.9-1
- * first release of JStatSoft version
- * full port to Scythe 1.0.2 (nearly all C++ has been radically changed)
- * MCMCpanel() deprecated
- * a number of minor fixes [thanks to anonymous reviewers]
-
-0.8-1 to 0.8-2
- * models with multivariate normal priors now check for symmetry and positive
- semi-definiteness of precision matrix
- * fixed bug in how models with marginal likelihood calculations check for
- prior propriety [thanks to Gary Rosner for spotting this]
- * fixed bug in how MCMCmetrop1R handled a singular Hessian
- [thanks to Piers Dunstan for spotting this]
-
-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:
- http://www.stat.columbia.edu/~cook/movabletype/archives/2006/03/ adaptive_metrop.html
- switched a < to a <= in the shrinkage procedures used in the various slice
- sampling implementations.
- * modified tomogplot() and dtomogplot() to handle situations in which
- r1[i] == 0. [thanks to David Hugh-Jones for making this suggestion].
- * allowed for more user control of the initial call to optim in
- MCMCmetrop1R [thanks to Luca La Rocca for this suggestion].
- * allowed users to pass the variance-covariance matrix of the Gaussian
- proposal directly without resorting to a call to optim
- [thanks to Luca La Rocca for this suggestion].
- * fixed minor bug in developer mode in automate.R [thanks to Ben Goodrich].
- * fixed a minor bug in the developer mode documentation template.
- * fixed minor bug in nonconst call in MCMCregress.
-
-
-
-0.7-1 to 0.7-2
- * added procrustes() to NAMESPACE so that it can be seen (function and
- documentation already there).
- * fixed the prior in the full conditional for theta_i in MCMCirtKdRob.
- Previously was uniform on the unit circle which was not consistent with
- the documentation and could cause problems given how the starting values
- were chosen.
- * fixed prior for lambda_j in MCMCirtKdRob (similar to point above)
- * fixed the the calculation of the full conditional for gamma[i] in
- MCMCSVDreg()
- * Fixed a bug in MCMCbetabinomial() [thanks to John Wood].
-
-0.6-6 to 0.7-1
- * Added robust k-dimensional IRT model MCMCirtKdRob().
- * Added SVD regression MCMCSVDreg().
- * Updated auto.Scythe.call() to have any number of non-constants args; this
- is necessary to return other things from the C++ code besides the
- posterior density sample (including log-marginal likelihoods,
- acceptance rates, etc.).
- * Updated form.mcmc.object() to add additional attributes to an
- mcmc object to hold other quantities of interest (such as
- log-marginal likelihoods, acceptance rates, data, etc.).
- * For linear regression model, added option to compute log-marginal
- likelihood via the method of Chib (1995) or via the Laplace Approximation.
- * For logistic regression model, added option to compute log-marginal
- likelihood via the Laplace Approximation.
- * For probit regression model, added option to compute log-marginal
- likelihood via the Laplace Approximation.
- * For Poisson regression model, added option to compute log-marginal
- likelihood via the Laplace Approximation.
- * Added methods to calculate Bayes factors and posterior probabilities
- of models and handle log-marginal likelihoods.
- * Added teaching model MCbinomialbeta().
- * Added teaching model MCpoissongamma().
- * Added teaching model MCnormalnormal().
- * Added teaching model MCmultinomdirichlet().
- * patched xpnd() function to provide more functionality. Thanks to
- Gregor Gorjan for the patches.
- * added the function procrustes() that performs a Procrustes
- transformation of a matrix.
- * made some minor changes to MCMCirt1d() to help conserve memory when
- dealing with large datasets (MORE OPTIMIZATION NEEDS TO BE DONE IN
- TERMS OF BOTH SPEED AND MEMORY)
-
-0.6-5 to 0.6-6
- * fixed the std::accumulate problem pointed out by Kurt Hornik. Thanks
- to Dan Pemstein for tracking the problem down and making the fix.
- * fixed up how the force.samp option works with MCMCmetrop1R-- previously
- if a diagonal element of the Hessian was 0 the correction wouldn't work.
- * changed how additional arguments are passed to user-defined function
- in MCMCmetrop1R. This breaks old code but provides a more standard
- interface.
- * when logfun==FALSE in MCMCmetrop1R the initial call to optim() now
- maximizes the log of the user-defined function
- * cleaned up the *.Rd files for the model fitting functions a bit
- (more still needs to be done)
- * MCMClogit now optionally accepts a user-defined prior density
-
-0.6-4 to 0.6-5
- * added a check so one can run MCMCirt1d without passing constraints
- * fixed the Senate dataset so there are no duplicate names in the
- member variables, and modified examples accordingly
- * slightly edited DESCRIPTION file to get some important keywords
- in the first sentence
- * compute dinvgamma() on the log scale so it does not blow up
- for large shape and scale parameters
- * changed how MCMCmetrop1R handles a non-negative definite Hessian
- from optim--sampling can now proceed if the user flips a switch
- * fixed minor error in how std.mean was passed to MCMCmixfactanal
- * fixed the error in the beta full conditional in MCMCpanel that resulted
- from a typo in the Chib and Carlin paper
- * changed factor.score.eigen.start() in hidden.R to use squared
- distances rather than distances so as to be compatible with
- Poole, 2005. _Spatial Models of Parliamentary Voting_
- * fixed a minor bug in how the verbose argument was handled in the C++
- code for MCMCmetrop1R and MCMCfactanal [thanks to Lingji Chen]
- * fixed a minor documentation bug in the MCMCirtKd model [thanks to
- Guillermo Rosas]
- * fixed the text echoed at start-up so it does not have to be updated
- year to year
-
-0.6-3 to 0.6-4
- * fixed bug with verbose in MCMCmetrop1R pointed out by a referee for
- Rnews
- * added a little clarification to the lecuyer.h file about the dual-
- licensing scheme Pierre L'Ecuyer agreed to on 7 August 2004.
-
-0.6-2 to 0.6-3
- * cleaned up the docs for MCMCirtKd
- * MCMCirtKd now fits a model with a difficulty parameter rather than
- a negative difficulty (easiness ???) parameter
- * MCMCirtKd and MCMCordfactanal now do a better job of dropping
- variables that have 0 variance.
- * fixed Senate examples for MCMCirt1d and MCMCirtKd. EVENTUALLY THE
- SENATE DATA SHOULD BE FIXED DIRECTLY SO THAT THERE ARE NO DUPLICATE
- NAMES IN THE MEMBER VARIABLE.
- * users can interrupt all estimation algorithms now with CTL-C
- * changed the behavior of the verbose switch. Now if verbose is greater
- than 0 output is printed every verboseth iteration.
- * changed the settings when mcmc objects are created so that
- "start=burnin+1" and "end=mcmc+burnin"
- * fixed documentation of MCMCmetrop1R to make clear how data are passed
- to the log-posterior function
- * fixed factor.score.start.check to not return an error when a vector of
- starting values is passed
- * improved the way factor.score.start.check enforces constraints
-
-0.6-1 to 0.6-2
- * fixed documentation for rinvgamma() and dinvgamma() so that it is clear
- that these functions take shape and scale parameters as arguments.
- * fixed documentation and R code for riwish(), diwish(), rwish(), and
- dwish() so that it is more clear how these distributions are parameterized.
- * fixed line endings in MCMCtobit.cc
- * changed how MCMCordfactanal and MCMCmixfactanal report MH acceptance
- rates. They now report separate rates for each manifest variable.
- * fixed a bug with starting values in the IRT and factor models. If no
- starting values were passed, those created in factor.score.start.check()
- did not include the constraints (hard or soft), so the check failed. Now
- the constructed starting values meet the constraints.
-
-0.5-2 to 0.6-1
- * added a tobit model for a linear model with censoring in MCMCtobit()
- * added multinomial logit model in MCMCmnl()
- * added vote data from the Netherlands to illustrate MCMCmnl()
- * added choicevar() function to specify choice-specific variables in
- multinomial choice models
- * fixed some data-handling issues in MCMCmixfactanal() and
- MCMCordfactanal() [thanks to Ben Goodrich for isolating these and
- providing patches]
- * fixed the Metropolis-Hastings step for the Cowles algorithm for
- cutpoints in MCMCoprobit(), MCMCordfactanal(), and MCMCmixfactanal()
- [thanks to Alexander Raach for isolating the problem and providing
- a patch]
- * removed gcc specific compilation flags from Makevars.in (per the
- request of Brian Ripley and Kurt Hornik)
- * added informative message for templates created in auto.Scythe.call()
- * modified vector.tune() hidden function
- * fixed Neal's shrinkage procedure in MCMChierEI.cc
- * a number of editorial fixes, updating for new calendar year, etc.
- * removed functions acosh, asinh, atanh, and expm1 from smath.h and
- smath.cc so cross-compilation will work
-
-//
-// Old Changes and Bug Fixes
-//
-
-MCMCpack 0.5-1 was a major revision of MCMCpack. The entire package was
-been essentially re-written using the new development environment (documented
-in the MCMCpack specification) and the new Scythe Statistical Library 1.0.
-This following list summarizes major changes, but is by no means
-exhaustive.
-
-0.5-1 to 0.5-2
- * C++ code for truncated normal draws optimized for speed
- * with the permission of Pierre L'Ecuyer licensing of RngStream code
- changed to a dual license setup that is GPL compatible. Thanks to Chris
- Lawrence for bringing the licensing issues to our attention and
- drafting a new licensing statement and to Pierre L'Ecuyer for
- agreeing to use the new licensing statement for his RngStream
- code.
- * Fixed serious bugs in MCMChierEI() and MCMCdynamicEI()
- * Implemented a new sampling scheme based on slice sampling for
- MCMChierEI() and MCMCdynamicEI().
- * Removed MCMCbaselineEI()
- * Added delay argument to dtomogplog()
-
-0.4-8 to 0.5-1
- * NAMESPACE implemented
- * hidden functions are now available to aid in development (see hidden.R)
- * a function is available to automate the C++ call and generate
- template C++ code for estimation (see automate.R)
- * all model functions have been updated to the new specification, and
- most use hidden functions and automate
- * added a general purpose Metropolis sampler that allows the user
- to sample from an arbitrary (log)-density.
- * C++ code now using Scythe 1.0 (now using the unedited Scythe codebase
- through IFDEFs)
- * support for arbitrary random number generators, including
- the L'Ecuyer RNG for parallel computation (the RNG helper functions
- are available in MCMCrng.cc)
- * many full conditional distributions are available in MCMfcnds.cc
- * documentation for density functions and RNGs have been made "R-like"
- * fixed some spelling errors and misnomers in the documentation
- * all documentation updated to reflect changes
- * MCMCirt1d() has a new interface with new types of constraints--
- sampling for this model is also now much faster.
-
-0.4-8 to 0.4-9
- * Fixed a minor Scythe issue to fix error found by gcc 3.4.
-
-0.4-7 to 0.4-8
- * Repaired Scythe_Simulate.*, for which an outdated version was included
- in the last release.
-
-0.4-6 to 0.4-7
- * Fixed some Scythe bugs, including a problem with memory allocation for
- matrix multiplication. See http://sourceforge.net/projects/scythestat/
- for the latest version of Scythe, which is now distributed with MCMCpack.
- The Scythe code differs slightly in the paths ../include and ../src are
- replace with the current path, and in pnorm2 the isnan() function is
- commented out to allow for cross-compilation.
- * Rolled out http://mcmcpack.wustl.edu website.
- * Mixed response factor code.
- * Fixed factanal.
- * Fixed irtKd.
-
-0.4-5 to 0.4-6
- * Fixed a bug in rnoncenhypergeom() [thanks to Tom LaFramboise]
- * Patched Scythe0.3 to fix an error in inv() [thanks to Donour Sizemore].
- Note that this function is not called in MCMCpack, so was causing no
- explicit errors.
-
-0.4-3 to 0.4-5
- * Fixed a bug in xpnd() [thanks to Michael Man]
- * Fixed some inconsistencies in documentation [thanks to Kurt Hornik]
-
-0.4-2 to 0.4-3
- * Fixed bug in Scythe truncated Normal generators (which had been
- fixed before but sneaked into the last release) -- this fixes
- a problem with MCMCirt1d
- * Cleaned up MCMCbaselineDA.cc (eliminated unused arguments)
- * Cleaned up MCMCbaseline.R (tuning argument)
- * Set seed in MCMClogit.cc fixed
- * Set seed in MCMCpoisson.cc fixed
- * Fixed all examples such that they work out of the box
-
-0.4-1 to 0.4-2
- * Optimized some of the Scythe 0.4 code, which provides faster
- computation for most models.
- * Corrected a permissions problem on cleanup [thanks to Kurt Hornik]
- * Added explicit licensing information and a text echo when loading
- MCMCpack.
-
-0.3-11 to 0.4-1
- * Ported to Scythe Version 0.4 (which will soon be publicly available)
- * Cleaned up the codebase and documentation (changes will
- soon be part of the specification)
- * Added vech() and xpnd() utility functions
- * Included data file of 106th Senate roll call votes for the
- MCMCirt1d() and MCMCirtKd() models
- * Added Dirichlet, Noncentral Hypergeometric, and Inverse Gamma
- generators and densities [with contributions from Kevin Rompala]
- * Added read.Scythe() function to read matrices written by Scythe
- [contributed by Kevin Rompala]
- * Added helper functions to make coding easier [contributed by
- Kevin Rompala]
- * Added three models: a K-dimensional item response theory
- model (MCMCirtKd), a linear factor model (MCMCfactanal), and
- an ordinal item response theory model (MCMCordfactanal)
- * Added a pre-processor command to handle ininf() compilation
- issues on SGI [thanks to Dave Henderson]
- * All MCMC* functions now only allow starting values for the
- first simulated block of parameters and use check.parameters()
- function.
- * Range checking is turned off in the compiled C++ code, yielding
- significant speed gains for most models.
-
-0.3-10 to 0.3-11
- * Fixed a bug in MCMCpoisson() re: non-negative counts
- * Included a data file of Supreme Court votes for the
- MCMCirt1d() model [thanks to Simon Jackman for the suggestion]
- * Fixed memory leak caused by Scythe_Matrix.cc [thanks to Dan Pemstein]
diff --git a/MD5 b/MD5
index d334de5..3fd574d 100644
--- a/MD5
+++ b/MD5
@@ -1,36 +1,34 @@
-6470104b1e1bf32d945b869f9811a91e *DESCRIPTION
-19a6a667aa62fb26c81531564e05298c *HISTORY
-b1382cb26ee38f7a860c379826362e50 *NAMESPACE
-1eaf24e8d468ac0b624a9aa90cded306 *R/BayesFactors.R
+f93e438a313020ce15e634db8c557b9d *DESCRIPTION
+d332312be2ebe4eecbb2a8d76647127a *NAMESPACE
+21f17b9451d7ec3760fe3fe69336fe0c *R/BayesFactors.R
9c52f72ee21a71bfbfd312cb56233125 *R/HMMpanelFE.R
ca0b3fff125085e13069fae3f237e5ed *R/HMMpanelRE.R
626c7d968d4f56fe081b595cd59d4976 *R/MCMCSVDreg.R
c93b534213a07cec26b71eaa48f02295 *R/MCMCbinaryChange.R
2736b2f4759f92e4e1d49d92fd0811e8 *R/MCMCdynamicEI.R
-87107048db50664880297536b41ce8e7 *R/MCMCdynamicIRT1d-b.R
-b0810a4be00979118db6e89ea9f2be46 *R/MCMCdynamicIRT1d.R
+fd4d6882374299e9905ee3d31fc5e6ba *R/MCMCdynamicIRT1d-b.R
+40579946831f6a692aa34e59d054f1ff *R/MCMCdynamicIRT1d.R
f99ecb9f321a578b1cdfa9e5b6e163c6 *R/MCMCfactanal.R
dcd83d013877c39ace4ca33c4357c779 *R/MCMChierBetaBinom.R
065bd9bb18b1bc5d6011137361dc3412 *R/MCMChierEI.R
4b95a74aa87e670a88106d4661b6e833 *R/MCMChlogit.R
e2b64d567aa2e22410251a753f0b55c6 *R/MCMChpoisson.R
0a759c18d9f3b280dbe6090c92e20e35 *R/MCMChregress.R
-e4ece8084610294f593ca5fc632dd432 *R/MCMCintervention.R
-a4db065464be685b498d341ecc3c2109 *R/MCMCirt1d.R
-78b6d6361c4e9b3056226d66c50a3479 *R/MCMCirtHier1d.R
+7661ca467359c2c29e74cbf4ab65224c *R/MCMCirt1d.R
+08f85f3cfecd3e35b749077954b5bac7 *R/MCMCirtHier1d.R
4b433d0cace650f7304ff9a94e69ad4a *R/MCMCirtKd.R
7582e3ad855f5ece89f1dfbcf4a2a856 *R/MCMCirtKdHet.R
-0293dbca0de9cb96b203f1739cda2c02 *R/MCMCirtKdRob.R
+e84797c832aa60879c1de72a407d6342 *R/MCMCirtKdRob.R
ba21d6d4a4dcaaa56b6784b749bcbc07 *R/MCMClogit.R
-7a83fa7762a8e777a5cacc21c14f2498 *R/MCMCmetrop1R.R
-33f53d5be86c9d714dd7b68a27285c2f *R/MCMCmixfactanal.R
+80c184ceb693fdf962f999ca526247d7 *R/MCMCmetrop1R.R
+06b616772aa5a050a7826ff9986a446f *R/MCMCmixfactanal.R
a6508dd38804705e750995ca27070cde *R/MCMCmnl.R
4b7cb2b7fa8eb3d0435c55dfdcbed40e *R/MCMCoprobit.R
6b63259803b2890eae6d5778592ed7dd *R/MCMCoprobitChange.R
-be0928d3cd74e5485b1e26a9f0a15c59 *R/MCMCordfactanal.R
+14e27069ba0fcfa26f3888f9466013e9 *R/MCMCordfactanal.R
d438fc9dabd9993da5d4a45a3df68c39 *R/MCMCpoisson.R
cbf40bf99270db5f9dd8f4dea8f5274d *R/MCMCpoissonChange.R
-b98bdefca76d3316c664bcb2ba2bcc38 *R/MCMCprobit.R
+0fb7f4befd7db6df2f852b0b6e472237 *R/MCMCprobit.R
31fdcb514ab53ba258a581536fcbda60 *R/MCMCprobitChange.R
e27dcf535f23c25018bb778b61640f22 *R/MCMCquantreg.R
fe5ffab10fce83348751c3bd7abaa83a *R/MCMCregress.R
@@ -49,12 +47,12 @@ fdfcfe9909cadaf15c47bbb604be2257 *R/automate.R
f799711e5d3fb7140cbc0995360339b3 *R/procrust.R
2e00b2b1593d22b5d1c579e92b294106 *R/scythe.R
4de1106e38b8ca3b12396ec0d84bba4c *R/testpanelGroupBreak.R
-0907ed9bbd42212d50c6208dce57d640 *R/testpanelSubjectBreak.R
+01de088e7f9a83c20fa0a82897cbf440 *R/testpanelSubjectBreak.R
c67ac5238de55cdc8c49d017e0226ed9 *R/tomog.R
6377da3453982f578df4a886f673122c *R/utility.R
7e2021def8a88edb41fe78d0a0ae642c *R/zzz.R
-7bf94529cc64482070ef07fa6805012c *README
-46e955baeb55718a380a11fb34d2ea67 *cleanup
+a49acfed057fd07c87e294015890c211 *README
+31f24d1e199bde63095dcfa4556e848b *cleanup
ce8d6bc2b702996ba91799bc33b82840 *configure
21e1a006fb388ed39daf16e6d2ec9549 *configure.ac
e17202d22ce5ab5ab7b01e3b248a4008 *data/Nethvote.rda
@@ -64,52 +62,51 @@ accb7b5b46207a02e49de2702a6faff4 *data/Senate.rda
ce79b1b47b561cb44dd4f93589d3d755 *data/SupremeCourt.rda
a5cfcf73e21c03eeaf4f3baa5c492c14 *inst/CITATION
2b42855343b5f0c51c7e9411aef57509 *man/BayesFactor.Rd
-c3e38d2c8a1adee2ff096056416ca6ec *man/HMMpanelFE.Rd
-2e322e903df0d9c20513670710537995 *man/HMMpanelRE.Rd
-293ff17f852d0e60dd1aac1f419353d2 *man/MCMCSVDreg.Rd
-9af751257fcacc3a57543e4ae6a5e6cc *man/MCMCbinaryChange.Rd
-5078ab3b9d86c7165909791acf1fc886 *man/MCMCdynamicEI.Rd
+873c9a9c8f8e4808c417f14a857a562b *man/HMMpanelFE.Rd
+c199d84ef2c1edef1858a3b25b8189e0 *man/HMMpanelRE.Rd
+1fbc51fe87be958a30d1ae52085b9a6b *man/MCMCSVDreg.Rd
+e253ed3b10b78fa171e7a987fba251bb *man/MCMCbinaryChange.Rd
+b7f1245bf640f205c7927ac8b4dd5e1b *man/MCMCdynamicEI.Rd
102a94ced183cf023da7696f948bb82f *man/MCMCdynamicIRT1d.Rd
-29e89e27c6b11794c5dc6f50618375e0 *man/MCMCfactanal.Rd
-5dabac071a3722205135ae2cbae0b3d9 *man/MCMChierEI.Rd
-00a29263d43f0e2006de60d4800f5f7f *man/MCMChlogit.Rd
-b894aa07a0e68f9ffdeb3b09bab0aee2 *man/MCMChpoisson.Rd
-8d9c22e8d965d60581ec837803790c80 *man/MCMChregress.Rd
-9afab254ca5f275f437f7b22ab798b5a *man/MCMCintervention.Rd
-52f2cb9bc1cde5e38d35aa5d9858be4a *man/MCMCirt1d.Rd
-336682b062ff43261f0e75a4b51a34d8 *man/MCMCirtHier1d.Rd
-0d6b55f1f6062483b894fd4f4e9cecfd *man/MCMCirtKd.Rd
-19142a1e2a1a9217001a2d19b377b4ce *man/MCMCirtKdHet.Rd
-1bcb750570641123acbbc149bdad93af *man/MCMCirtKdRob.Rd
-1ee97a34364e07ec83542d78fce4d823 *man/MCMClogit.Rd
-8c4a033b4e5fca56699e9cbd96a56222 *man/MCMCmetrop1R.Rd
-751e56edabc12da2d8246f7c34cd7613 *man/MCMCmixfactanal.Rd
-b7dfcf274126008379b0ee903f134b25 *man/MCMCmnl.Rd
-cf07e2c9f307215c8fe841f5381b62f8 *man/MCMCoprobit.Rd
-1556620313c5a4ca175fbe0b68f84161 *man/MCMCoprobitChange.Rd
-a491294e97dd0b7aa5d0de525542ea8f *man/MCMCordfactanal.Rd
-03f2e93e6876a4f0e87b3469bfe38d2f *man/MCMCpoisson.Rd
-d754c71781dbdbcdfc4fecd6ea807056 *man/MCMCpoissonChange.Rd
-1917e296040b80c92525e5cb8ad02e71 *man/MCMCprobit.Rd
-c3631c612ccc2d8975a7e164c8506331 *man/MCMCprobitChange.Rd
-2f9caf8983e405f912eb5b25f29784eb *man/MCMCquantreg.Rd
-94f027a39b5a875cf5c54757a093f144 *man/MCMCregress.Rd
-65de16850240fd21bfe6e62c63fc5c10 *man/MCMCregressChange.Rd
+c217fc089f74217445d253851fc1914b *man/MCMCfactanal.Rd
+a2caad5ef711c10d14cdaf78f22eedc1 *man/MCMChierEI.Rd
+ff610bd17282d6410a06b39a6aef4894 *man/MCMChlogit.Rd
+8e4ae1e291d954cbe901c973b9559978 *man/MCMChpoisson.Rd
+cd31668e1de44486d2b4239d7ce0aca6 *man/MCMChregress.Rd
+9f93a019fac0c29b9509a3b963f52e5a *man/MCMCirt1d.Rd
+25b810b4cf59c1a11a49029ca916eb7c *man/MCMCirtHier1d.Rd
+a06d0a36b4288dac99ff3af14699a26b *man/MCMCirtKd.Rd
+3caeaa4d9321ec9631d97609f114465f *man/MCMCirtKdHet.Rd
+9483a7a5ca0d62bd6a8c48a24182715e *man/MCMCirtKdRob.Rd
+e857d4c2a3b7de7d70e1f26c9e8f6a38 *man/MCMClogit.Rd
+d2fc8449cd90706631ea2a763ef369d8 *man/MCMCmetrop1R.Rd
+076a43d4ed1e5c65ec06bfe27e8a3b32 *man/MCMCmixfactanal.Rd
+9ba1d602d6a9c90ddab157f9e3e49990 *man/MCMCmnl.Rd
+dea3128fd1118249d8c8525c67ad4889 *man/MCMCoprobit.Rd
+0ecd8d01cf09a72edd4e67a06b17027f *man/MCMCoprobitChange.Rd
+8014411538c970d9629a8488ce0b412d *man/MCMCordfactanal.Rd
+db8012e84acf5ae131e29042bf51856f *man/MCMCpoisson.Rd
+f7faabb370a5a432ce36871860c426bb *man/MCMCpoissonChange.Rd
+88e8572e82876f0c840818baf7f46a61 *man/MCMCprobit.Rd
+937ec3f54951e90191795373a490c221 *man/MCMCprobitChange.Rd
+1210ba8e3b786f3595be1d9b7b5201a4 *man/MCMCquantreg.Rd
+a9bb9b3e56b45f4bc8ae2ab8e675dfd5 *man/MCMCregress.Rd
+c6d6af6891968992b97bec70eaa3a31e *man/MCMCregressChange.Rd
02144cc6cca72312f4eafb013cfe23c3 *man/MCMCresidualBreakAnalysis.Rd
-c314e36afdd60449cf4f8ec08e00f80d *man/MCMCtobit.Rd
+e705eb63d2894d739bd565c5fbdde62c *man/MCMCtobit.Rd
a96eacae6b64827e40505661fb636cd4 *man/MCbinomialbeta.Rd
193401fb83df78594bb965c7bbac968f *man/MCmultinomdirichlet.Rd
4f2164ebeaaf83bde6a7f6c78b718d32 *man/MCnormalnormal.Rd
ecc40a241e849aacfa5ee3f17b8dd97d *man/MCpoissongamma.Rd
c04a5372b5a75cf67ef78d7e94d121e9 *man/Nethvote.Rd
-1aae2198de60fb0e4d568b0f2f4d994c *man/PErisk.Rd
+f5d931ef328d5fd614ea8547a3f0f62a *man/PErisk.Rd
034a4e9a54b8eff7dddcee5d0a923a02 *man/PostProbMod.Rd
aa99f8aa2b507bc8e1465a2b6e3993bd *man/QRSSVSplot.Rd
594c242058a0e3f9f2b8b5d4067c5b72 *man/QRSSVSsummary.Rd
-3ad4d7cb7b7b9fea628a16b8a9ab4053 *man/Rehnquist.Rd
-43eae86dc52b891ed0a9fc9d1c6887bb *man/SSVSquantreg.Rd
-e6849fc07eb2dd45dcc07300e7cd149d *man/Senate.Rd
-03bbb0adb292bf12b04fe466816d20a6 *man/SupremeCourt.Rd
+4ec36dfc16b0a856aa96dab7fd4e9075 *man/Rehnquist.Rd
+f76137bc031458ca844bb958d4c396bc *man/SSVSquantreg.Rd
+d1227e10d63f88c27fce979a25fc6299 *man/Senate.Rd
+121aaffb6bef26d405a56017759a8fca *man/SupremeCourt.Rd
1d8a9684290a7afeb7202f3797f5da71 *man/choicevar.Rd
50b0fec7baf798a609432b1578f65b41 *man/dirichlet.Rd
5203f2019ccb8b2dbf032c1689335ec3 *man/dtomog.Rd
@@ -118,13 +115,12 @@ b6f23a0e12fa9bb16f571261b265921f *man/invgamma.Rd
a5378f49cfb5f7ae58439b828e10fb84 *man/make.breaklist.Rd
629023d17bc9299449eec0b74a73340d *man/mptable.Rd
96ab1db8e6a5fb7c6640f06dd95a162c *man/noncenhypergeom.Rd
-170bc2168de62e08f7157a8cbc635e06 *man/plotChangepoint.Rd
-34999350c3ac48a2b5dbbb7766e5d463 *man/plotIntervention.Rd
-c64abc923f2c6f56a167a815c881ffb1 *man/plotState.Rd
+f01aec7464613309a58bcac6734d9321 *man/plotChangepoint.Rd
+ce50584a41db18441158e7e81c98a7b5 *man/plotState.Rd
6a88d877e88b5abfd7a66d659c9b4e6a *man/procrust.Rd
41feb17dda3b00cffa2f6176ba584b74 *man/readscythe.Rd
-f151473ede68e2cff21c06afd3060373 *man/testpanelGroupBreak.Rd
-24a4ef81a3573532c757c0c1dad393eb *man/testpanelSubjectBreak.Rd
+ef4daaf98db7ca5d55247575f47b67e2 *man/testpanelGroupBreak.Rd
+61444c0f42c5c15f7bf6cd72c4d3157c *man/testpanelSubjectBreak.Rd
254168bb5258fd00a7b766fb7d3bbfd9 *man/tomog.Rd
f5ba8eb700c67f188505a4585dd39f6b *man/topmodels.Rd
a4c560360ba912bb342a429e434bcc96 *man/vech.Rd
@@ -146,7 +142,6 @@ d4f36358e358d81d2d432c8cf23cc33d *src/MCMChierEI.cc
c938af28648eb358458f111d52119052 *src/MCMChlogit.cc
ab7fc8d513b7c66d9b777af3fda19133 *src/MCMChpoisson.cc
ea2d9f60541c419765dde67824ca9139 *src/MCMChregress.cc
-d0c3136bef3a5243dd984dab7f2df5ab *src/MCMCintervention.cc
fa7404ccc4ce053d13b2e0b9a86819b8 *src/MCMCirt1d.cc
abdf98660af4bed11d89e04074d9f1d1 *src/MCMCirtHier1d.cc
ac9b6d3f588b46cf009816bf028ee5ed *src/MCMCirtKdHet.cc
@@ -168,12 +163,12 @@ e2239f7d7503dc4eca812034f4bc1da8 *src/MCMCprobitChange.cc
a6a01c3906a8b92b9f1a12cfa92a1a63 *src/MCMCprobitres.cc
806af936660056856d4b908f8fb0cfa4 *src/MCMCquantreg.cc
c89523a7da37782992fe91d8364b1cd3 *src/MCMCregress.cc
-3fa981b1ba947770ff834c45237be2b1 *src/MCMCregressChange.cc
+67f37ba3bbf30e3f000e9ab885b7f2b5 *src/MCMCregressChange.cc
cfbd489b6fd8886e262a57ca9629c6da *src/MCMCresidualBreakAnalysis.cc
ef5566f824887d9bbef5b76f1b64d839 *src/MCMCrng.h
f28dc81e54a6ca54a2202a0d6f52ed40 *src/MCMCtobit.cc
-b3e224bb1a8452c3da347348b8e510cc *src/Makevars
43f4124319734194009441b15aff3695 *src/Makevars.in
+b3e224bb1a8452c3da347348b8e510cc *src/Makevars.win
6ba1ab7daac8f7149663686fd7abfd53 *src/SSVSquantreg.cc
5362ed41a42ce34cf7817f6bfd988935 *src/algorithm.h
80eab6b8072ae34de6c192af41324093 *src/datablock.h
diff --git a/NAMESPACE b/NAMESPACE
index 0eb1506..17b65ab 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,10 +1,18 @@
useDynLib(MCMCpack)
-import(coda)
-import(MASS)
-import(stats)
+import("coda")
+import("MASS")
+import("stats")
+import("utils")
+import("grDevices")
+import("lattice")
+import("methods")
+import("mcmc")
+import("quantreg")
+importFrom("graphics", "plot", "lines", "abline", "axTicks", "axis", "box", "layout",
+ "lcm", "legend", "par", "plot.new", "plot.window", "points",
+ "rect")
-export(
- BayesFactor,
+export(BayesFactor,
choicevar,
ddirichlet,
dinvgamma,
@@ -26,11 +34,9 @@ export(
MCMChlogit,
MCMChpoisson,
MCMChregress,
- MCMCintervention,
MCMCirt1d,
MCMCirtHier1d,
MCMCirtKd,
- MCMCirtKdHet,
MCMCirtKdRob,
MCMClogit,
MCMCmetrop1R,
@@ -52,7 +58,6 @@ export(
mptable,
plotState,
plotChangepoint,
- plotIntervention,
PostProbMod,
procrustes,
rdirichlet,
diff --git a/R/BayesFactors.R b/R/BayesFactors.R
index 3720fca..d24fbf2 100644
--- a/R/BayesFactors.R
+++ b/R/BayesFactors.R
@@ -149,8 +149,12 @@
BF.mat <- matrix(NA, M, M)
BF.log.mat <- matrix(NA, M, M)
rownames(BF.mat) <- colnames(BF.mat) <-
- rownames(BF.log.mat) <- colnames(BF.log.mat) <- model.names
- BF.logmarglike <- array(NA, M, dimnames=model.names)
+ rownames(BF.log.mat) <- colnames(BF.log.mat) <- model.names
+
+ BF.logmarglike <- array(NA, c(1, M), dimnames=list("logmarglike", model.names))
+ ## Bill Dunlap found that R did not warn about illegal dimnames in array() but rather would just disregard them.
+ ## So based on the patch by Martin Maechler, JHP changed the code. "Thu Feb 4 10:35:15 2016"
+ ## BF.logmarglike <- array(NA, M, dimnames=model.names)
BF.call <- NULL
for (i in 1:M){
@@ -163,7 +167,8 @@
BF.mat[i,j] <- exp(BF.log.mat[i,j])
}
else{
- warning(paste(model.names[i], " and ", model.names[j], " do not have exactly identical y data.\nBayes factors are not defined.\n", sep=""))
+ warning(paste(model.names[i], " and ", model.names[j],
+ " do not have exactly identical y data.\nBayes factors are not defined.\n", sep=""))
}
}
}
diff --git a/R/MCMCdynamicIRT1d-b.R b/R/MCMCdynamicIRT1d-b.R
index 78a25b1..703973f 100644
--- a/R/MCMCdynamicIRT1d-b.R
+++ b/R/MCMCdynamicIRT1d-b.R
@@ -120,7 +120,7 @@
}
else {
ab.starts[i,] <- coef(suppressWarnings(glm(local.y~theta.start,
- family=binomial(probit),
+ family=binomial(link="probit"),
control=glm.control(
maxit=8, epsilon=1e-3)
)))
diff --git a/R/MCMCdynamicIRT1d.R b/R/MCMCdynamicIRT1d.R
index 590aa48..8db5df5 100644
--- a/R/MCMCdynamicIRT1d.R
+++ b/R/MCMCdynamicIRT1d.R
@@ -118,7 +118,7 @@
}
else {
ab.starts[i,] <- coef(suppressWarnings(glm(local.y~theta.start,
- family=binomial(probit),
+ family=binomial(link="probit"),
control=glm.control(
maxit=8, epsilon=1e-3)
)))
diff --git a/R/MCMCintervention.R b/R/MCMCintervention.R
deleted file mode 100644
index 64e66cd..0000000
--- a/R/MCMCintervention.R
+++ /dev/null
@@ -1,231 +0,0 @@
-#########################################################
-## Internvetion Analysis using Changepoint Model
-#########################################################
-
-"MCMCintervention"<-
- function(y, data=parent.frame(), m = 1, intervention = 1,
- prediction.type=c("trend", "ar"),
- change.type = c("fixed", "random", "all"),
- b0 = 0, B0 = 0, c0 = 0.001, d0 = 0.001, sigma.mu = NA, sigma.var = NA,
- a = NULL, b = NULL,
- mcmc = 1000, burnin = 1000, thin = 1, verbose = 0,
- seed = NA, beta.start = NA, P.start = NA,
- marginal.likelihood = c("none", "Chib95"), ...){
-
- ## form response and model matrices
- y <- as.vector(y)
- n <- length(y)
- ns <- m + 1 # number of states
- if (prediction.type == "trend"){
- X <- matrix(cbind(1, c(1:n)), n, 2)
- xnames <- c("constant", "trend")
- }
- else if (prediction.type == "ar"){
- y1 <- y[1:(n/2)]
- ar1 <- arima(y1, c(1,0,0))
- muy0 <- as.numeric(ar1$coef[2])
- sigmay0 <- sqrt(as.numeric(as.numeric(ar1[2])/(1 - ar1$coef[1]^2)))
- y0 <- rnorm(1, muy0, sigmay0)
- X <- matrix(cbind(1, c(y0, y[-n])), n, 2)
- xnames <- c("constant", "lag.y")
- }
- else {
- X <- matrix(1, n, 1)
- xnames <- c("constant")
- }
- k <- ncol(X) # number of covariates
-
- ## check iteration parameters
- check.mcmc.parameters(burnin, mcmc, thin)
- totiter <- mcmc + burnin
- nstore <- mcmc/thin
- 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)
-
- ## prior
- mvn.prior <- form.mvn.prior(b0, B0, k)
- b0 <- mvn.prior[[1]]
- B0 <- mvn.prior[[2]]
- if (prediction.type == "ar"){
- if (b0[2] > 1|b0[2] < -1){
- stop("The prior of AR coefficient ",b0[2], " is outside the stationary region! \n")
- }
- }
- if (is.na(sigma.mu)|is.na(sigma.var)) {
- check.ig.prior(c0, d0)
- }
- else {
- d0 <- 2*(sigma.mu + sigma.mu^3/sigma.var)
- c0 <- 2*(1 + (d0/2)/sigma.mu)
- }
-
- ## get marginal likelihood argument
- marginal.likelihood <- match.arg(marginal.likelihood)
-
- ## following MCMCregress, set chib as binary
- logmarglike <- loglik <- NULL
- chib <- 0
- if (marginal.likelihood == "Chib95"){
- chib <- 1
- }
-
- ## initial values
- Y <- matrix(y, n, 1)
- if (m == 0){
- output <- MCMCregress(formula = Y ~ X-1, mcmc=mcmc,
- burnin=burnin, verbose=verbose, thin=thin,
- b0 = b0, B0 = solve(B0), c0 = c0, d0 = d0,
- marginal.likelihood = marginal.likelihood)
- attr(output, "y") <- y
- attr(output, "intervention") <- intervention
- yhatout <- output[, 1:2]%*%t(X)
- attr(output, "yhat") <- matrix(yhatout, nstore, n) ## X_j*beta_j
- attr(output, "yforepred") <- matrix(yhatout, nstore, n) ## yt|Y_t-1, theta
- attr(output, "ybackpred") <- matrix(yhatout, nstore, n) ## yt|Y_t+1, theta
-
- }
- else{ ## if m > 0
- A0 <- trans.mat.prior(m=m, n=n, a=a, b=b)
- Pstart <- check.P(P.start, m, a=a, b=b)
- ols <- lm(y~X-1)
- bols <- coef(ols)
- betastart <- matrix(rep(bols, ns), ns, k, byrow = TRUE)
- Sigmastart <- rep(summary(ols)$sigma^2, ns)
- statestart <- sort(sample(1:ns, n, replace=T))
- AR <- 0
- if (prediction.type == "ar"){
- AR <- 1
- }
- change <- 0
- betaout.N <- nstore*ns*k
- Sigmaout.N <- nstore*ns
- if (change.type == "fixed"){
- change <- 1
- betaout.N <- nstore*ns*k
- Sigmaout.N <- nstore
- }
- if (change.type == "random"){
- change <- 2
- betaout.N <- nstore*k
- Sigmaout.N <- nstore*ns
- }
-
- ## call C++ code to draw sample
- posterior <- .C("MCMCintervention",
- accept = as.double(0.0),
- betaout = as.double(rep(0.0, betaout.N)),
- Sigmaout = as.double(rep(0.0, Sigmaout.N)),
- 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)),
-
- yhatout = as.double(rep(0.0, nstore*n)),
- yerrorout = as.double(rep(0.0, nstore*n)),
- yforepredout = as.double(rep(0.0, nstore*n)),
- ybackpredout = 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),
- intervention = as.integer(intervention),
- burnin = as.integer(burnin),
- mcmc = as.integer(mcmc),
- thin = as.integer(thin),
- verbose = as.integer(verbose),
-
- lecuyer=as.integer(lecuyer),
- seedarray=as.integer(seed.array),
- lecuyerstream=as.integer(lecuyer.stream),
-
- betastart = as.double(betastart),
- Sigmastart = as.double(Sigmastart),
- Pstart = as.double(Pstart),
- statestart = as.integer(statestart),
-
- a = as.double(a),
- b = as.double(b),
- b0data = as.double(b0),
- B0data = as.double(B0),
- c0 = as.double(c0),
- d0 = as.double(d0),
- A0data = as.double(A0),
- logmarglikeholder = as.double(0.0),
- loglikeholder = as.double(0.0),
- ar = as.integer(AR),
- change = as.integer(change),
- 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
- beta.holder <- matrix(posterior$betaout, nstore, )
- Sigma.holder <- matrix(posterior$Sigmaout, nstore, )
- P.holder <- matrix(posterior$Pout, nstore, )
- s.holder <- matrix(posterior$sout, nstore, )
- ps.holder <- matrix(posterior$psout, n, )
-
- output1 <- mcmc(data=beta.holder, start=burnin+1, end=burnin + mcmc, thin=thin)
- if (change == 0){
- varnames(output1) <- sapply(c(1:ns),
- function(i){
- paste(c(xnames), "_regime", i, sep = "")
- })
- output2 <- mcmc(data=Sigma.holder, start=burnin+1, end=burnin + mcmc, thin=thin)
- varnames(output2) <- sapply(c(1:ns),
- function(i){
- paste(c("sigma2"), "_regime", i, sep = "")
- })
- }
- else if (change == 1){
- varnames(output1) <- sapply(c(1:ns),
- function(i){
- paste(c(xnames), "_regime", i, sep = "")
- })
- output2 <- mcmc(data=Sigma.holder, start=burnin+1, end=burnin + mcmc, thin=thin)
- names(output2) <- c("sigma2")
- }
- else{
- output2 <- mcmc(data=Sigma.holder, start=burnin+1, end=burnin + mcmc, thin=thin)
- varnames(output2) <- sapply(c(1:ns),
- function(i){
- paste(c("sigma2"), "_regime", i, sep = "")
- })
- }
- ## To check the acceptance rate (if AR == 1)
- accept <- posterior$accept
-
- output <- as.mcmc(cbind(output1, output2))
- attr(output, "title") <- "MCMCintervention Posterior Sample"
- attr(output, "intervention") <- intervention
- attr(output, "accept") <- accept
- attr(output, "y") <- y
- attr(output, "X") <- X
- attr(output, "m") <- m
- attr(output, "call") <- cl
- attr(output, "logmarglike") <- logmarglike
- attr(output, "loglik") <- loglik
- attr(output, "prob.state") <- ps.holder/nstore
- attr(output, "s.store") <- s.holder
- attr(output, "yhat") <- matrix(posterior$yhatout, nstore, n)## X_j*beta_j
- attr(output, "yerror") <- matrix(posterior$yerrorout, nstore, n)## y_j - X_j*beta_j
- attr(output, "yforepred") <- matrix(posterior$yforepredout, nstore, n)## yt|Y_t-1, theta
- attr(output, "ybackpred") <- matrix(posterior$ybackpredout, nstore, n)## yt|Y_t+1, theta
- }
- return(output)
-
-}## end of MCMC function
-
diff --git a/R/MCMCirt1d.R b/R/MCMCirt1d.R
index 935aca7..934a6f0 100644
--- a/R/MCMCirt1d.R
+++ b/R/MCMCirt1d.R
@@ -88,7 +88,7 @@
}
else {
ab.starts[i,] <- coef(suppressWarnings(glm(local.y~theta.start,
- family=binomial(probit),
+ family=binomial(link="probit"),
control=glm.control(
maxit=8, epsilon=1e-3)
)))
diff --git a/R/MCMCirtHier1d.R b/R/MCMCirtHier1d.R
index 713afa2..687de4c 100644
--- a/R/MCMCirtHier1d.R
+++ b/R/MCMCirtHier1d.R
@@ -89,7 +89,7 @@
}
else {
ab.starts[i,] <- coef(suppressWarnings(glm(local.y~theta.start,
- family=binomial(probit),
+ family=binomial(link="probit"),
control=glm.control(
maxit=8, epsilon=1e-3)
)))
diff --git a/R/MCMCirtKdRob.R b/R/MCMCirtKdRob.R
index 903c7db..8a1dd1b 100644
--- a/R/MCMCirtKdRob.R
+++ b/R/MCMCirtKdRob.R
@@ -214,7 +214,7 @@
if(Lambda.ineq.constraints[i,j]==0){
if (j==1){
probit.out <- glm(as.factor(X[X[,i]!=-999,i])~1,
- family=binomial(link=logit))
+ family=binomial(link="logit"))
probit.beta <- coef(probit.out)
Lambda[i,j] <- -1 * probit.beta[1]
}
diff --git a/R/MCMCmetrop1R.R b/R/MCMCmetrop1R.R
index 2073536..c1755d6 100644
--- a/R/MCMCmetrop1R.R
+++ b/R/MCMCmetrop1R.R
@@ -32,12 +32,19 @@
REPORT=10, maxit=500),
...){
+ ## following block creates an explicit copy of theta.init so that theta.init
+ ## is not modified in place by the code below
+ theta.init.0 <- rep(NA, length(theta.init))
+ for (i in 1:length(theta.init)){
+ theta.init.0[i] <- theta.init[i]
+ }
+
## error checking here
check.offset(list(...))
check.mcmc.parameters(burnin, mcmc, thin)
## form the tuning vector
- tune <- vector.tune(tune, length(theta.init))
+ tune <- vector.tune(tune, length(theta.init.0))
## form seed
seeds <- form.seeds(seed)
@@ -66,7 +73,7 @@
if (is.null(V)){
## find approx mode and Hessian using optim()
- opt.out <- optim(theta.init, maxfun,
+ opt.out <- optim(theta.init.0, maxfun,
control=optim.control,
lower=optim.lower, upper=optim.upper,
method=optim.method, hessian=TRUE, ...)
@@ -111,7 +118,7 @@
V <- tune %*% solve(-1*hess.new) %*% tune
}
else{ ## V is non NULL
- if (nrow(V) != ncol(V) || nrow(V) != length(theta.init)){
+ if (nrow(V) != ncol(V) || nrow(V) != length(theta.init.0)){
cat("V not of appropriate dimension.\n")
stop("Check V and theta.init and call MCMCmetrop1R() again. \n",
call.=FALSE)
@@ -127,7 +134,7 @@
}
## Call the C++ function to do the MCMC sampling
- sample <- .Call("MCMCmetrop1R_cc", userfun, as.double(theta.init),
+ sample <- .Call("MCMCmetrop1R_cc", userfun, as.double(theta.init.0),
my.env, as.integer(burnin), as.integer(mcmc),
as.integer(thin),
as.integer(verbose),
diff --git a/R/MCMCmixfactanal.R b/R/MCMCmixfactanal.R
index f686c77..d6d5aac 100644
--- a/R/MCMCmixfactanal.R
+++ b/R/MCMCmixfactanal.R
@@ -165,7 +165,7 @@
}
if (ncat[i] == 2){
probit.out <- glm(as.factor(X[X[,i]!=-999,i])~1,
- family=binomial(link=probit))
+ family=binomial(link="probit"))
probit.beta <- coef(probit.out)
Lambda[i,j] <- probit.beta[1]
}
diff --git a/R/MCMCordfactanal.R b/R/MCMCordfactanal.R
index 721ed63..655c432 100644
--- a/R/MCMCordfactanal.R
+++ b/R/MCMCordfactanal.R
@@ -156,7 +156,7 @@
if (j==1){
if (ncat[i] == 2){
probit.out <- glm(as.factor(X[X[,i]!=-999,i])~1,
- family=binomial(link=probit))
+ family=binomial(link="probit"))
probit.beta <- coef(probit.out)
Lambda[i,j] <- probit.beta[1]
}
diff --git a/R/MCMCprobit.R b/R/MCMCprobit.R
index a43b685..203f091 100644
--- a/R/MCMCprobit.R
+++ b/R/MCMCprobit.R
@@ -45,7 +45,7 @@
## starting values and priors
beta.start <- coef.start(beta.start, K, formula,
- family=binomial(link=probit), data)
+ family=binomial(link="probit"), data)
mvn.prior <- form.mvn.prior(b0, B0, K)
b0 <- mvn.prior[[1]]
B0 <- mvn.prior[[2]]
diff --git a/R/testpanelSubjectBreak.R b/R/testpanelSubjectBreak.R
index 7e4766c..9965bf8 100644
--- a/R/testpanelSubjectBreak.R
+++ b/R/testpanelSubjectBreak.R
@@ -56,7 +56,7 @@
out <- as.list(rep(NA, max.break))
if(nk > minimum){
for (k in 0:max.break){
- out[[k+1]] <- MCMCpack:::MCMCresidualBreakAnalysis(residual, m=k,
+ out[[k+1]] <- MCMCresidualBreakAnalysis(residual, m=k,
b0=b0, B0=B0, c0=c0, d0=d0, a=a, b=b,
burnin=burnin, mcmc=mcmc, thin=thin, verbose=verbose,
marginal.likelihood="Chib95")
diff --git a/README b/README
index b083a3c..fd100f1 100644
--- a/README
+++ b/README
@@ -4,9 +4,9 @@
// Authors
-Andrew D. Martin <admartin at wustl.edu>
+Andrew D. Martin <admart at umich.edu>
Kevin M. Quinn <kquinn at law.berkeley.edu>
-Jong Hee Park <jhp at uchicago.edu>
+Jong Hee Park <jongheepark at snu.ac.kr>
// Compilation
@@ -39,18 +39,14 @@ Neither the National Science Foundation, Washington University, or
Harvard University bear any responsibility for the content of this
package.
-Please contact Andrew D. Martin <admartin at wustl.edu> if
+Please contact Jong Hee Park <jongheepark at snu.ac.kr> if
you have any problems or questions.
--
-Andrew D. Martin, Ph.D.
-Charles Nagel Chair of Constitutional Law and Political Science
-Vice Dean and CERL Director, School of Law
-Washington University in St. Louis
+Jong Hee Park, Ph.D.
+Associate Professor in Dept. Political Science and International Relations
+Seoul National University
-(314) 935-5863 (Office)
-(314) 935-3836 (Fax)
-
-Email: admartin at wustl.edu
-WWW: http://adm.wustl.edu
+Email: jongheepark at snu.ac.kr
+WWW: http://jhp.snu.ac.kr
diff --git a/cleanup b/cleanup
index bc672f4..9418c29 100755
--- a/cleanup
+++ b/cleanup
@@ -1,3 +1,3 @@
#! /bin/sh
-rm -f config.log
+rm -f config.log src/Makevars
diff --git a/man/HMMpanelFE.Rd b/man/HMMpanelFE.Rd
index 8c6913a..8bad4a1 100644
--- a/man/HMMpanelFE.Rd
+++ b/man/HMMpanelFE.Rd
@@ -49,7 +49,7 @@
\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
+ 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.}
@@ -85,10 +85,7 @@
\details{
- \code{HMMpanelFE} simulates from the fixed-effect hidden Markov panel model introduced by Park (2011).
- The model takes the following form:
- \deqn{y_{it} = x'_{it} \beta + \varepsilon_{it}\;\; m = 1, \ldots, M}{y_it = x'_it * beta + epsilon_it, m = 1,...,M.}
- Unlike conventional fixed-effects models, individual effects and variances are assumed to be time-varying at the subject level:
+ \code{HMMpanelFE} simulates from the fixed-effect hidden Markov pbject level:
\deqn{\varepsilon_{it} \sim \mathcal{N}(\alpha_{im}, \sigma^2_{im})}{epsilon_it ~ N(alpha_im, sigma^2_im)}
We assume standard, semi-conjugate priors:
@@ -109,8 +106,6 @@
}
-\author{Jong Hee Park, \email{jhp at uchicago.edu}, \url{http://home.uchicago.edu/~jhp/}.}
-
\value{
An mcmc object that contains the posterior sample. This
object can be summarized by functions provided by the coda package.
@@ -119,8 +114,8 @@
}
\references{
- Jong Hee Park, 2011. ``A Unified Method for Dynamic and Cross-Sectional Heterogeneity:
- Introducing Hidden Markov Panel Models." Working Paper.
+ Jong Hee Park, 2012. ``Unified Method for Dynamic and Cross-Sectional Heterogeneity:
+ Introducing Hidden Markov Panel Models.'' \emph{American Journal of Political Science}.56: 1040-1054.
Siddhartha Chib. 1998. ``Estimation and comparison of multiple change-point models.''
\emph{Journal of Econometrics}. 86: 221-241.
@@ -169,8 +164,9 @@
weight <- Weight[j:(j+T-1)]
true.alpha1[i] <- rnorm(1, 0, true.sigma.alpha)
true.alpha2[i] <- -1*true.alpha1[i]
- y[j:(j+T-1)] <- ((1-weight)*true.mean + (1-weight)*rnorm(T, 0, true.sigma) + (1-weight)*true.alpha1[i]) +
- (weight*true.mean + weight*rnorm(T, 0, true.sigma) + weight*true.alpha2[i])
+ y[j:(j+T-1)] <- ((1-weight)*true.mean + (1-weight)*rnorm(T, 0, true.sigma) +
+ (1-weight)*true.alpha1[i]) +
+ (weight*true.mean + weight*rnorm(T, 0, true.sigma) + weight*true.alpha2[i])
j <- j + T
}
diff --git a/man/HMMpanelRE.Rd b/man/HMMpanelRE.Rd
index 50f0a5c..82ecaa9 100644
--- a/man/HMMpanelRE.Rd
+++ b/man/HMMpanelRE.Rd
@@ -147,16 +147,14 @@ HMMpanelRE(subject.id, time.id, y, X, W, m=1,
}
\references{
- Jong Hee Park, 2011. ``A Unified Method for Dynamic and Cross-Sectional Heterogeneity:
- Introducing Hidden Markov Panel Models." Working Paper.
+ Jong Hee Park, 2012. ``Unified Method for Dynamic and Cross-Sectional Heterogeneity:
+ Introducing Hidden Markov Panel Models.'' \emph{American Journal of Political Science}.56: 1040-1054.
Siddhartha Chib. 1998. ``Estimation and comparison of multiple change-point models.''
\emph{Journal of Econometrics}. 86: 221-241.
}
-\author{Jong Hee Park, \email{jhp at uchicago.edu}, \url{http://home.uchicago.edu/~jhp/}.}
-
\examples{
\dontrun{
## data generating
diff --git a/man/MCMCSVDreg.Rd b/man/MCMCSVDreg.Rd
index 6e96604..8e78954 100644
--- a/man/MCMCSVDreg.Rd
+++ b/man/MCMCSVDreg.Rd
@@ -149,9 +149,11 @@ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011.
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/}.
+ Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2006.
+ ``Output Analysis and Diagnostics for MCMC (CODA)'',
+ \emph{R News}. 6(1): 7-11.
+ \url{http://CRAN.R-project.org/doc/Rnews/Rnews_2006-1.pdf}.
+
}
diff --git a/man/MCMCbinaryChange.Rd b/man/MCMCbinaryChange.Rd
index 3678c76..496a116 100755
--- a/man/MCMCbinaryChange.Rd
+++ b/man/MCMCbinaryChange.Rd
@@ -84,8 +84,6 @@
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{
Jong Hee Park. 2011. ``Changepoint Analysis of Binary and Ordinal Probit Models:
An Application to Bank Rate Policy Under the Interwar Gold Standard." \emph{Political Analysis}. 19: 188-204.
@@ -95,10 +93,6 @@
\emph{Journal of Statistical Software}. 42(9): 1-21.
\url{http://www.jstatsoft.org/v42/i09/}.
- Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2002.
- \emph{Output Analysis and Diagnostics for MCMC (CODA)}.
- \url{http://www-fis.iarc.fr/coda/}.
-
Siddhartha Chib. 1995. ``Marginal Likelihood from the Gibbs Output.''
\emph{Journal of the American Statistical Association}. 90:
1313-1321.
@@ -116,12 +110,18 @@
y3 <- rbinom(120, 1, true.phi[3])
y <- as.ts(c(y1, y2, y3))
- model0 <- MCMCbinaryChange(y, m=0, c0=2, d0=2, mcmc=1000, burnin=1000, verbose=500, marginal.likelihood = "Chib95")
- 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")
+ model0 <- MCMCbinaryChange(y, m=0, c0=2, d0=2, mcmc=100, burnin=100, verbose=50,
+ marginal.likelihood = "Chib95")
+ model1 <- MCMCbinaryChange(y, m=1, c0=2, d0=2, mcmc=100, burnin=100, verbose=50,
+ marginal.likelihood = "Chib95")
+ model2 <- MCMCbinaryChange(y, m=2, c0=2, d0=2, mcmc=100, burnin=100, verbose=50,
+ marginal.likelihood = "Chib95")
+ model3 <- MCMCbinaryChange(y, m=3, c0=2, d0=2, mcmc=100, burnin=100, verbose=50,
+ marginal.likelihood = "Chib95")
+ model4 <- MCMCbinaryChange(y, m=4, c0=2, d0=2, mcmc=100, burnin=100, verbose=50,
+ marginal.likelihood = "Chib95")
+ model5 <- MCMCbinaryChange(y, m=5, c0=2, d0=2, mcmc=100, burnin=100, verbose=50,
+ marginal.likelihood = "Chib95")
print(BayesFactor(model0, model1, model2, model3, model4, model5))
diff --git a/man/MCMCdynamicEI.Rd b/man/MCMCdynamicEI.Rd
index b9f6e5e..8f9ab82 100644
--- a/man/MCMCdynamicEI.Rd
+++ b/man/MCMCdynamicEI.Rd
@@ -140,7 +140,7 @@ Strategies}. Gary King, Ori Rosen, and Martin A. Tanner (eds.). New
York: Cambridge University Press.
Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011.
- ``MCMCpack: Markov Chain Monte Carlo in R.'',
+ ``MCMCpack: Markov Chain Monte Carlo in R'',
\emph{Journal of Statistical Software}. 42(9): 1-21.
\url{http://www.jstatsoft.org/v42/i09/}.
@@ -150,12 +150,13 @@ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011.
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/}.
+ Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2006.
+ ``Output Analysis and Diagnostics for MCMC (CODA)'',
+ \emph{R News}. 6(1): 7-11.
+ \url{http://CRAN.R-project.org/doc/Rnews/Rnews_2006-1.pdf}.
-Jonathan C. Wakefield. 2004. ``Ecological Inference for 2 x 2 Tables.'' \emph{Journal
- of the Royal Statistical Society, Series A}. 167(3): 385445.
+ Jonathan C. Wakefield. 2004. ``Ecological Inference for 2 x 2 Tables.''
+ \emph{Journal of the Royal Statistical Society, Series A}. 167(3): 385445.
}
diff --git a/man/MCMCfactanal.Rd b/man/MCMCfactanal.Rd
index 2b215ba..83412bf 100644
--- a/man/MCMCfactanal.Rd
+++ b/man/MCMCfactanal.Rd
@@ -165,9 +165,10 @@ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011.
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/}.
+ Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2006.
+ ``Output Analysis and Diagnostics for MCMC (CODA)'',
+ \emph{R News}. 6(1): 7-11.
+ \url{http://CRAN.R-project.org/doc/Rnews/Rnews_2006-1.pdf}.
}
diff --git a/man/MCMChierEI.Rd b/man/MCMChierEI.Rd
index e8ca4cc..6846b91 100644
--- a/man/MCMChierEI.Rd
+++ b/man/MCMChierEI.Rd
@@ -153,9 +153,10 @@ MCMChierEI(r0, r1, c0, c1, burnin=5000, mcmc=50000, thin=1,
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/}.
+ Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2006.
+ ``Output Analysis and Diagnostics for MCMC (CODA)'',
+ \emph{R News}. 6(1): 7-11.
+ \url{http://CRAN.R-project.org/doc/Rnews/Rnews_2006-1.pdf}.
}
\examples{
diff --git a/man/MCMChlogit.Rd b/man/MCMChlogit.Rd
index 5265153..3d0cc33 100644
--- a/man/MCMChlogit.Rd
+++ b/man/MCMChlogit.Rd
@@ -182,9 +182,10 @@ nu=0.001, delta=0.001, FixOD=0, ...)}
Political Science Panel Data.'' Paper presented at the 2002 Annual Meeting
of the American Political Science Association.
- 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/}.
+ Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2006.
+ ``Output Analysis and Diagnostics for MCMC (CODA)'',
+ \emph{R News}. 6(1): 7-11.
+ \url{http://CRAN.R-project.org/doc/Rnews/Rnews_2006-1.pdf}.
Diggle P., Heagerty P., Liang K., and Zeger S. 2004. ``Analysis of Longitudinal Data.''
\emph{Oxford University Press}, 2sd Edition.
diff --git a/man/MCMChpoisson.Rd b/man/MCMChpoisson.Rd
index cc63275..ccb7bfb 100644
--- a/man/MCMChpoisson.Rd
+++ b/man/MCMChpoisson.Rd
@@ -181,9 +181,10 @@ nu=0.001, delta=0.001, FixOD=0, ...)}
Political Science Panel Data.'' Paper presented at the 2002 Annual Meeting
of the American Political Science Association.
- 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/}.
+ Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2006.
+ ``Output Analysis and Diagnostics for MCMC (CODA)'',
+ \emph{R News}. 6(1): 7-11.
+ \url{http://CRAN.R-project.org/doc/Rnews/Rnews_2006-1.pdf}.
}
\author{
diff --git a/man/MCMChregress.Rd b/man/MCMChregress.Rd
index 84fe37e..54361d0 100644
--- a/man/MCMChregress.Rd
+++ b/man/MCMChregress.Rd
@@ -164,9 +164,10 @@ nu=0.001, delta=0.001, ...)}
Political Science Panel Data.'' Paper presented at the 2002 Annual Meeting
of the American Political Science Association.
- 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/}.
+ Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2006.
+ ``Output Analysis and Diagnostics for MCMC (CODA)'',
+ \emph{R News}. 6(1): 7-11.
+ \url{http://CRAN.R-project.org/doc/Rnews/Rnews_2006-1.pdf}.
}
\author{
diff --git a/man/MCMCintervention.Rd b/man/MCMCintervention.Rd
deleted file mode 100644
index 0cc4d58..0000000
--- a/man/MCMCintervention.Rd
+++ /dev/null
@@ -1,209 +0,0 @@
-\name{MCMCintervention}
-\alias{MCMCintervention}
-
-\title{Markov Chain Monte Carlo for a linear Gaussian Multiple Changepoint Model}
-\description{
- This function generates a sample from the posterior distribution
- of a linear Gaussian model with multiple changepoints. The function uses
- the Markov chain Monte Carlo method of Chib (1998).
- The user supplies data and priors, and
- a sample from the posterior distribution is returned as an mcmc
- object, which can be subsequently analyzed with functions
- provided in the coda package.
-}
-
-\usage{MCMCintervention(y, data=parent.frame(), m = 1, intervention,
- prediction.type=c("trend","ar"), change.type=c("fixed", "random", "all"),
- b0 = 0, B0 = 0, c0 = 0.001, d0 = 0.001, sigma.mu = NA, sigma.var = NA,
- a = NULL, b = NULL, mcmc = 1000, burnin = 1000, thin = 1, verbose = 0,
- seed = NA, beta.start = NA, P.start = NA,
- marginal.likelihood = c("none", "Chib95"), ...)}
-
-\arguments{
- \item{y}{data.}
-
- \item{data}{Data frame.}
-
- \item{m}{The number of changepoints.}
-
- \item{intervention}{The timing of intervention measured by its place in the response vector.
- It should be larger than 0 and smaller than the length of the response vector. No default value is provided.}
-
- \item{prediction.type}{The type of local process. "trend" denotes the linear trend model and "ar"
- denotes AR(1) process. By default, MCMCintervention uses the linear trend model.}
-
- \item{change.type}{The tyep of parameteric breaks. "all" denotes that all parameters have breaks,
- "fixed" denotes that only the intercept and the slope have breaks, and "random" denotes that
- only the variance has breaks. By default, MCMCintervetnion assumes that all parameters have breaks.}
-
- \item{b0}{The prior mean of \eqn{\beta}{beta}. This can either be a
- scalar or a
- column vector with dimension equal to the number of betas. If this
- takes a scalar value, then that value will serve as the prior
- mean for all of the betas.}
-
- \item{B0}{The prior precision of \eqn{\beta}{beta}. This can either be a
- scalar or a square matrix with dimensions equal to the number of betas.
- If this
- takes a scalar value, then that value times an identity matrix serves
- as the prior precision of beta. Default value of 0 is equivalent to
- an improper uniform prior for beta.}
-
- \item{c0}{\eqn{c_0/2}{c0/2} is the shape parameter for the inverse
- Gamma prior on \eqn{\sigma^2}{sigma^2} (the variance of the
- disturbances). The amount of information in the inverse Gamma prior
- is something like that from \eqn{c_0}{c0} pseudo-observations.}
-
- \item{d0}{\eqn{d_0/2}{d0/2} is the scale parameter for the
- inverse Gamma prior on \eqn{\sigma^2}{sigma^2} (the variance of the
- disturbances). In constructing the inverse Gamma prior,
- \eqn{d_0}{d0} acts like the sum of squared errors from the
- \eqn{c_0}{c0} pseudo-observations.}
-
- \item{sigma.mu}{The mean of the
- inverse Gamma prior on \eqn{\sigma^2}{sigma^2}.
- \eqn{sigma.mu}{sigma.mu} and \eqn{sigma.var}{sigma.var} allow users to choose the inverse Gamma prior
- by choosing its mean and variance. }
-
- \item{sigma.var}{The variacne of the
- inverse Gamma prior on \eqn{\sigma^2}{sigma^2}.
- \eqn{sigma.mu}{sigma.mu} and \eqn{sigma.var}{sigma.var} allow users to choose the inverse Gamma prior
- by choosing its mean and variance. }
-
- \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 burnin.}
-
- \item{thin}{The thinning interval used in the simulation. The number of
- MCMC iterations must be divisible by this value.}
-
- \item{verbose}{A switch which determines whether or not the progress of
- the sampler is printed to the screen. If \code{verbose} is greater
- than 0 the iteration number, the
- \eqn{\beta}{beta} vector, and the error variance are printed to
- the screen every \code{verbose}th iteration.}
-
- \item{seed}{The seed for the random number generator. If NA, the Mersenne
- Twister generator is used with default seed 12345; if an integer is
- passed it is used to seed the Mersenne twister. The user can also
- pass a list of length two to use the L'Ecuyer random number generator,
- which is suitable for parallel computation. The first element of the
- list is the L'Ecuyer seed, which is a vector of length six or NA (if NA
- a default seed of \code{rep(12345,6)} is used). The second element of
- list is a positive substream number. See the MCMCpack
- specification for more details.}
-
- \item{beta.start}{The starting values for the \eqn{\beta}{beta} vector.
- This can either be a scalar or a
- column vector with dimension equal to the number of betas.
- The default value of of NA will use the MLE
- estimate of \eqn{\beta}{beta} as the starting value. If this is a
- scalar, that value will serve as the starting value
- mean for all of the betas.}
-
- \item{P.start}{The starting values for the transition matrix.
- A user should provide a square matrix with dimension equal to the number of states.
- By default, draws from the \code{Beta(0.9, 0.1)} are used to construct
- a proper transition matrix for each raw except the last raw.}
-
- \item{marginal.likelihood}{How should the marginal likelihood be
- calculated? Options are: \code{none} in which case the marginal
- likelihood will not be calculated, and
- \code{Chib95} in which case the method of Chib (1995) is used.}
-
- \item{...}{further arguments to be passed}
-}
-
-\value{
- An mcmc object that contains the posterior sample. This
- object can be summarized by functions provided by the coda package.
- The object contains an attribute \code{prob.state} storage matrix that contains the probability of \eqn{state_i}{state_i} for each period, the
- log-likelihood of the model (\code{loglike}), and
- the log-marginal likelihood of the model (\code{logmarglike}).
-}
-
-\details{
- \code{MCMCintervention} simulates from the posterior distribution of
- a binary model with multiple changepoints.
-
- The model takes the following form:
- \deqn{y_t = x_t ' \beta_i + I(s_t = i)\varepsilon_{t},\;\; i = 1, \ldots, k}{y_t = x_t'beta_i + I(s_t = i)epsilon_t, i = 1,...,k.}
- Where \eqn{k}{k} is the number of states and \eqn{I(s_t = i)}{I(s_t = i)} is an indicator function
- that becomes 1 when a state at \eqn{t}{t} is \eqn{i}{i} and otherwise 0.
-
- The errors are assumed to be Gaussian in each regime:
- \deqn{I(s_t = i)\varepsilon_{t} \sim \mathcal{N}(0, \sigma^2_i)}{I(s_t = i)epsilon_t ~ N(0,
- sigma^2_i)}
-
- We assume standard, semi-conjugate priors:
- \deqn{\beta_i \sim \mathcal{N}(b_0,B_0^{-1}),\;\; i = 1, \ldots, k}{beta_i ~ N(b0,B0^(-1)), i = 1,...,k.}
- And:
- \deqn{\sigma^{-2}_i \sim \mathcal{G}amma(c_0/2, d_0/2),\;\; i = 1, \ldots, k}{sigma^(-2)_i ~
- Gamma(c0/2, d0/2), i = 1,...,k.}
- Where \eqn{\beta_i}{beta_i} and \eqn{\sigma^{-2}_i}{sigma^(-2)_i} are assumed
- \emph{a priori} independent.
-
- The simulation proper is done in compiled C++ code to maximize efficiency. }
-
-\references{
- Jong Hee Park. 2012. "A Change-point Approach to Intervention Analysis Using Bayesian Inference"
- Presented at the 2012 Annual Meeting of Korean Statistical Society.
-
- Siddhartha Chib. 1998. "Estimation and comparison of multiple change-point models."
- \emph{Journal of Econometrics}. 86: 221-241.
-}
-
-\examples{
- \dontrun{
- Nile.std <- (Nile - mean(Nile))/sd(Nile)
- set.seed(1973)
- b0 <- matrix(c(0, 0), 2, 1);
- B0 <- diag(rep(0.25, 2))
- c0 = 2; d0 = 1
-
- ## Model Comparison
- ar0 <- MCMCintervention(Nile.std, m=0, prediction.type="ar", change.type = "all",
- b0=b0, B0=B0, c0=c0, d0=d0, mcmc = 1000, burnin = 1000, verbose = 500,
- intervention = 29, marginal.likelihood = "Chib95")
- ar1all <- MCMCintervention(Nile.std, m=1, prediction.type="ar", change.type = "all",
- b0=b0, B0=B0, c0=c0, d0=d0, mcmc = 1000, burnin = 1000, verbose = 500,
- intervention = 29, marginal.likelihood = "Chib95")
- ar1fixed <- MCMCintervention(Nile.std, m=1, prediction.type="ar", change.type = "fixed",
- b0=b0, B0=B0, c0=c0, d0=d0, mcmc = 1000, burnin = 1000, verbose = 500,
- intervention = 29, marginal.likelihood = "Chib95")
- ar1random <- MCMCintervention(Nile.std, m=1, prediction.type="ar", change.type = "random",
- b0=b0, B0=B0, c0=c0, d0=d0, mcmc = 1000, burnin = 1000, verbose = 500,
- intervention = 29, marginal.likelihood = "Chib95")
- tr0 <- MCMCintervention(Nile.std, m=0, prediction.type="trend", change.type = "all",
- b0=b0, B0=B0, c0=c0, d0=d0, mcmc = 1000, burnin = 1000, verbose = 500,
- intervention = 29, marginal.likelihood = "Chib95")
- tr1all <- MCMCintervention(Nile.std, m=1, prediction.type="trend", change.type = "all",
- b0=b0, B0=B0, c0=c0, d0=d0, mcmc = 1000, burnin = 1000, verbose = 500,
- intervention = 29, marginal.likelihood = "Chib95")
- tr1fixed <- MCMCintervention(Nile.std, m=1, prediction.type="trend", change.type = "fixed",
- b0=b0, B0=B0, c0=c0, d0=d0, mcmc = 1000, burnin = 1000, verbose = 500,
- intervention = 29, marginal.likelihood = "Chib95")
- tr1random <- MCMCintervention(Nile.std, m=1, prediction.type="trend", change.type = "random",
- b0=b0, B0=B0, c0=c0, d0=d0, mcmc = 1000, burnin = 1000, verbose = 500,
- intervention = 29, marginal.likelihood = "Chib95")
-
- BayesFactor(ar0, ar1all, ar1fixed, ar1random, tr0, tr1all, tr1fixed, tr1random)
-
- par(mfrow=c(1,3))
- plotState(ar1fixed, start=1871, main="Hidden Regime Change")
- plotIntervention(ar1fixed, start=1871, main="Forward Analysis", alpha= 0.5, ylab="Nile River flow", xlab="Year")
- plotIntervention(ar1fixed, forward=FALSE, start=1871, main="Backward Analysis", alpha= 0.5, ylab="Nile River flow", xlab="Year")
- }
-}
-
-\keyword{models}
-
-\seealso{\code{\link{plotIntervention}}}
diff --git a/man/MCMCirt1d.Rd b/man/MCMCirt1d.Rd
index 10f2ed7..e33150d 100644
--- a/man/MCMCirt1d.Rd
+++ b/man/MCMCirt1d.Rd
@@ -192,9 +192,10 @@ MCMCirt1d(datamatrix, theta.constraints=list(), burnin = 1000,
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/}.
+ Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2006.
+ ``Output Analysis and Diagnostics for MCMC (CODA)'',
+ \emph{R News}. 6(1): 7-11.
+ \url{http://CRAN.R-project.org/doc/Rnews/Rnews_2006-1.pdf}.
Douglas Rivers. 2004. ``Identification of Multidimensional Item-Response
Models." Stanford University, typescript.
diff --git a/man/MCMCirtHier1d.Rd b/man/MCMCirtHier1d.Rd
index 5d66e13..fb8e1d1 100644
--- a/man/MCMCirtHier1d.Rd
+++ b/man/MCMCirtHier1d.Rd
@@ -224,7 +224,7 @@ The parameters of interest are
As is the case with all measurement models, make sure that you have plenty
of free memory, especially when storing the item parameters.
}
- \author{Michael Malecki, \email{malecki at wustl.edu}, \url{http://malecki.wustl.edu}.}
+ \author{Michael Malecki, \email{mike at crunch.io}, \url{https://github.com/malecki}.}
\references{
James H. Albert. 1992. ``Bayesian Estimation of Normal Ogive Item Response
@@ -248,10 +248,11 @@ The parameters of interest are
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/}.
-
+ Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2006.
+ ``Output Analysis and Diagnostics for MCMC (CODA)'',
+ \emph{R News}. 6(1): 7-11.
+ \url{http://CRAN.R-project.org/doc/Rnews/Rnews_2006-1.pdf}.
+
Douglas Rivers. 2004. ``Identification of Multidimensional Item-Response
Models." Stanford University, typescript.
}
@@ -290,8 +291,9 @@ Xjdata <- data.frame(presparty= c(1,1,0,1,1,1,1,0,0),
## Randomly 10% Missing -- this affects the expansion parameter, increasing
## the variance of the (unidentified) latent parameter alpha.
-scMiss <- SupremeCourt
-scMiss[matrix(as.logical(rbinom(nrow(SupremeCourt)*ncol(SupremeCourt), 1, .1)), dim(SupremeCourt))] <- NA
+ scMiss <- SupremeCourt
+ scMiss[matrix(as.logical(rbinom(nrow(SupremeCourt)*ncol(SupremeCourt), 1, .1)),
+ dim(SupremeCourt))] <- NA
posterior1.miss <- MCMCirtHier1d(t(scMiss),
burnin=80000, mcmc=10000, thin=20,
diff --git a/man/MCMCirtKd.Rd b/man/MCMCirtKd.Rd
index ca72491..a62e585 100644
--- a/man/MCMCirtKd.Rd
+++ b/man/MCMCirtKd.Rd
@@ -203,9 +203,10 @@ MCMCirtKd(datamatrix, dimensions, item.constraints=list(),
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/}.
+ Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2006.
+ ``Output Analysis and Diagnostics for MCMC (CODA)'',
+ \emph{R News}. 6(1): 7-11.
+ \url{http://CRAN.R-project.org/doc/Rnews/Rnews_2006-1.pdf}.
Douglas Rivers. 2004. ``Identification of Multidimensional Item-Response
Models." Stanford University, typescript.
diff --git a/man/MCMCirtKdHet.Rd b/man/MCMCirtKdHet.Rd
index 5cdedcf..19b35e9 100644
--- a/man/MCMCirtKdHet.Rd
+++ b/man/MCMCirtKdHet.Rd
@@ -145,8 +145,8 @@ store.sigma = TRUE, drop.constant.items = TRUE)
\emph{Political Analysis.} 18: 151-171.
}
\author{
-Benjamin E. Lauderdale, \email{blauderd at princeton.edu},
-\url{http://www.princeton.edu/~blauderd/}.
+Benjamin E. Lauderdale, \email{b.e.lauderdale at lse.ac.uk},
+\url{http://www.benjaminlauderdale.net}.
Modified from \code{\link[MCMCpack]{MCMCirtKd}} and \code{\link[MCMCpack]{MCMCordfactanal}}. Suggestions for additional options are welcome.
}
@@ -183,7 +183,8 @@ for (i in 1:102){
legend(x="topleft",legend=c("Point sizes proportional to estimated legislator",
"variance under heteroskedastic model.","Some legislators with large variance have",
"more extreme estimated ideal points under the","heteroskedastic model because their",
-"deviations from the party line are attributable","to idiosyncrasy rather than moderation."),cex=0.5)
+"deviations from the party line are attributable","to idiosyncrasy rather than moderation."),
+cex=0.5)
}
}
diff --git a/man/MCMCirtKdRob.Rd b/man/MCMCirtKdRob.Rd
index 9ee7406..fe68dc8 100644
--- a/man/MCMCirtKdRob.Rd
+++ b/man/MCMCirtKdRob.Rd
@@ -307,11 +307,12 @@ MCMCirtKdRob(datamatrix, dimensions, item.constraints=list(),
Radford Neal. 2003. ``Slice Sampling'' (with discussion). \emph{Annals of
Statistics}, 31: 705-767.
-
- 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/}.
-
+
+ Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2006.
+ ``Output Analysis and Diagnostics for MCMC (CODA)'',
+ \emph{R News}. 6(1): 7-11.
+ \url{http://CRAN.R-project.org/doc/Rnews/Rnews_2006-1.pdf}.
+
Douglas Rivers. 2004. ``Identification of Multidimensional Item-Response
Models." Stanford University, typescript.
}
diff --git a/man/MCMClogit.Rd b/man/MCMClogit.Rd
index d916c82..514ff79 100644
--- a/man/MCMClogit.Rd
+++ b/man/MCMClogit.Rd
@@ -143,9 +143,10 @@ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011.
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/}.
+ Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2006.
+ ``Output Analysis and Diagnostics for MCMC (CODA)'',
+ \emph{R News}. 6(1): 7-11.
+ \url{http://CRAN.R-project.org/doc/Rnews/Rnews_2006-1.pdf}.
}
diff --git a/man/MCMCmetrop1R.Rd b/man/MCMCmetrop1R.Rd
index edcd04e..9b95b54 100644
--- a/man/MCMCmetrop1R.Rd
+++ b/man/MCMCmetrop1R.Rd
@@ -217,9 +217,10 @@ mode. This last calculation is done via an initial call to
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/}.
+ Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2006.
+ ``Output Analysis and Diagnostics for MCMC (CODA)'',
+ \emph{R News}. 6(1): 7-11.
+ \url{http://CRAN.R-project.org/doc/Rnews/Rnews_2006-1.pdf}.
Christian P. Robert and George Casella. 2004. \emph{Monte Carlo
Statistical Methods}. 2nd Edition. New York: Springer.
diff --git a/man/MCMCmixfactanal.Rd b/man/MCMCmixfactanal.Rd
index 8ae66cb..92a2788 100644
--- a/man/MCMCmixfactanal.Rd
+++ b/man/MCMCmixfactanal.Rd
@@ -234,9 +234,10 @@ MCMCmixfactanal(x, factors, lambda.constraints=list(),
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/}.
+ Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2006.
+ ``Output Analysis and Diagnostics for MCMC (CODA)'',
+ \emph{R News}. 6(1): 7-11.
+ \url{http://CRAN.R-project.org/doc/Rnews/Rnews_2006-1.pdf}.
}
diff --git a/man/MCMCmnl.Rd b/man/MCMCmnl.Rd
index 5c3f670..5ff37db 100644
--- a/man/MCMCmnl.Rd
+++ b/man/MCMCmnl.Rd
@@ -163,9 +163,10 @@ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011.
Radford Neal. 2003. ``Slice Sampling'' (with discussion). \emph{Annals of
Statistics}, 31: 705-767.
- 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/}.
+ Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2006.
+ ``Output Analysis and Diagnostics for MCMC (CODA)'',
+ \emph{R News}. 6(1): 7-11.
+ \url{http://CRAN.R-project.org/doc/Rnews/Rnews_2006-1.pdf}.
Siddhartha Chib, Edward Greenberg, and Yuxin Chen. 1998.
``MCMC Methods for Fitting and Comparing Multinomial Response Models."
diff --git a/man/MCMCoprobit.Rd b/man/MCMCoprobit.Rd
index 5859650..9ed6844 100644
--- a/man/MCMCoprobit.Rd
+++ b/man/MCMCoprobit.Rd
@@ -143,9 +143,10 @@ that can be used to analyze the posterior sample.
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/}
+ Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2006.
+ ``Output Analysis and Diagnostics for MCMC (CODA)'',
+ \emph{R News}. 6(1): 7-11.
+ \url{http://CRAN.R-project.org/doc/Rnews/Rnews_2006-1.pdf}.
}
diff --git a/man/MCMCoprobitChange.Rd b/man/MCMCoprobitChange.Rd
index 55bc1ae..cc2e394 100644
--- a/man/MCMCoprobitChange.Rd
+++ b/man/MCMCoprobitChange.Rd
@@ -140,7 +140,8 @@
\references{
Jong Hee Park. 2011. ``Changepoint Analysis of Binary and Ordinal Probit Models:
- An Application to Bank Rate Policy Under the Interwar Gold Standard." \emph{Political Analysis}. 19: 188-204.
+ An Application to Bank Rate Policy Under the Interwar Gold Standard."
+ \emph{Political Analysis}. 19: 188-204.
Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011.
``MCMCpack: Markov Chain Monte Carlo in R.'',
@@ -173,17 +174,15 @@ formula <- y ~ x1
## fit multiple models with a varying number of breaks
out1 <- MCMCoprobitChange(formula, m=1,
- mcmc=1000, burnin=1000, thin=1, tune=c(.5, .5), verbose=1000,
+ mcmc=100, burnin=100, thin=1, tune=c(.5, .5), verbose=100,
b0=0, B0=10, marginal.likelihood = "Chib95")
out2 <- MCMCoprobitChange(formula, m=2,
- mcmc=1000, burnin=1000, thin=1, tune=c(.5, .5, .5), verbose=1000,
- b0=0, B0=10, marginal.likelihood = "Chib95")
-out3 <- MCMCoprobitChange(formula, m=3,
- mcmc=1000, burnin=1000, thin=1, tune=c(.5, .5, .5, .5), verbose=1000,
+ mcmc=100, burnin=100, thin=1, tune=c(.5, .5, .5), verbose=100,
b0=0, B0=10, marginal.likelihood = "Chib95")
-## find the most reasonable one
-BayesFactor(out1, out2, out3)
+## Do model comparison
+## NOTE: the chain should be run longer than this example!
+BayesFactor(out1, out2)
## draw plots using the "right" model
plotState(out1)
diff --git a/man/MCMCordfactanal.Rd b/man/MCMCordfactanal.Rd
index 8c8d662..7c8ec32 100644
--- a/man/MCMCordfactanal.Rd
+++ b/man/MCMCordfactanal.Rd
@@ -200,9 +200,10 @@ MCMCordfactanal(x, factors, lambda.constraints=list(),
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/}.
+ Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2006.
+ ``Output Analysis and Diagnostics for MCMC (CODA)'',
+ \emph{R News}. 6(1): 7-11.
+ \url{http://CRAN.R-project.org/doc/Rnews/Rnews_2006-1.pdf}.
}
diff --git a/man/MCMCpoisson.Rd b/man/MCMCpoisson.Rd
index 055763c..fd23da9 100644
--- a/man/MCMCpoisson.Rd
+++ b/man/MCMCpoisson.Rd
@@ -117,9 +117,11 @@ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011.
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/}.
+ Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2006.
+ ``Output Analysis and Diagnostics for MCMC (CODA)'',
+ \emph{R News}. 6(1): 7-11.
+ \url{http://CRAN.R-project.org/doc/Rnews/Rnews_2006-1.pdf}.
+
}
diff --git a/man/MCMCpoissonChange.Rd b/man/MCMCpoissonChange.Rd
index bf4bfd9..46c1fbb 100644
--- a/man/MCMCpoissonChange.Rd
+++ b/man/MCMCpoissonChange.Rd
@@ -114,8 +114,6 @@
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{
Jong Hee Park. 2010. ``Structural Change in the U.S. Presidents' Use of Force Abroad.''
\emph{American Journal of Political Science} 54: 766-782.
@@ -126,7 +124,7 @@
Siddhartha Chib. 1998. ``Estimation and comparison of multiple change-point models.''
\emph{Journal of Econometrics}. 86: 221-241.
-Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011.
+ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011.
``MCMCpack: Markov Chain Monte Carlo in R.'',
\emph{Journal of Statistical Software}. 42(9): 1-21.
\url{http://www.jstatsoft.org/v42/i09/}.
diff --git a/man/MCMCprobit.Rd b/man/MCMCprobit.Rd
index 22b5808..40d6037 100644
--- a/man/MCMCprobit.Rd
+++ b/man/MCMCprobit.Rd
@@ -123,9 +123,11 @@ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011.
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/}.
+ Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2006.
+ ``Output Analysis and Diagnostics for MCMC (CODA)'',
+ \emph{R News}. 6(1): 7-11.
+ \url{http://CRAN.R-project.org/doc/Rnews/Rnews_2006-1.pdf}.
+
}
diff --git a/man/MCMCprobitChange.Rd b/man/MCMCprobitChange.Rd
index 090b937..1c60e6d 100644
--- a/man/MCMCprobitChange.Rd
+++ b/man/MCMCprobitChange.Rd
@@ -118,8 +118,6 @@
}
-\author{Jong Hee Park, \email{jhp at uchicago.edu}, \url{http://home.uchicago.edu/~jhp/}.}
-
\references{
Jong Hee Park. 2011. ``Changepoint Analysis of Binary and Ordinal Probit Models:
An Application to Bank Rate Policy Under the Interwar Gold Standard." \emph{Political Analysis}. 19: 188-204.
diff --git a/man/MCMCquantreg.Rd b/man/MCMCquantreg.Rd
index 80e387a..e5ee711 100644
--- a/man/MCMCquantreg.Rd
+++ b/man/MCMCquantreg.Rd
@@ -111,9 +111,11 @@ The likelihood is formed based on assuming independent Asymmetric Laplace distri
Keming Yu and Jin Zhang. 2005. ``A Three Parameter Asymmetric Laplace Distribution and it's extensions.''
\emph{Communications in Statistics - Theory and Methods}, 34, 1867-1879.
- 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/}.}
+ Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2006.
+ ``Output Analysis and Diagnostics for MCMC (CODA)'',
+ \emph{R News}. 6(1): 7-11.
+ \url{http://CRAN.R-project.org/doc/Rnews/Rnews_2006-1.pdf}.
+}
\examples{
\dontrun{
diff --git a/man/MCMCregress.Rd b/man/MCMCregress.Rd
index 42b9f36..5ebd289 100644
--- a/man/MCMCregress.Rd
+++ b/man/MCMCregress.Rd
@@ -143,16 +143,18 @@ MCMCregress(formula, data = NULL, burnin = 1000, mcmc = 10000,
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/}.
+ Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2006.
+ ``Output Analysis and Diagnostics for MCMC (CODA)'',
+ \emph{R News}. 6(1): 7-11.
+ \url{http://CRAN.R-project.org/doc/Rnews/Rnews_2006-1.pdf}.
}
\examples{
\dontrun{
line <- list(X = c(-2,-1,0,1,2), Y = c(1,3,3,3,5))
-posterior <- MCMCregress(Y~X, b0=0, B0 = 0.1, sigma.mu = 5, sigma.var = 25, data=line, verbose=1000)
+posterior <- MCMCregress(Y~X, b0=0, B0 = 0.1,
+ sigma.mu = 5, sigma.var = 25, data=line, verbose=1000)
plot(posterior)
raftery.diag(posterior)
summary(posterior)
diff --git a/man/MCMCregressChange.Rd b/man/MCMCregressChange.Rd
index b919fb2..d7da0b2 100644
--- a/man/MCMCregressChange.Rd
+++ b/man/MCMCregressChange.Rd
@@ -156,9 +156,11 @@
\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/}.
+ Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2006.
+ ``Output Analysis and Diagnostics for MCMC (CODA)'',
+ \emph{R News}. 6(1): 7-11.
+ \url{http://CRAN.R-project.org/doc/Rnews/Rnews_2006-1.pdf}.
+
}
\examples{
@@ -187,11 +189,16 @@ sigma.mu=sd(y)
sigma.var=var(y)
## models
-model1 <- MCMCregressChange(formula, m=1, b0=b0, B0=B0, sigma.mu=sigma.mu, sigma.var=sigma.var, marginal.likelihood="Chib95")
-model2 <- MCMCregressChange(formula, m=2, b0=b0, B0=B0, sigma.mu=sigma.mu, sigma.var=sigma.var, marginal.likelihood="Chib95")
-model3 <- MCMCregressChange(formula, m=3, b0=b0, B0=B0, sigma.mu=sigma.mu, sigma.var=sigma.var, marginal.likelihood="Chib95")
-model4 <- MCMCregressChange(formula, m=4, b0=b0, B0=B0, sigma.mu=sigma.mu, sigma.var=sigma.var, marginal.likelihood="Chib95")
-model5 <- MCMCregressChange(formula, m=5, b0=b0, B0=B0, sigma.mu=sigma.mu, sigma.var=sigma.var, marginal.likelihood="Chib95")
+model1 <- MCMCregressChange(formula, m=1, b0=b0, B0=B0, mcmc=100, burnin=100,
+ sigma.mu=sigma.mu, sigma.var=sigma.var, marginal.likelihood="Chib95")
+model2 <- MCMCregressChange(formula, m=2, b0=b0, B0=B0, mcmc=100, burnin=100,
+ sigma.mu=sigma.mu, sigma.var=sigma.var, marginal.likelihood="Chib95")
+model3 <- MCMCregressChange(formula, m=3, b0=b0, B0=B0, mcmc=100, burnin=100,
+ sigma.mu=sigma.mu, sigma.var=sigma.var, marginal.likelihood="Chib95")
+model4 <- MCMCregressChange(formula, m=4, b0=b0, B0=B0, mcmc=100, burnin=100,
+ sigma.mu=sigma.mu, sigma.var=sigma.var, marginal.likelihood="Chib95")
+model5 <- MCMCregressChange(formula, m=5, b0=b0, B0=B0, mcmc=100, burnin=100,
+ sigma.mu=sigma.mu, sigma.var=sigma.var, marginal.likelihood="Chib95")
print(BayesFactor(model1, model2, model3, model4, model5))
plotState(model1)
diff --git a/man/MCMCtobit.Rd b/man/MCMCtobit.Rd
index bcd32d7..8e4dfa8 100644
--- a/man/MCMCtobit.Rd
+++ b/man/MCMCtobit.Rd
@@ -132,9 +132,10 @@ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011.
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/}.
+ Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2006.
+ ``Output Analysis and Diagnostics for MCMC (CODA)'',
+ \emph{R News}. 6(1): 7-11.
+ \url{http://CRAN.R-project.org/doc/Rnews/Rnews_2006-1.pdf}.
Siddhartha Chib. 1992. ``Bayes inference in the Tobit censored regression model."
\emph{Journal of Econometrics}. 51:79-99.
@@ -144,7 +145,7 @@ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011.
}
\author{Ben Goodrich, \email{goodrich.ben at gmail.com},
- \url{http://www.people.fas.harvard.edu/~goodrich/}}
+ \url{http://www.columbia.edu/~bg2382}}
\seealso{
\code{\link[coda]{plot.mcmc}},
diff --git a/man/PErisk.Rd b/man/PErisk.Rd
index d032452..433b6f1 100644
--- a/man/PErisk.Rd
+++ b/man/PErisk.Rd
@@ -33,7 +33,7 @@ Political Economic Risk Data from 62 Countries in 1987.
\source{
Mike Alvarez, Jose Antonio Cheibub, Fernando Limongi, and Adam
Przeworski. 1999. ``ACLP Political and Economic Database.''
- \url{http://www.ssc.upenn.edu/~cheibub/data/}.
+ \url{https://sites.google.com/site/joseantoniocheibub/datasets/democracy-and-development-aclp}.
Witold J. Henisz. 2002. ``The Political Constraint Index (POLCON)
Dataset.'' \\
@@ -41,7 +41,7 @@ Political Economic Risk Data from 62 Countries in 1987.
Monty G. Marshall, Ted Robert Gurr, and Barbara Harff. 2002. ``State
Failure Task Force Problem Set.''
- \url{http://www.cidcm.umd.edu/inscr/stfail/index.htm}.
+
}
\references{
Kevin M. Quinn. 2004. ``Bayesian Factor Analysis for Mixed Ordinal
diff --git a/man/Rehnquist.Rd b/man/Rehnquist.Rd
index d147dd8..fb9d9da 100644
--- a/man/Rehnquist.Rd
+++ b/man/Rehnquist.Rd
@@ -26,7 +26,7 @@ data(SupremeCourt)
\source{
Harold J. Spaeth. 2005. \emph{Original United States Supreme Court Database:
- 1953-2004 Terms.} \url{http://www.as.uky.edu/polisci/ulmerproject/sctdata.htm}.
+ 1953-2004 Terms.} \url{http://facweb.knowlton.ohio-state.edu/pviton/support/codebook-c.html}.
}
\keyword{datasets}
diff --git a/man/SSVSquantreg.Rd b/man/SSVSquantreg.Rd
index a24ef01..4844438 100644
--- a/man/SSVSquantreg.Rd
+++ b/man/SSVSquantreg.Rd
@@ -102,10 +102,12 @@ If it is certain that a predictor should be included, all predictors specified a
Keming Yu and Jin Zhang. 2005. "A Three Parameter Asymmetric Laplace Distribution and it's extensions."
\emph{Communications in Statistics - Theory and Methods}, 34, 1867-1879.
- 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/}.}
-
+ Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2006.
+ ``Output Analysis and Diagnostics for MCMC (CODA)'',
+ \emph{R News}. 6(1): 7-11.
+ \url{http://CRAN.R-project.org/doc/Rnews/Rnews_2006-1.pdf}.
+}
+
\examples{
\dontrun{
diff --git a/man/Senate.Rd b/man/Senate.Rd
index 4c17dee..396b012 100644
--- a/man/Senate.Rd
+++ b/man/Senate.Rd
@@ -23,7 +23,7 @@ data(Senate)
\source{
Keith Poole. 2005. \emph{106th Roll Call Vote Data}.
- \url{http://voteview.uh.edu/}.
+ \url{http://www.voteview.com/}.
}
\keyword{datasets}
diff --git a/man/SupremeCourt.Rd b/man/SupremeCourt.Rd
index 5520340..d8b4ad4 100644
--- a/man/SupremeCourt.Rd
+++ b/man/SupremeCourt.Rd
@@ -25,7 +25,7 @@ data(SupremeCourt)
\source{
Harold J. Spaeth. 2005. \emph{Original United States Supreme Court Database:
- 1953-2004 Terms.} \url{http://www.as.uky.edu/polisci/ulmerproject/sctdata.htm}.
+ 1953-2004 Terms.} \url{http://supremecourtdatabase.org}.
}
\keyword{datasets}
diff --git a/man/plotChangepoint.Rd b/man/plotChangepoint.Rd
index 1f58d92..8b84439 100755
--- a/man/plotChangepoint.Rd
+++ b/man/plotChangepoint.Rd
@@ -25,8 +25,6 @@
\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{MCMCpoissonChange}}, \code{\link{MCMCbinaryChange}}}
diff --git a/man/plotIntervention.Rd b/man/plotIntervention.Rd
deleted file mode 100644
index 18c1582..0000000
--- a/man/plotIntervention.Rd
+++ /dev/null
@@ -1,37 +0,0 @@
-\name{plotIntervention}
-\alias{plotIntervention}
-\title{Plot of intervention analysis}
-\description{Plot the results of changepoint internvetion analysis}
-
-\usage{
- plotIntervention(mcmcout, forward = TRUE, start = 1,
- alpha = 0.05, col = "red", main="", ylab="", xlab="")
-}
-
-\arguments{
-
-\item{mcmcout}{The \code{mcmcout} object containing the posterior density sample from changepoint intervention analysis.}
-
-\item{forward}{If TRUE, draw the result of forward prediction. If FALSE,
-draw the result of backward prediction.}
-
-\item{start}{The time of the first observation to be shown in the time series plot.}
-
-\item{alpha}{The size of the prediction error band in the plot. By default, alpha = 0.05, which means the 95 percent Bayesian credible interval.}
-
-\item{col}{The color of predicted data.}
-
-\item{main}{Title of the plot}
-
-\item{ylab}{Label for the y-axis.}
-
-\item{xlab}{Label for the x-axis.}
-
-}
-
-\references{
- Jong Hee Park. 2012. "A Change-point Approach to Intervention Analysis Using Bayesian Inference"
- Presented at the 2012 Annual Meeting of Korean Statistical Society.
-}
-
-\seealso{\code{\link{MCMCintervention}}}
diff --git a/man/plotState.Rd b/man/plotState.Rd
index 2c55e49..975b3b0 100755
--- a/man/plotState.Rd
+++ b/man/plotState.Rd
@@ -4,7 +4,8 @@
\description{Plot the posterior probability that each time point is in each state.}
\usage{
- plotState(mcmcout, main="Posterior Regime Probability", ylab=expression(paste("Pr(", S[t], "= k |", Y[t], ")")),
+ 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)
}
@@ -26,8 +27,6 @@
}
-\author{Jong Hee Park, \email{jhp at uchicago.edu}, \url{http://home.uchicago.edu/~jhp/}.}
-
\keyword{hplot}
\seealso{\code{\link{MCMCpoissonChange}}, \code{\link{MCMCbinaryChange}}}
diff --git a/man/testpanelGroupBreak.Rd b/man/testpanelGroupBreak.Rd
index 02650dd..7ff505d 100644
--- a/man/testpanelGroupBreak.Rd
+++ b/man/testpanelGroupBreak.Rd
@@ -94,16 +94,14 @@ testpanelGroupBreak(subject.id, time.id, resid, m=1,
}
\references{
- Jong Hee Park, 2011. ``A Unified Method for Dynamic and Cross-Sectional Heterogeneity:
- Introducing Hidden Markov Panel Models." Working Paper.
+ Jong Hee Park, 2012. ``Unified Method for Dynamic and Cross-Sectional Heterogeneity:
+ Introducing Hidden Markov Panel Models.'' \emph{American Journal of Political Science}.56: 1040-1054.
Siddhartha Chib. 1998. ``Estimation and comparison of multiple change-point models.''
\emph{Journal of Econometrics}. 86: 221-241.
}
-\author{Jong Hee Park, \email{jhp at uchicago.edu}, \url{http://home.uchicago.edu/~jhp/}.}
-
\examples{
\dontrun{
## data generating
diff --git a/man/testpanelSubjectBreak.Rd b/man/testpanelSubjectBreak.Rd
index 98cb96e..955c99b 100644
--- a/man/testpanelSubjectBreak.Rd
+++ b/man/testpanelSubjectBreak.Rd
@@ -89,8 +89,6 @@
}
-\author{Jong Hee Park, \email{jhp at uchicago.edu}, \url{http://home.uchicago.edu/~jhp/}.}
-
\value{
The returned object is a matrix containing log marginal likelihoods for all HMMs.
The dimension of the returned object is the number of panel subjects by max.break + 1.
@@ -98,8 +96,8 @@
}
\references{
- Jong Hee Park, 2011. ``A Unified Method for Dynamic and Cross-Sectional Heterogeneity:
- Introducing Hidden Markov Panel Models." Working Paper.
+ Jong Hee Park, 2012. ``Unified Method for Dynamic and Cross-Sectional Heterogeneity:
+ Introducing Hidden Markov Panel Models.'' \emph{American Journal of Political Science}.56: 1040-1054.
Siddhartha Chib. 1998. ``Estimation and comparison of multiple change-point models.''
\emph{Journal of Econometrics}. 86: 221-241.
@@ -147,8 +145,9 @@
weight <- Weight[j:(j+T-1)]
true.alpha1[i] <- rnorm(1, 0, true.sigma.alpha)
true.alpha2[i] <- -1*true.alpha1[i]
- y[j:(j+T-1)] <- ((1-weight)*true.mean + (1-weight)*rnorm(T, 0, true.sigma) + (1-weight)*true.alpha1[i]) +
- (weight*true.mean + weight*rnorm(T, 0, true.sigma) + weight*true.alpha2[i])
+ y[j:(j+T-1)] <- ((1-weight)*true.mean + (1-weight)*rnorm(T, 0, true.sigma) +
+ (1-weight)*true.alpha1[i]) +
+ (weight*true.mean + weight*rnorm(T, 0, true.sigma) + weight*true.alpha2[i])
j <- j + T
}
diff --git a/src/MCMCintervention.cc b/src/MCMCintervention.cc
deleted file mode 100644
index 4d316b4..0000000
--- a/src/MCMCintervention.cc
+++ /dev/null
@@ -1,1664 +0,0 @@
-////////////////////////////////////////////////////////////////////
-// MCMCintervention.cc is a C++ code to estimate
-// linear regression changepoint model
-//
-// Jong Hee Park
-// Department of Political Science and International Relations
-// Seoul National University
-// jongheepark at snu.ac.kr
-//
-// Written 03/03/2009
-////////////////////////////////////////////////////////////////////
-
-#ifndef MCMCREGRESSCHANGE_CC
-#define MCMCREGRESSCHANGE_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;
-
-
-double lndinvgamma_pjh (const double x, const double shape, const double scale){
- double log_density = shape *::log(scale) - lngammafn(shape) - (shape + 1) * ::log(x) - (scale/x);
- return (log_density);
-}
-
-template <typename RNGTYPE>
-Matrix<double> gaussian_state_fixedBeta_sampler(rng<RNGTYPE>& stream,
- const int m,
- const Matrix<double>& Y,
- const Matrix<double>& X,
- const Matrix<double>& beta,
- const Matrix<double>& Sigma,
- const Matrix<double>& P){
-
- const int ns = m + 1;
- const int n = Y.rows();
- Matrix<double> F(n, ns);
- Matrix<double> pr1(ns, 1);
- pr1[0] = 1;
- Matrix<double> py(ns, 1);
- Matrix<double> pstyt1(ns, 1);
- Matrix<int> s(n, 1); // holder for state variables
- Matrix<double> ps = Matrix<double>(n, ns); // holder for state probabilities
- for (int tt=0; tt<n ; ++tt){
- Matrix<double> mu = X(tt,_)*beta; //k by 1 vector
- for (int j = 0; j< ns; ++j){
- py[j] = dnorm(Y[tt], mu[0], sqrt(Sigma[j]));
- }
- if (tt==0) pstyt1 = pr1;
- else {
- pstyt1 = ::t(F(tt-1,_)*P); // make it an ns by 1 matrix
- }
- Matrix<double> unnorm_pstyt = pstyt1%py;
- const Matrix<double> pstyt = unnorm_pstyt/sum(unnorm_pstyt); // pstyt = Pr(st|Yt)
- for (int j=0; j<ns ; ++j) F(tt,j) = pstyt(j);
-
- }// end of F matrix filtering
- ps(n-1,_) = F(n-1,_);
- s(n-1) = ns;
- Matrix<double> pstyn = Matrix<double>(ns, 1);
- double pone = 0.0;
- int tt = n-2;
- while (tt >= 0){
- int st = s(tt+1);
- Matrix<double> Pst_1 = ::t(P(_,st-1));
- Matrix<double> unnorm_pstyn = F(tt,_)%Pst_1;
- pstyn = unnorm_pstyn/sum(unnorm_pstyn);
- if (st==1) s(tt) = 1;
- else{
- pone = pstyn(st-2);
- if(stream.runif() < pone) s(tt) = st-1;
- else s(tt) = st;
- }
- ps(tt,_) = pstyn;
- --tt;
- }// end of while loop
- Matrix<double> 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
-
-
-template <typename RNGTYPE>
-Matrix<double> gaussian_state_sampler(rng<RNGTYPE>& stream,
- const int m,
- const Matrix<double>& Y,
- const Matrix<double>& X,
- const Matrix<double>& beta,
- const Matrix<double>& Sigma,
- const Matrix<double>& P){
-
- const int ns = m + 1;
- const int n = Y.rows();
-
- // P is a (m+1 by m+1) transition matrix
- // F matrix contains all the information of Pr(st|Yt)
- Matrix<double> F(n, ns);
- Matrix<double> pr1(ns, 1);
- pr1[0] = 1;
- Matrix<double> py(ns, 1);
- Matrix<double> pstyt1(ns, 1);
- Matrix<int> s(n, 1); // holder for state variables
- Matrix<double> ps = Matrix<double>(n, ns); // holder for state probabilities
-
- //
- // Forward sampling: update F matrix
- //
- for (int tt=0; tt<n ; ++tt){
- Matrix<double> mu = X(tt,_)*::t(beta); //k by 1 vector
- for (int j = 0; j< ns; ++j){
- py[j] = dnorm(Y[tt], mu[j], sqrt(Sigma[j]));
- }
- if (tt==0) pstyt1 = pr1;
- else {
- pstyt1 = ::t(F(tt-1,_)*P); // make it an ns by 1 matrix
- }
- Matrix<double> unnorm_pstyt = pstyt1%py;
- const Matrix<double> pstyt = unnorm_pstyt/sum(unnorm_pstyt); // pstyt = Pr(st|Yt)
- for (int j=0; j<ns ; ++j) F(tt,j) = pstyt(j);
-
- }// end of F matrix filtering
-
- //
- // Backward recursions: sample s using F and a transition matrix P
- //
- ps(n-1,_) = F(n-1,_); // we know last elements of ps and s
- s(n-1) = ns; // the last period is s(*n-1) but one down in C location
-
- Matrix<double> pstyn = Matrix<double>(ns, 1);
- double pone = 0.0;
- int tt = n-2;
- while (tt >= 0){
- int st = s(tt+1);
- Matrix<double> Pst_1 = ::t(P(_,st-1)); // prob of being at a previous state
- Matrix<double> unnorm_pstyn = F(tt,_)%Pst_1;
- pstyn = unnorm_pstyn/sum(unnorm_pstyn); // normalize into a prob. density
- if (st==1) s(tt) = 1; // If this is the first period, state should be 1.
- // Otherwise, draw a state from a discrete prob. distribution("pstyn")
- // using the inverse CDF method.
- else{
- pone = pstyn(st-2);
- if(stream.runif() < pone) s(tt) = st-1;
- else s(tt) = st;
- }
- ps(tt,_) = pstyn;
- --tt;
- }// end of while loop
-
- // name and report outputs
- Matrix<double> 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
-
-////////////////////////////////////////////
-// MCMCintervention random effect changes only
-////////////////////////////////////////////
-template <typename RNGTYPE>
-void MCMCintervention_random_impl(rng<RNGTYPE>& stream,
- const double m, const int intervention,
- const Matrix<>& Y, const Matrix<>& X,
- Matrix<>& beta, Matrix<>& Sigma, Matrix<>& P, Matrix<int>& s,
- Matrix<>& b0, Matrix<>& B0,
- const double c0, const double d0,
- const Matrix<>& A0,
- unsigned int burnin, unsigned int mcmc, unsigned int thin,
- unsigned int verbose, bool chib, bool ar,
- Matrix<>& beta_store, Matrix<>& Sigma_store,
- Matrix<>& P_store, Matrix<>& ps_store, Matrix<int>& s_store,
- double& logmarglike,
- Matrix<>& yhat_mat,
- Matrix<>& yerror_mat,
- Matrix<>& yfore_pred_mat,
- Matrix<>& yback_pred_mat)
-{
- // define constants and form cross-product matrices
- const int tot_iter = burnin + mcmc; //total iterations
- const int nstore = mcmc / thin; // number of draws to store
- const int n = Y.rows();
- const int ns = m + 1; // number of states
- const int k = X.cols();
- const Matrix<> B0inv = invpd(B0);
- Matrix<> sigma(ns, 1);
-
- //MCMC loop
- unsigned int count = 0;
- unsigned int reject = 0;
- double rejectionrate = 0;
- Matrix <> sum_e(ns, 1);
-
- for (int iter = 0; iter < tot_iter; ++iter){
- //////////////////////
- // 1. Sample beta and Sigma
- //////////////////////
- int beta_count = 0;
- Matrix<int> nstate(ns, 1);
- Matrix<double> XtX(k, k);
- Matrix<double> XtY(k, 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
- nstate[j] = nstate[j] + 1;
- }// end of if
- }// end of int i<n
- beta_count = beta_count + nstate[j];
-
- Matrix<double> yj = Y((beta_count - nstate[j]), 0, (beta_count - 1), 0);
- Matrix<double> Xj = X((beta_count - nstate[j]), 0, (beta_count - 1), k-1);
-
- // SIGMA UPDATE
- double shape = (c0 + (double)nstate[j])/2;
- Matrix<> yhat_j = Xj*beta;
- Matrix<double> ej = yj - yhat_j;
- Matrix<double> sum_ej = t(ej)*ej;
- double scale =(d0 + sum_ej[0])/2;
-
- Sigma[j] = 1/stream.rgamma(shape, scale);
- sigma[j] = sqrt(Sigma[j]);
-
- if (iter >= burnin && ((iter % thin)==0)){
- yhat_mat(count, (beta_count - nstate[j]), count, (beta_count - 1)) = yhat_j(_,0);
- yerror_mat(count, (beta_count - nstate[j]), count, (beta_count - 1)) = ej(_,0);
- }
-
- // CONSTANT BETA UPDATE
- XtX = XtX + (crossprod(Xj))/Sigma[j];
- XtY = XtY + ::t(Xj)*yj/Sigma[j];
- }
- Matrix<double> Bn = invpd(B0 + XtX);
- Matrix<double> bn = Bn*(B0*b0 + XtY);
- if (ar == 1){
- Matrix<double> beta_can = stream.rmvnorm(bn, Bn);
- if (beta_can(1) > 1 | beta_can(1) < -1){
- // Rprintf("\n AR coefficient %10.5f is outside the stationary region! \n", beta_can(1));
- ++reject;
- }
- else{
- beta = beta_can;
- }
- }
- else{
- beta = stream.rmvnorm(bn, Bn);
- }
- //////////////////////
- // 2. 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) + (double)nstate[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);
- }
-
- //////////////////////
- // 3. Sample s
- //////////////////////
- Matrix<double> F(n, ns);
- Matrix<double> pr1(ns, 1);
- pr1[0] = 1;
- Matrix<double> py(ns, 1);
- Matrix<double> pstyt1(ns, 1);
- Matrix<double> ps = Matrix<double>(n, ns); // holder for state probabilities
-
- //
- // Forward sampling: update F matrix
- //
- for (int tt=0; tt<n ; ++tt){
- Matrix<double> mu = X(tt,_)*beta; //k by 1 vector
- for (int j = 0; j< ns; ++j){
- py[j] = dnorm(Y[tt], mu[0], sigma[j]);
- }
- if (tt==0) pstyt1 = pr1;
- else {
- pstyt1 = ::t(F(tt-1,_)*P); // make it an ns by 1 matrix
- }
- Matrix<double> unnorm_pstyt = pstyt1%py;
-
- /////////////////////////////////////////////////////////////////////
- // Prediction of future outcomes based on pre-intervention state
- /////////////////////////////////////////////////////////////////////
- if (tt==(intervention - 1)&&iter >= burnin && ((iter % thin)==0)){
- // Forward prediction
- Matrix <> yfore_pred(1, n);
- for (int ttt=tt; ttt<n ; ++ttt){
- int ss = s(tt-1); // 1 unit before the intervention
- mu = X(ttt,_)*beta; //k by 1 vector
- yfore_pred(ttt) = stream.rnorm(mu(0), sigma[ss-1]);
- }
- yfore_pred_mat(count, _) = yfore_pred(0, _);
- // backward prediction
- Matrix <> yback_pred(1, n);
- for (int ttt=tt; ttt>=0 ; --ttt){
- int ss = s(tt+1);
- mu = X(ttt,_)*beta; //k by 1 vector
- yback_pred(ttt) = stream.rnorm(mu(0), sigma[ss-1]);
- }
- yback_pred_mat(count, _) = yback_pred(0, _);
- }
-
- const Matrix<double> pstyt = unnorm_pstyt/sum(unnorm_pstyt); // pstyt = Pr(st|Yt)
- for (int j=0; j<ns ; ++j) F(tt,j) = pstyt(j);
-
- }// end of F matrix filtering
-
- //
- // Backward recursions: sample s using F and a transition matrix P
- //
- ps(n-1,_) = F(n-1,_); // we know last elements of ps and s
- s(n-1) = ns; // the last period is s(*n-1) but one down in C location
-
- Matrix<double> pstyn = Matrix<double>(ns, 1);
- double pone = 0.0;
- int tt = n-2;
- while (tt >= 0){
- int st = s(tt+1);
- Matrix<double> Pst_1 = ::t(P(_,st-1)); // prob of being at a previous state
- Matrix<double> unnorm_pstyn = F(tt,_)%Pst_1;
- pstyn = unnorm_pstyn/sum(unnorm_pstyn); // normalize into a prob. density
- if (st==1) s(tt) = 1; // If this is the first period, state should be 1.
- // Otherwise, draw a state from a discrete prob. distribution("pstyn")
- // using the inverse CDF method.
- else{
- pone = pstyn(st-2);
- if(stream.runif() < pone) s(tt) = st-1;// jump from tt-1 to tt
- else s(tt) = st;// stay
- }
- ps(tt,_) = pstyn;
- --tt;
- }// end of while loop
-
-
- // load draws into sample array
- if (iter >= burnin && ((iter % thin)==0)){
- for (int i=0; i<k; ++i)
- beta_store(count,i) = beta[i];// stored by the order of (11, 12, 13, 21, 22, 23)
- for (int i=0; i<ns; ++i)
- Sigma_store(count, i) = Sigma[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,_); // add two matrices
- s_store(count,_) = s;
-
- ++count; // count when (iter % *thin)==0
-
- } // end of if(iter...)
-
-
- // print output to stdout
- if(verbose > 0 && iter % verbose == 0){
- Rprintf("\nMCMCintervention iteration %i of %i \n", (iter+1), tot_iter);
- if (ar == 1 ){
- double rejectionrate = (double)reject/(double)(iter+1);
- Rprintf("\n\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n");
- Rprintf("The acceptance rate was %3.5f", 1 - rejectionrate);
- Rprintf("\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n");
- }
- for (int j = 0;j<ns; ++j){
- Rprintf("\n nstate[j] %10.5f", static_cast<double>(nstate[j]));
- }
- Rprintf("\n beta \n");
- for (int i = 0; i<k; ++i){
- Rprintf("%10.5f\n", beta(i));
- }
- Rprintf("\n sigma^2 \n");
- for (int i = 0; i<ns; ++i){
- Rprintf("%10.5f\n", Sigma(i));
- }
- }
-
- }// end MCMC loop
-
-
- if(chib ==1){\
- double rejectionrate = (double)reject/(double)nstore;
- if(rejectionrate > .05){
- Rprintf("The acceptance rate is too low.\n");
- }
- else {
- Matrix<double> beta_st(k, 1);
- Matrix<double> betast = meanc(beta_store); //meanc(beta_store)=(11, 12, 13, 21, 22, 23)
- for (int i = 0; i<k ; ++i){
- beta_st(i) = betast(i);
- //Rprintf("beta_st(i) %10.5f\n", beta_st(i));
- }
- Matrix<double> Sigma_st = meanc(Sigma_store);
- Matrix<double> sigma_st(ns, 1);
- for (int j = 0; j<ns ; ++j){
- sigma_st(j) = sqrt(Sigma_st(j));
- //Rprintf("sigma_st(j) %10.5f\n", sigma_st(j));
- }
-
- Matrix<double> P_vec_st = meanc(P_store);
- const Matrix<double> P_st(ns, ns);
- for (int j = 0; j< ns*ns; ++j){
- P_st[j] = P_vec_st[j];
- }
-
- //////////////////////
- // 1. pdf.beta
- //////////////////////
- Matrix<double> density_beta(nstore, 1);
- for (int iter = 0; iter<nstore; ++iter){
- Matrix<double> XtX(k, k);
- Matrix<double> XtY(k, 1);
- Matrix<int> nstate(ns, 1); // contains total numbers of each state
- int beta_count = 0;
-
- for (int j = 0; j<ns ; ++j){
- for (int i = 0; i<n; ++i){
- if (s_store(iter, i) == (j+1)) { // since j starts from 0
- nstate[j] = nstate[j] + 1;
- }// end of if
- }// end of int i<n
- beta_count = beta_count + nstate[j];
-
- const Matrix<double> yj = Y((beta_count - nstate[j]), 0, (beta_count - 1), 0);
- const Matrix<double> Xj = X((beta_count - nstate[j]), 0, (beta_count - 1), k-1);
- const double precision = 1.0/Sigma_store(iter, j);
- XtX = XtX + (crossprod(Xj))*precision;
- XtY = XtY + ::t(Xj)*yj*precision;
- }
- Matrix<double> Bn = invpd(B0 + XtX);
- Matrix<double> bn = Bn*(B0*b0 + XtY);
- if (k == 1){
- density_beta(iter) = dnorm(beta_st(0), bn(0), sqrt(Bn(0)));
- }
- else{
- density_beta(iter) = ::exp(lndmvn(beta_st, bn, Bn));
- }
- }
-
- double pdf_beta = log(prod(meanc(density_beta)));
-
- ////////////////////// ////////////////////// //////////////////////
- // 2. pdf.Sigma|beta_st, S, P
- ////////////////////// ////////////////////// //////////////////////
- Matrix<double> density_Sigma(nstore, ns);
- for (int iter = 0; iter<nstore; ++iter){
- // STEP 2.1 S|y, beta.st, Sigma, P
- Matrix <double> Sout = gaussian_state_fixedBeta_sampler(stream, m, Y, X, beta_st, Sigma, P);
- Matrix <double> s = Sout(_, 0);
-
- // STEP 2.2 Sigma|y, beta.st, S, P
- int beta_count = 0;
- Matrix<int> nstate(ns, 1); // contains total numbers of each state
-
- for (int j = 0; j <ns ; ++j){
- for (int i = 0; i<n; ++i){
- if (s[i] == (j+1)) { // since j starts from 0
- nstate[j] = nstate[j] + 1;
- }// end of if
- }// end of int i<n
- beta_count = beta_count + nstate[j];
-
- Matrix<double> yj = Y((beta_count - nstate[j]), 0, (beta_count - 1), 0);
- Matrix<double> Xj = X((beta_count - nstate[j]), 0, (beta_count - 1), k-1);
- Matrix<double> ej = yj - Xj*beta_st;
- Matrix<double> sum_ej = ::t(ej)*ej;
- double shape = (c0 + (double)nstate[j])/2;
- double scale =(d0 + sum_ej[0])/2;
-
- Sigma[j] = stream.rigamma(shape, scale);
- density_Sigma(iter, j) = ::exp(lndinvgamma_pjh(Sigma_st[j], shape, scale));
- }// end of sampling beta and Sigma
-
- // STEP 2.3 P|S
- 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) + (double)nstate[j] - 1;
- shape2 = A0(j,j+1) + 1; //
- P(j,j) = stream.rbeta(shape1, shape2);
- P(j,j+1) = 1 - P(j,j);
- }
- }// end of pdf.Sigma
-
- double pdf_Sigma = log(prod(meanc(density_Sigma)));
-
- //////////////////////
- // 3. pdf.P|beta_st, Sigma_st, S
- //////////////////////
- Matrix<double> density_P(nstore, ns);
- for (int iter = 0; iter < nstore; ++iter){
-
- // STEP 2.1 S|y, beta.st, Sigma, P
- Matrix <double> Sout = gaussian_state_fixedBeta_sampler(stream, m, Y, X, beta_st, Sigma_st, P);
- Matrix <double> s = Sout(_, 0);
-
- double shape1 = 0;
- double shape2 = 0;
- P(ns-1, ns-1) = 1; //no jump at the last state
-
- // compute addN
- Matrix <double> P_addN(ns, 1);
- for (int j = 0; j<ns; ++j){
- for (int i = 0; i<n; ++i){
- if (s(i) == (j+1)) { // since j starts from 0
- P_addN[j] = P_addN[j] + 1;
- }// end of if
- } // end of for i
- } // end of for j
-
- // evaluate P
- 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<double> F = Matrix<double>(n, ns);
- Matrix<double> like(n, 1);
- Matrix<double> pr1 = Matrix<double>(ns, 1);
- pr1[0] = 1;
- Matrix<double> py(ns, 1);
- Matrix<double> pstyt1(ns, 1);
-
- for (int t=0; t<n ; ++t){
- int yt = (int) Y[t];
- Matrix<double> mu = X(t,_)*beta_st;
- for (int j = 0; j< ns; ++j){
- py[j] = dnorm(Y[t], mu[0], sigma_st[j]);
- }
- if (t==0) pstyt1 = pr1;
- else {
- pstyt1 = ::t(F(t-1,_)*P_st);
- }
- Matrix<double> unnorm_pstyt = pstyt1%py;
- Matrix<double> pstyt = unnorm_pstyt/sum(unnorm_pstyt);
- for (int j=0; j<ns ; ++j){
- F(t,j) = pstyt(j);
- }
- like[t] = sum(unnorm_pstyt);
- }// end of likelihood computation
-
- const double loglike = sum(log(like));
-
- //////////////////////
- // log prior ordinate
- //////////////////////
- double density_beta_prior = 0.0;
- Matrix<double> density_Sigma_prior(ns, 1);
- Matrix<double> density_P_prior(ns, 1);
- density_P[ns-1] = 1; //
- if (k == 1){
- density_beta_prior = log(dnorm(beta_st(0), b0(0), sqrt(B0inv(0))));
- }
- else{
- density_beta_prior = lndmvn(beta_st, b0, B0inv);
- }
- for (int j=0; j<ns ; ++j){
- density_Sigma_prior[j] = lndinvgamma_pjh(Sigma_st[j], c0/2, d0/2);
- }
-
- 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 = density_beta_prior + sum(density_Sigma_prior) + sum(density_P_prior);
- logmarglike = (loglike + logprior) - (pdf_beta + pdf_Sigma + pdf_P);
- if(verbose > 0){
- Rprintf("logmarglike = %10.5f\n", logmarglike);
- Rprintf("loglike = %10.5f\n", loglike);
- Rprintf("logprior = %10.5f\n", logprior);
- Rprintf("pdf_beta = %10.5f\n", pdf_beta);
- Rprintf("pdf_Sigma = %10.5f\n", pdf_Sigma);
- Rprintf("pdf_P = %10.5f\n", pdf_P);
- }
- }
- }// end of marginal likelihood computation
-}
-
-
-////////////////////////////////////////////
-// MCMCintervention fixed effect changes only
-////////////////////////////////////////////
-template <typename RNGTYPE>
-void MCMCintervention_fixed_impl(rng<RNGTYPE>& stream,
- const double m, const int intervention,
- const Matrix<>& Y, const Matrix<>& X,
- Matrix<>& beta, Matrix<>& Sigma, Matrix<>& P, Matrix<int>& s,
- Matrix<>& b0, Matrix<>& B0,
- const double c0, const double d0,
- const Matrix<>& A0,
- unsigned int burnin, unsigned int mcmc, unsigned int thin,
- unsigned int verbose, bool chib, bool ar,
- Matrix<>& beta_store, Matrix<>& Sigma_store,
- Matrix<>& P_store, Matrix<>& ps_store, Matrix<int>& s_store,
- double& logmarglike,
- Matrix<>& yhat_mat,
- Matrix<>& yerror_mat,
- Matrix<>& yfore_pred_mat,
- Matrix<>& yback_pred_mat,
- double acceptance)
-{
- // define constants and form cross-product matrices
- const int tot_iter = burnin + mcmc; //total iterations
- const int nstore = mcmc / thin; // number of draws to store
- const int n = Y.rows();
- const int ns = m + 1; // number of states
- const int k = X.cols();
- const Matrix<> B0inv = invpd(B0);
- double sigma2 = Sigma(0);
- double sigma = sqrt(sigma2);
-
- //MCMC loop
- unsigned int count = 0;
- unsigned int reject = 0;
- Matrix <> sum_e(ns, 1);
-
- for (int iter = 0; iter < tot_iter; ++iter){
- //////////////////////
- // 1. Sample beta and Sigma
- //////////////////////
- int beta_count = 0;
- Matrix<int> nstate(ns, 1); // contains total numbers of each state
- for (int j = 0; j<ns; ++j){
- for (int i = 0; i<n; ++i){
- if (s[i] == j + 1) { // since j starts from 0
- nstate[j] = nstate[j] + 1;
- }// end of if
- }// end of int i<n
- beta_count = beta_count + nstate[j];
-
- // BETA UPDATE
- Matrix<double> yj = Y((beta_count - nstate[j]), 0, (beta_count - 1), 0);
- Matrix<double> Xj = X((beta_count - nstate[j]), 0, (beta_count - 1), k-1);
- Matrix<double> tXj = ::t(Xj);
- Matrix<double> Bn = invpd(B0 + crossprod(Xj)/sigma2);
- Matrix<double> bn = Bn*(B0*b0 + tXj*yj/sigma2);
- if (ar == 1){
- Matrix<double> beta_can = stream.rmvnorm(bn, Bn);
- if (beta_can(1) > 1 | beta_can(1) < -1){
- // Rprintf("\n AR coefficient %10.5f is outside the stationary region! \n", beta_can(1));
- ++reject;
- }
- else{
- for (int kk = 0; kk<k; ++kk){
- beta(j,kk) = beta_can(kk);
- }
- }
- }
- else{
- beta(j,_) = stream.rmvnorm(bn, Bn);
- }
- Matrix<> yhat_j = Xj*::t(beta(j,_));
- Matrix<double> ej = yj - yhat_j;
- Matrix<double> sum_ej = t(ej)*ej;
- sum_e(j) = sum_ej[0];
- if (iter >= burnin && ((iter % thin)==0)){
- yhat_mat(count, (beta_count - nstate[j]), count, (beta_count - 1)) = yhat_j(_,0);
- yerror_mat(count, (beta_count - nstate[j]), count, (beta_count - 1)) = ej(_,0);
- }
-
- }// end of sampling beta and Sigma
-
- // SIGMA UPDATE
- double shape = (c0 + (double)n)/2;
- double scale =(d0 + sum(sum_e))/2;
-
- sigma2 = 1/stream.rgamma(shape, scale);
- sigma = sqrt(sigma2);
-
- //////////////////////
- // 2. 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) + (double)nstate[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);
- }
-
- //////////////////////
- // 3. Sample s
- //////////////////////
- Matrix<double> F(n, ns);
- Matrix<double> pr1(ns, 1);
- pr1[0] = 1;
- Matrix<double> py(ns, 1);
- Matrix<double> pstyt1(ns, 1);
- Matrix<double> ps = Matrix<double>(n, ns); // holder for state probabilities
-
- //
- // Forward sampling: update F matrix
- //
- for (int tt=0; tt<n ; ++tt){
- Matrix<double> mu = X(tt,_)*::t(beta); //k by 1 vector
- for (int j = 0; j< ns; ++j){
- py[j] = dnorm(Y[tt], mu[j], sigma);
- }
- if (tt==0) pstyt1 = pr1;
- else {
- pstyt1 = ::t(F(tt-1,_)*P); // make it an ns by 1 matrix
- }
- Matrix<double> unnorm_pstyt = pstyt1%py;
-
- /////////////////////////////////////////////////////////////////////
- // Prediction of future outcomes based on pre-intervention state
- /////////////////////////////////////////////////////////////////////
- if (tt==(intervention - 1)&&iter >= burnin && ((iter % thin)==0)){
- // Forward prediction
- Matrix <> yfore_pred(1, n);
- for (int ttt=tt; ttt<n ; ++ttt){
- int ss = s(tt-1); // 1 unit before the intervention
- mu = X(ttt,_)*::t(beta); //k by 1 vector
- yfore_pred(ttt) = stream.rnorm(mu[ss-1], sigma);
- }
- yfore_pred_mat(count, _) = yfore_pred(0, _);
- // backward prediction
- Matrix <> yback_pred(1, n);
- for (int ttt=tt; ttt>=0 ; --ttt){
- int ss = s(tt+1);
- mu = X(ttt,_)*::t(beta); //k by 1 vector
- yback_pred(ttt) = stream.rnorm(mu[ss-1], sigma);
- }
- yback_pred_mat(count, _) = yback_pred(0, _);
- }
-
- const Matrix<double> pstyt = unnorm_pstyt/sum(unnorm_pstyt); // pstyt = Pr(st|Yt)
- for (int j=0; j<ns ; ++j) F(tt,j) = pstyt(j);
-
- }// end of F matrix filtering
-
- //
- // Backward recursions: sample s using F and a transition matrix P
- //
- ps(n-1,_) = F(n-1,_); // we know last elements of ps and s
- s(n-1) = ns; // the last period is s(*n-1) but one down in C location
-
- Matrix<double> pstyn = Matrix<double>(ns, 1);
- double pone = 0.0;
- int tt = n-2;
- while (tt >= 0){
- int st = s(tt+1);
- Matrix<double> Pst_1 = ::t(P(_,st-1)); // prob of being at a previous state
- Matrix<double> unnorm_pstyn = F(tt,_)%Pst_1;
- pstyn = unnorm_pstyn/sum(unnorm_pstyn); // normalize into a prob. density
- if (st==1) s(tt) = 1; // If this is the first period, state should be 1.
- // Otherwise, draw a state from a discrete prob. distribution("pstyn")
- // using the inverse CDF method.
- else{
- pone = pstyn(st-2);
- if(stream.runif() < pone) s(tt) = st-1;// jump from tt-1 to tt
- else s(tt) = st;// stay
- }
- ps(tt,_) = pstyn;
- --tt;
- }// end of while loop
-
-
- // load draws into sample array
- if (iter >= burnin && ((iter % thin)==0)){
- Matrix<double> tbeta = ::t(beta); //transpose beta for R output
- for (int i=0; i<(ns*k); ++i)
- beta_store(count,i) = tbeta[i];// stored by the order of (11, 12, 13, 21, 22, 23)
- Sigma_store(count) = sigma2;
- 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,_); // add two matrices
- s_store(count,_) = s;
-
- ++count; // count when (iter % *thin)==0
-
- } // end of if(iter...)
-
-
- // print output to stdout
- if(verbose > 0 && iter % verbose == 0){
- Rprintf("\nMCMCintervention iteration %i of %i \n", (iter+1), tot_iter);
- if (ar == 1 ){
- double rejectionrate = (double)reject/(double)(iter+1);
- Rprintf("\n\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n");
- Rprintf("The acceptance rate was %3.5f", 1 - rejectionrate);
- Rprintf("\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n");
- }
- for (int j = 0;j<ns; ++j){
- Rprintf("\n nstate[j] %10.5f", static_cast<double>(nstate[j]));
- }
- Rprintf("\n beta \n");
- for (int i = 0; i<ns; ++i){
- for (int j = 0; j<k; ++j){
- Rprintf("%10.5f\n", beta(i, j));
- }
- }
- Rprintf("\n sigma^2 \n");
- Rprintf("%10.5f\n", sigma2);
- }
-
- }// end MCMC loop
-
- acceptance = 1 - (double)reject/(double)nstore;
-
- if(chib ==1){
- // 0.05 is used for an arbitrary threshold to avoid the computation of marg like
- // Later, the threshold should be set by users
- if(acceptance < .95){
- Rprintf("The acceptance rate is too low.\n");
- }
- else{
- Matrix<double> betast = meanc(beta_store); //meanc(beta_store)=(11, 12, 13, 21, 22, 23)
- Matrix<double, Row> beta_st(ns, k);
- for (int j = 0; j<ns*k; ++j){
- beta_st[j] = betast[j];
- }
-
- Matrix<double> Sigma_st(ns, 1);
- for (int j = 0; j<ns ; ++j){
- Sigma_st(j) = mean(Sigma_store);
- }
- Matrix<double> P_vec_st = meanc(P_store);
- const Matrix<double> P_st(ns, ns);
- for (int j = 0; j< ns*ns; ++j){
- P_st[j] = P_vec_st[j];
- }
-
- //////////////////////
- // 1. pdf.beta
- //////////////////////
- Matrix<double> density_beta(nstore, ns);
- for (int iter = 0; iter<nstore; ++iter){
-
- Matrix<int> nstate(ns, 1); // contains total numbers of each state
- int beta_count = 0;
-
- for (int j = 0; j<ns ; ++j){
- for (int i = 0; i<n; ++i){
- if (s_store(iter, i) == (j+1)) { // since j starts from 0
- nstate[j] = nstate[j] + 1;
- }// end of if
- }// end of int i<n
- beta_count = beta_count + nstate[j];
-
- // BETA UPDATE
- const Matrix<double> yj = Y((beta_count - nstate[j]), 0, (beta_count - 1), 0);
- const Matrix<double> Xj = X((beta_count - nstate[j]), 0, (beta_count - 1), k-1);
- const double precision = 1.0/Sigma_store(iter);
- const Matrix<double> XpX = crossprod(Xj);
- const Matrix<double> XpY = (::t(Xj)*yj);
- const Matrix<double> Bn = invpd(B0 + XpX*precision);
- const Matrix<double> bn = Bn*gaxpy(B0, b0, XpY*precision);
- if (k == 1){
- density_beta(iter, j) = dnorm(beta_st(j), bn(0), sqrt(Bn(0)));
- }
- else{
- density_beta(iter, j) = ::exp(lndmvn(::t(beta_st(j,_)), bn, Bn));
- }
- }// end of sampling beta and Sigma
-
- }// end of pdf.beta
-
- double pdf_beta = log(prod(meanc(density_beta)));
-
- //////////////////////
- // 2. pdf.Sigma|beta_st, S, P
- //////////////////////
- Matrix<double> density_Sigma(nstore, 1);
- for (int iter = 0; iter<nstore; ++iter){
-
- // STEP 2.1 S|y, beta.st, Sigma, P
- Matrix <double> Sout = gaussian_state_sampler(stream, m, Y, X, beta_st, Sigma, P);
- Matrix <double> s = Sout(_, 0);
-
- // STEP 2.2 Sigma|y, beta.st, S, P
- int beta_count = 0;
- Matrix<int> nstate(ns, 1);
- Matrix<> sum_e(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
- nstate[j] = nstate[j] + 1;
- }// end of if
- }// end of int i<n
- beta_count = beta_count + nstate[j];
-
- Matrix<double> yj = Y((beta_count - nstate[j]), 0, (beta_count - 1), 0);
- Matrix<double> Xj = X((beta_count - nstate[j]), 0, (beta_count - 1), k-1);
- Matrix<double> ej = yj - Xj*::t(beta_st(j,_));
- Matrix<double> sum_ej = ::t(ej)*ej;
- sum_e(j) = sum_ej(0);
- }// end of sampling beta and Sigma
- double shape = (c0 + (double)n)/2;
- double scale =(d0 + sum(sum_e))/2;
-
- sigma2 = stream.rigamma(shape, scale);
- density_Sigma(iter) = ::exp(lndinvgamma_pjh(Sigma_st(0), shape, scale));
-
- for (int j = 0; j <ns ; ++j){
- Sigma[j] = sigma2;// to use gaussian_state_sampler
- }
-
- // STEP 2.3 P|S
- 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) + (double)nstate[j] - 1;
- shape2 = A0(j,j+1) + 1; //
- P(j,j) = stream.rbeta(shape1, shape2);
- P(j,j+1) = 1 - P(j,j);
-
- // not to change gaussian_state_sampler, I use Sigma matrix
- Sigma[j] = sigma2;
- }
-
- }// end of pdf.Sigma
-
- double pdf_Sigma = log(prod(meanc(density_Sigma)));
-
- //////////////////////
- // 3. pdf.P|beta_st, Sigma_st, S
- //////////////////////
- Matrix<double> density_P(nstore, ns);
- for (int iter = 0; iter < nstore; ++iter){
-
- // STEP 2.1 S|y, beta.st, Sigma, P
- Matrix <double> Sout = gaussian_state_sampler(stream, m, Y, X, beta_st, Sigma_st, P);
- Matrix <double> s = Sout(_, 0);
-
- double shape1 = 0;
- double shape2 = 0;
- P(ns-1, ns-1) = 1; //no jump at the last state
-
- // compute addN
- Matrix <double> P_addN(ns, 1);
- for (int j = 0; j<ns; ++j){
- for (int i = 0; i<n; ++i){
- if (s(i) == (j+1)) { // since j starts from 0
- P_addN[j] = P_addN[j] + 1;
- }// end of if
- } // end of for i
- } // end of for j
-
- // evaluate P
- 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<double> F = Matrix<double>(n, ns);
- Matrix<double> like(n, 1);
- Matrix<double> pr1 = Matrix<double>(ns, 1);
- pr1[0] = 1;
- Matrix<double> py(ns, 1);
- Matrix<double> pstyt1(ns, 1);
-
- for (int t=0; t<n ; ++t){
- int yt = (int) Y[t];
- Matrix<double> mu = X(t,_)*::t(beta_st); //k by 1 vector
- for (int j = 0; j< ns; ++j){
- py[j] = dnorm(Y[t], mu[j], sqrt(Sigma_st[0]));
- }
- 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
-
- const double loglike = sum(log(like));
-
- //////////////////////
- // log prior ordinate
- //////////////////////
- Matrix<double> density_beta_prior(ns, 1);
- double density_Sigma_prior = 0.0;
- Matrix<double> density_P_prior(ns, 1);
- density_P[ns-1] = 1; //
- if (k == 1){
- for (int j=0; j<ns ; ++j){
- density_beta_prior[j] = log(dnorm(beta_st(j), b0(0), sqrt(B0inv(0))));
- }
- }
- else{
- for (int j=0; j<ns ; ++j){
- density_beta_prior[j] = lndmvn(::t(beta_st(j,_)), b0, B0inv);
- }
- }
-
- density_Sigma_prior = lndinvgamma_pjh(Sigma_st(0), c0/2, d0/2);
-
- for (int j =0; j< (ns-1); ++j){
- density_P_prior[j] = log(dbeta(P_st(j,j), A0(j,j), A0(j,j+1)));
- }
-
- // compute marginal likelihood
- double logprior = sum(density_beta_prior) + density_Sigma_prior + sum(density_P_prior);
- logmarglike = (loglike + logprior) - (pdf_beta + pdf_Sigma + pdf_P);
- if(verbose > 0){
- Rprintf("logmarglike = %10.5f\n", logmarglike);
- Rprintf("loglike = %10.5f\n", loglike);
- Rprintf("logprior = %10.5f\n", logprior);
- Rprintf("pdf_beta = %10.5f\n", pdf_beta);
- Rprintf("pdf_Sigma = %10.5f\n", pdf_Sigma);
- Rprintf("pdf_P = %10.5f\n", pdf_P);
- }
- }
- }// end of marginal likelihood computation
-}
-
-template <typename RNGTYPE>
-void MCMCintervention_impl(rng<RNGTYPE>& stream,
- const double m, const int intervention,
- const Matrix<>& Y, const Matrix<>& X,
- Matrix<>& beta, Matrix<>& Sigma, Matrix<>& P, Matrix<int>& s,
- Matrix<>& b0, Matrix<>& B0,
- const double c0, const double d0,
- const Matrix<>& A0,
- unsigned int burnin, unsigned int mcmc, unsigned int thin,
- unsigned int verbose, bool chib, bool ar,
- Matrix<>& beta_store, Matrix<>& Sigma_store,
- Matrix<>& P_store, Matrix<>& ps_store, Matrix<int>& s_store,
- double& logmarglike,
- Matrix<>& yhat_mat,
- Matrix<>& yerror_mat,
- Matrix<>& yfore_pred_mat,
- Matrix<>& yback_pred_mat,
- double acceptance)
-{
- // define constants and form cross-product matrices
- const int tot_iter = burnin + mcmc; //total iterations
- const int nstore = mcmc / thin; // number of draws to store
- const int n = Y.rows();
- const int ns = m + 1; // number of states
- const int k = X.cols();
- const Matrix<> B0inv = invpd(B0);
- Matrix<double> sigma(ns, 1);
-
- //MCMC loop
- unsigned int count = 0;
- int reject = 0;
- for (int iter = 0; iter < tot_iter; ++iter){
- //////////////////////
- // 1. Sample beta and Sigma
- //////////////////////
- int beta_count = 0;
- Matrix<int> nstate(ns, 1); // contains total numbers of each state
- for (int j = 0; j<ns; ++j){
- for (int i = 0; i<n; ++i){
- if (s[i] == j + 1) { // since j starts from 0
- nstate[j] = nstate[j] + 1;
- }// end of if
- }// end of int i<n
- beta_count = beta_count + nstate[j];
-
- // BETA UPDATE
- Matrix<double> yj = Y((beta_count - nstate[j]), 0, (beta_count - 1), 0);
- Matrix<double> Xj = X((beta_count - nstate[j]), 0, (beta_count - 1), k-1);
- Matrix<double> tXj = ::t(Xj);
- Matrix<double> Bn = invpd(B0 + crossprod(Xj)/Sigma[j]);
- Matrix<double> bn = Bn*(B0*b0 + tXj*yj/Sigma[j]);
- if (ar == 1){
- Matrix<double> beta_can = stream.rmvnorm(bn, Bn);
- if (beta_can(1) > 1 | beta_can(1) < -1){
- // Rprintf("\n AR coefficient %10.5f is outside the stationary region! \n", beta_can(1));
- ++reject;
- }
- else{
- for (int kk = 0; kk<k; ++kk){
- beta(j,kk) = beta_can(kk);
- }
- }
- }
- else{
- beta(j,_) = stream.rmvnorm(bn, Bn);
- }
-
- // SIGMA UPDATE
- double shape = (c0 + (double)nstate[j])/2;
- Matrix<> yhat_j = Xj*::t(beta(j,_));
- Matrix<double> ej = yj - yhat_j;
- Matrix<double> sum_ej = t(ej)*ej;
- double scale =(d0 + sum_ej[0])/2;
-
- Sigma[j] = 1/stream.rgamma(shape, scale);
- sigma(j) = sqrt(Sigma(j));
-
-
- if (iter >= burnin && ((iter % thin)==0)){
- yhat_mat(count, (beta_count - nstate[j]), count, (beta_count - 1)) = yhat_j(_,0);
- yerror_mat(count, (beta_count - nstate[j]), count, (beta_count - 1)) = ej(_,0);
- }
-
- }// end of sampling beta and Sigma
-
- //////////////////////
- // 2. 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) + (double)nstate[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);
- }
-
- //////////////////////
- // 3. Sample s
- //////////////////////
- Matrix<double> F(n, ns);
- Matrix<double> pr1(ns, 1);
- pr1[0] = 1;
- Matrix<double> py(ns, 1);
- Matrix<double> pstyt1(ns, 1);
- Matrix<double> ps = Matrix<double>(n, ns); // holder for state probabilities
-
- //
- // Forward sampling: update F matrix
- //
- for (int tt=0; tt<n ; ++tt){
- Matrix<double> mu = X(tt,_)*::t(beta); //k by 1 vector
- for (int j = 0; j< ns; ++j){
- py[j] = dnorm(Y[tt], mu[j], sigma[j]);
- }
- if (tt==0) pstyt1 = pr1;
- else {
- pstyt1 = ::t(F(tt-1,_)*P); // make it an ns by 1 matrix
- }
- Matrix<double> unnorm_pstyt = pstyt1%py;
-
- /////////////////////////////////////////////////////////////////////
- // Prediction of future outcomes based on pre-intervention state
- /////////////////////////////////////////////////////////////////////
- if (tt==(intervention - 1)&&iter >= burnin && ((iter % thin)==0)){
- // Forward prediction
- Matrix <> yfore_pred(1, n);
- for (int ttt=tt; ttt<n ; ++ttt){
- int ss = s(tt-1); // 1 unit before the intervention
- mu = X(ttt,_)*::t(beta); //k by 1 vector
- yfore_pred(ttt) = stream.rnorm(mu[ss-1], sigma[ss-1]);
- }
- yfore_pred_mat(count, _) = yfore_pred(0, _);
- // backward prediction
- Matrix <> yback_pred(1, n);
- for (int ttt=tt; ttt>=0 ; --ttt){
- int ss = s(tt+1);
- mu = X(ttt,_)*::t(beta); //k by 1 vector
- yback_pred(ttt) = stream.rnorm(mu[ss-1], sigma[ss-1]);
- }
- yback_pred_mat(count, _) = yback_pred(0, _);
- }
-
- const Matrix<double> pstyt = unnorm_pstyt/sum(unnorm_pstyt); // pstyt = Pr(st|Yt)
- for (int j=0; j<ns ; ++j) F(tt,j) = pstyt(j);
-
- }// end of F matrix filtering
-
- //
- // Backward recursions: sample s using F and a transition matrix P
- //
- ps(n-1,_) = F(n-1,_); // we know last elements of ps and s
- s(n-1) = ns; // the last period is s(*n-1) but one down in C location
-
- Matrix<double> pstyn = Matrix<double>(ns, 1);
- double pone = 0.0;
- int tt = n-2;
- while (tt >= 0){
- int st = s(tt+1);
- Matrix<double> Pst_1 = ::t(P(_,st-1)); // prob of being at a previous state
- Matrix<double> unnorm_pstyn = F(tt,_)%Pst_1;
- pstyn = unnorm_pstyn/sum(unnorm_pstyn); // normalize into a prob. density
- if (st==1) s(tt) = 1; // If this is the first period, state should be 1.
- // Otherwise, draw a state from a discrete prob. distribution("pstyn")
- // using the inverse CDF method.
- else{
- pone = pstyn(st-2);
- if(stream.runif() < pone) s(tt) = st-1;// jump from tt-1 to tt
- else s(tt) = st;// stay
- }
- ps(tt,_) = pstyn;
- --tt;
- }// end of while loop
-
-
- // load draws into sample array
- if (iter >= burnin && ((iter % thin)==0)){
- Matrix<double> tbeta = ::t(beta); //transpose beta for R output
- for (int i=0; i<(ns*k); ++i)
- beta_store(count,i) = tbeta[i];// stored by the order of (11, 12, 13, 21, 22, 23)
- for (int i=0; i<ns; ++i)
- Sigma_store(count,i) = Sigma[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,_); // add two matrices
- s_store(count,_) = s;
-
- ++count; // count when (iter % *thin)==0
-
- } // end of if(iter...)
-
-
- // print output to stdout
- if(verbose > 0 && iter % verbose == 0){
- Rprintf("\nMCMCintervention iteration %i of %i \n", (iter+1), tot_iter);
- if (ar == 1 ){
- double rejectionrate = (double)reject/(double)(iter+1);
- Rprintf("\n\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n");
- Rprintf("The acceptance rate was %3.5f", 1 - rejectionrate);
- Rprintf("\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n");
- }
- for (int j = 0;j<ns; ++j){
- Rprintf("\n nstate[j] %10.5f", static_cast<double>(nstate[j]));
- }
- Rprintf("\n beta \n");
- for (int i = 0; i<ns; ++i){
- for (int j = 0; j<k; ++j){
- Rprintf("%10.5f\n", beta(i, j));
- }
- }
- Rprintf("\n sigma^2 \n");
- for (int i = 0; i<ns; ++i){
- Rprintf("%10.5f\n", Sigma(i));
- }
- }
-
- }// end MCMC loop
-
- acceptance = 1 - (double)reject/(double)nstore;
-
- if(chib ==1){
- // 0.05 is used for an arbitrary threshold to avoid the computation of marg like
- // Later, the threshold should be set by users
- if(acceptance < .95){
- Rprintf("The acceptance rate is too low.\n");
- }
- else{
- Matrix<double> betast = meanc(beta_store);
- Matrix<double, Row> beta_st(ns, k);
- for (int j = 0; j<ns*k; ++j){
- beta_st[j] = betast[j];
- }
-
- Matrix<double> Sigma_st = meanc(Sigma_store);
- Matrix<double> sigma_st(ns, 1);
- for (int j = 0; j<ns ; ++j){
- sigma_st(j) = sqrt(Sigma_st(j));
- }
-
- Matrix<double> P_vec_st = meanc(P_store);
- const Matrix<double> P_st(ns, ns);
- for (int j = 0; j< ns*ns; ++j){
- P_st[j] = P_vec_st[j];
- }
-
- //////////////////////
- // 1. pdf.beta
- //////////////////////
- Matrix<double> density_beta(nstore, ns);
- for (int iter = 0; iter<nstore; ++iter){
-
- Matrix<int> nstate(ns, 1);
- int beta_count = 0;
-
- 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;
- }// end of if
- }// end of int i<n
- beta_count = beta_count + nstate[j];
-
- // BETA UPDATE
- const Matrix<double> yj = Y((beta_count - nstate[j]), 0, (beta_count - 1), 0);
- const Matrix<double> Xj = X((beta_count - nstate[j]), 0, (beta_count - 1), k-1);
- const double precision = 1.0/Sigma_store(iter, j);
- const Matrix<double> XpX = crossprod(Xj);
- const Matrix<double> XpY = (::t(Xj)*yj);
- const Matrix<double> Bn = invpd(B0 + XpX*precision);
- const Matrix<double> bn = Bn*gaxpy(B0, b0, XpY*precision);
- if (k == 1){
- density_beta(iter, j) = log(dnorm(beta_st(j), bn(0), sqrt(Bn(0))));
- }
- else{
- density_beta(iter, j) = lndmvn(::t(beta_st(j,_)), bn, Bn);
- }
- }// end of sampling beta and Sigma
-
- }// end of pdf.beta
-
- double pdf_beta = sum(meanc(density_beta));
-
- //////////////////////
- // 2. pdf.Sigma|beta_st, S, P
- //////////////////////
- Matrix<double> density_Sigma(nstore, ns);
- for (int iter = 0; iter<nstore; ++iter){
-
- // STEP 2.1 S|y, beta.st, Sigma, P
- Matrix <double> Sout = gaussian_state_sampler(stream, m, Y, X, beta_st, Sigma, P);
- Matrix <double> s = Sout(_, 0);
-
- // STEP 2.2 Sigma|y, beta.st, S, P
- int beta_count = 0;
- Matrix<int> nstate(ns, 1); // contains total numbers of each state
-
- for (int j = 0; j <ns ; ++j){
- for (int i = 0; i<n; ++i){
- if (s[i] == (j+1)) { // since j starts from 0
- nstate[j] = nstate[j] + 1;
- }// end of if
- }// end of int i<n
- beta_count = beta_count + nstate[j];
-
- Matrix<double> yj = Y((beta_count - nstate[j]), 0, (beta_count - 1), 0);
- Matrix<double> Xj = X((beta_count - nstate[j]), 0, (beta_count - 1), k-1);
- Matrix<double> ej = yj - Xj*::t(beta_st(j,_));
- Matrix<double> sum_ej = ::t(ej)*ej;
- double shape = (c0 + (double)nstate[j])/2;
- double scale =(d0 + sum_ej[0])/2;
-
- Sigma[j] = stream.rigamma(shape, scale);
- density_Sigma(iter, j) = ::exp(lndinvgamma_pjh(Sigma_st[j], shape, scale));
-
- }// end of sampling beta and Sigma
-
- // STEP 2.3 P|S
- 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) + (double)nstate[j] - 1;
- shape2 = A0(j,j+1) + 1; //
- P(j,j) = stream.rbeta(shape1, shape2);
- P(j,j+1) = 1 - P(j,j);
- }
-
- }// end of pdf.Sigma
-
- double pdf_Sigma = log(prod(meanc(density_Sigma)));
-
- //////////////////////
- // 3. pdf.P|beta_st, Sigma_st, S
- //////////////////////
- Matrix<double> density_P(nstore, ns);
- for (int iter = 0; iter < nstore; ++iter){
-
- // STEP 2.1 S|y, beta.st, Sigma, P
- Matrix <double> Sout = gaussian_state_sampler(stream, m, Y, X, beta_st, Sigma_st, P);
- Matrix <double> s = Sout(_, 0);
-
- double shape1 = 0;
- double shape2 = 0;
- P(ns-1, ns-1) = 1;
-
- // compute addN
- Matrix <double> P_addN(ns, 1);
- for (int j = 0; j<ns; ++j){
- for (int i = 0; i<n; ++i){
- if (s(i) == (j+1)) {
- P_addN[j] = P_addN[j] + 1;
- }// end of if
- } // end of for i
- } // end of for j
-
- // evaluate P
- 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<double> F = Matrix<double>(n, ns);
- Matrix<double> like(n, 1);
- Matrix<double> pr1 = Matrix<double>(ns, 1);
- pr1[0] = 1;
- Matrix<double> py(ns, 1);
- Matrix<double> pstyt1(ns, 1);
-
- for (int t=0; t<n ; ++t){
- int yt = (int) Y[t];
- Matrix<double> mu = X(t,_)*::t(beta_st);
- for (int j = 0; j< ns; ++j){
- py[j] = dnorm(Y[t], mu[j], sigma_st[j]);
- }
- if (t==0) pstyt1 = pr1;
- else {
- pstyt1 = ::t(F(t-1,_)*P_st);
- }
- Matrix<double> unnorm_pstyt = pstyt1%py;
- Matrix<double> pstyt = unnorm_pstyt/sum(unnorm_pstyt);
- for (int j=0; j<ns ; ++j){
- F(t,j) = pstyt(j);
- }
- like[t] = sum(unnorm_pstyt);
- }// end of likelihood computation
-
- const double loglike = sum(log(like));
-
- //////////////////////
- // log prior ordinate
- //////////////////////
- Matrix<double> density_beta_prior(ns, 1);
- Matrix<double> density_Sigma_prior(ns, 1);
- Matrix<double> density_P_prior(ns, 1);
- density_P[ns-1] = 1; //
- if (k == 1){
- for (int j=0; j<ns ; ++j){
- density_beta_prior[j] = log(dnorm(beta_st(j), b0(0), sqrt(B0inv(0))));
- }
- }
- else{
- for (int j=0; j<ns ; ++j){
- density_beta_prior[j] = lndmvn(::t(beta_st(j,_)), b0, B0inv);
- }
- }
-
- for (int j=0; j<ns ; ++j){
- density_Sigma_prior[j] = lndinvgamma_pjh(Sigma_st[j], c0/2, d0/2);
- }
-
- for (int j =0; j< (ns-1); ++j){
- density_P_prior[j] = log(dbeta(P_st(j,j), A0(j,j), A0(j,j+1)));
- }
-
- // compute marginal likelihood
- double logprior = sum(density_beta_prior) + sum(density_Sigma_prior) + sum(density_P_prior);
- logmarglike = (loglike + logprior) - (pdf_beta + pdf_Sigma + pdf_P);
- if(verbose > 0){
- Rprintf("logmarglike = %10.5f\n", logmarglike);
- Rprintf("loglike = %10.5f\n", loglike);
- Rprintf("logprior = %10.5f\n", logprior);
- Rprintf("pdf_beta = %10.5f\n", pdf_beta);
- Rprintf("pdf_Sigma = %10.5f\n", pdf_Sigma);
- Rprintf("pdf_P = %10.5f\n", pdf_P);
- }
- }// end of marginal likelihood computation
- }
-}
-
-////////////////////////////////////////////
-// Start MCMCinterventionpoint function
-///////////////////////////////////////////
-extern "C"{
- void MCMCintervention(double *accept, double *betaout, double *Sigmaout, double *Pout,
- double *psout, double *sout,
- double *yhatout, double *yerrorout, double *yforepredout, double *ybackpredout,
- const double *Ydata, const int *Yrow, const int *Ycol,
- const double *Xdata, const int *Xrow, const int *Xcol,
- const int *m, const int *intervention,
- const int *burnin, const int *mcmc, const int *thin, const int *verbose,
- const int *uselecuyer, const int *seedarray, const int *lecuyerstream,
- const double *betastart, const double *Sigmastart, const double *Pstart,
- const int *statestart,
- const double *a, const double *b,
- const double *b0data, const double *B0data,
- const double *c0, const double *d0,
- const double *A0data,
- double *logmarglikeholder, double *loglikeholder,
- const int *ar, const int *change, const int *chib){
-
- // pull together Matrix objects
- const Matrix <double> Y(*Yrow, *Ycol, Ydata);
- const Matrix <double> X(*Xrow, *Xcol, Xdata);
- 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 k = X.cols();
- const int ns = *m + 1; // number of states
-
- // generate starting values
- Matrix <> Sigma(ns, 1, Sigmastart);
- Matrix <> P(ns, ns, Pstart);
- Matrix <int> s(n, 1, statestart);
- Matrix <> b0(k, 1, b0data);
- Matrix <> B0(k, k, B0data);
- const Matrix <> A0(ns, ns, A0data);
- double logmarglike;
- double acceptance;
-
- // storage matrices
- Matrix<> P_store(nstore, ns*ns);
- Matrix<> ps_store(n, ns);
- Matrix<int> s_store(nstore, n);
-
- Matrix<> yhat_mat(nstore, n);
- Matrix<> yerror_mat(nstore, n);
- Matrix<> yfore_pred_mat(nstore, n);
- Matrix<> yback_pred_mat(nstore, n);
-
- if (*change == 1){
- // fixed effects change only
- Matrix <> beta(ns, k, betastart);
- Matrix<> beta_store(nstore, ns*k);
- Matrix<> Sigma_store(nstore, 1);
- MCMCPACK_PASSRNG2MODEL(MCMCintervention_fixed_impl,
- *m,
- *intervention,
- Y, X,
- beta, Sigma, P, s, b0, B0,
- *c0, *d0, A0,
- *burnin, *mcmc, *thin, *verbose,
- *chib, *ar,
- beta_store, Sigma_store,
- P_store, ps_store, s_store,
- logmarglike,
- yhat_mat,
- yerror_mat,
- yfore_pred_mat,
- yback_pred_mat,
- acceptance);
- for (int i = 0; i<(nstore*ns*k); ++i){
- betaout[i] = beta_store[i];
- }
- for (int i = 0; i<(nstore); ++i){
- Sigmaout[i] = Sigma_store[i];
- }
- }
- else if (*change == 2){
- // random effects change only
- Matrix<> beta(k, 1, betastart);
- Matrix<> beta_store(nstore, k);
- Matrix<> Sigma_store(nstore, ns);
- MCMCPACK_PASSRNG2MODEL(MCMCintervention_random_impl,
- *m,
- *intervention,
- Y, X,
- beta, Sigma, P, s, b0, B0,
- *c0, *d0, A0,
- *burnin, *mcmc, *thin, *verbose,
- *chib, *ar,
- beta_store, Sigma_store,
- P_store, ps_store, s_store,
- logmarglike,
- yhat_mat,
- yerror_mat,
- yfore_pred_mat,
- yback_pred_mat) ;
- for (int i = 0; i<(nstore*k); ++i){
- betaout[i] = beta_store[i];
- }
- for (int i = 0; i<(nstore*ns); ++i){
- Sigmaout[i] = Sigma_store[i];
- }
- acceptance = 1;
- }
- else {
- Matrix <> beta(ns, k, betastart);
- Matrix<> beta_store(nstore, ns*k);
- Matrix<> Sigma_store(nstore, ns);
- MCMCPACK_PASSRNG2MODEL(MCMCintervention_impl,
- *m,
- *intervention,
- Y, X,
- beta, Sigma, P, s, b0, B0,
- *c0, *d0, A0,
- *burnin, *mcmc, *thin, *verbose,
- *chib, *ar,
- beta_store, Sigma_store,
- P_store, ps_store, s_store,
- logmarglike,
- yhat_mat,
- yerror_mat,
- yfore_pred_mat,
- yback_pred_mat,
- acceptance);
- // return output
- for (int i = 0; i<(nstore*ns*k); ++i){
- betaout[i] = beta_store[i];
- }
- for (int i = 0; i<(nstore*ns); ++i){
- Sigmaout[i] = Sigma_store[i];
- }
- }
-
- logmarglikeholder[0] = logmarglike;
- accept[0] = acceptance;
-
- 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];
- yhatout[i] = yhat_mat[i];
- yerrorout[i] = yerror_mat[i];
- yforepredout[i] = yfore_pred_mat[i];
- ybackpredout[i] = yback_pred_mat[i];
- }
-
- }// end of MCMCpoissonChange
-}//end extern "C"
-
-#endif
-
-
-
-
diff --git a/src/MCMCregressChange.cc b/src/MCMCregressChange.cc
index 6174ef2..7642cea 100755
--- a/src/MCMCregressChange.cc
+++ b/src/MCMCregressChange.cc
@@ -370,7 +370,7 @@ void MCMCregressChange_impl(rng<RNGTYPE>& stream,
Matrix<double> pstyt1(ns, 1);
for (int t=0; t<n ; ++t){
- int yt = (int) Y[t];
+ //int yt = (int) Y[t];
Matrix<double> mu = X(t,_)*::t(beta_st); //k by 1 vector
for (int j = 0; j< ns; ++j){
py[j] = dnorm(Y[t], mu[j], sqrt(Sigma_st[j]));
@@ -452,7 +452,7 @@ extern "C"{
// pull together Matrix objects
const Matrix <double> Y(*Yrow, *Ycol, Ydata);
const Matrix <double> X(*Xrow, *Xcol, Xdata);
- const unsigned int tot_iter = *burnin + *mcmc; //total iterations
+ // 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 k = X.cols();
diff --git a/src/Makevars b/src/Makevars.win
similarity index 100%
rename from src/Makevars
rename to src/Makevars.win
--
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