[r-cran-mcmcpack] 62/90: Imported Upstream version 1.2-2
Andreas Tille
tille at debian.org
Fri Dec 16 09:07:48 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 8c0a04599cd81db0bba984771480914baf73978a
Author: Andreas Tille <tille at debian.org>
Date: Fri Dec 16 08:07:21 2016 +0100
Imported Upstream version 1.2-2
---
DESCRIPTION | 10 +-
HISTORY | 20 +
LICENSE | 6 +-
MD5 | 61 +-
NAMESPACE | 6 +-
R/HMMpanelFE.R | 168 +++
R/HMMpanelRE.R | 223 +++
R/MCMCbinaryChange.R | 0
R/MCMCresidualBreakAnalysis.R | 128 ++
R/btsutil.R | 30 +
R/distn.R | 2 +-
R/testpanelGroupBreak.R | 173 +++
R/testpanelSubjectBreak.R | 103 ++
R/zzz.R | 16 +-
data/PErisk.rda | Bin 4653 -> 1966 bytes
data/SupremeCourt.rda | Bin 3896 -> 470 bytes
man/HMMpanelFE.Rd | 205 +++
man/HMMpanelRE.Rd | 226 ++++
man/MCMCbinaryChange.Rd | 12 +-
man/MCMCoprobitChange.Rd | 13 +-
man/MCMCpoissonChange.Rd | 3 +-
man/MCMCprobitChange.Rd | 13 +-
man/MCMCregress.Rd | 2 +-
...MCMCregress.Rd => MCMCresidualBreakAnalysis.Rd} | 142 +-
man/make.breaklist.Rd | 30 +
man/testpanelGroupBreak.Rd | 179 +++
man/testpanelSubjectBreak.Rd | 179 +++
src/HMMmultivariateGaussian.cc | 996 ++++++++++++++
src/HMMpanelFE.cc | 434 ++++++
src/HMMpanelRE.cc | 1413 ++++++++++++++++++++
src/MCMCirtKdRob.cc | 24 +-
src/MCMCregress.cc | 10 +-
src/MCMCresidualBreakAnalysis.cc | 517 +++++++
src/algorithm.h | 5 +
src/datablock.h | 75 +-
src/error.h | 21 +-
src/matrix.h | 15 +-
src/optimize.h | 14 +-
src/rng.h | 19 +-
src/stat.h | 4 +-
40 files changed, 5325 insertions(+), 172 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index aa4033b..1dc4999 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
Package: MCMCpack
-Version: 1.1-4
-Date: 2011-9-21
+Version: 1.2-2
+Date: 2012-3-2
Title: Markov chain Monte Carlo (MCMC) Package
Author: Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park
Maintainer: Jong Hee Park <jhp at uchicago.edu>
@@ -15,9 +15,9 @@ Description: This package contains functions to perform Bayesian
and pseudo-random number generators for statistical
distributions, a general purpose Metropolis sampling algorithm,
and tools for visualization.
-License: GPL-2
+License: GPL-3
SystemRequirements: gcc (>= 4.0)
URL: http://mcmcpack.wustl.edu
-Packaged: 2011-09-26 13:42:11 UTC; jongheepark
+Packaged: 2012-03-02 21:27:17 UTC; adm
Repository: CRAN
-Date/Publication: 2011-09-26 14:23:46
+Date/Publication: 2012-03-03 06:49:19
diff --git a/HISTORY b/HISTORY
index 8ec5fba..795db05 100644
--- a/HISTORY
+++ b/HISTORY
@@ -1,6 +1,26 @@
//
// Changes and Bug Fixes
//
+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()
diff --git a/LICENSE b/LICENSE
index 6f9c100..0ea20ef 100644
--- a/LICENSE
+++ b/LICENSE
@@ -3,9 +3,9 @@ Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
and Jong Hee Park
-This program is free software; you can redistribute it and/or modify
+This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
-the Free Software Foundation; either version 2 of the License, or
+the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
@@ -20,8 +20,8 @@ Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
Please contact:
Andrew D. Martin, Ph.D.
-Professor and Chair, Political Science, Arts & Sciences
Professor and CERL Director, School of Law
+Professor, Political Science, Arts & Sciences
Washington University in St. Louis
(314) 935-5863 (Office)
diff --git a/MD5 b/MD5
index ae0a22f..130294a 100644
--- a/MD5
+++ b/MD5
@@ -1,8 +1,10 @@
-c1f2fefa751f851fe593cc8bfed0c466 *DESCRIPTION
-510bf0941463c026351c1697e7516b8b *HISTORY
-c795674882accfc42204f5c71cb40180 *LICENSE
-3910bb9b28917e4e2473a2c83c4d24fa *NAMESPACE
+dd3a6642c422701839fb4c4b6fc003dc *DESCRIPTION
+ed15423ad364f55ae6092d7ab7db5683 *HISTORY
+45b74f9b9cce80e9215266776a1be85c *LICENSE
+4851e0ecf23d44b58fcd53e597210689 *NAMESPACE
1eaf24e8d468ac0b624a9aa90cded306 *R/BayesFactors.R
+9c52f72ee21a71bfbfd312cb56233125 *R/HMMpanelFE.R
+ca0b3fff125085e13069fae3f237e5ed *R/HMMpanelRE.R
626c7d968d4f56fe081b595cd59d4976 *R/MCMCSVDreg.R
c93b534213a07cec26b71eaa48f02295 *R/MCMCbinaryChange.R
2736b2f4759f92e4e1d49d92fd0811e8 *R/MCMCdynamicEI.R
@@ -32,33 +34,38 @@ b98bdefca76d3316c664bcb2ba2bcc38 *R/MCMCprobit.R
31fdcb514ab53ba258a581536fcbda60 *R/MCMCprobitChange.R
e27dcf535f23c25018bb778b61640f22 *R/MCMCquantreg.R
4b0822eb63ca0687af1c99ddc8f9c01d *R/MCMCregress.R
+510f86a50f94be5b4bbd76cd361025f4 *R/MCMCresidualBreakAnalysis.R
720834bbb59d6c4ce7c610acd22558be *R/MCMCtobit.R
4d2218bcdb4715ac9d791f26d5469f9b *R/MCmodels.R
3a98b4d919a265b46c29d14310bb90d4 *R/SSVSquantreg.R
555c87acdcc2f8572b14151e29022e20 *R/SSVSquantregsummary.R
fdfcfe9909cadaf15c47bbb604be2257 *R/automate.R
-188a57ae832e3cc3bd447271d03a6996 *R/btsutil.R
-5aca34f825d2ec38f932094b6bbe601b *R/distn.R
+9eb132545aaf170132109e198464b69d *R/btsutil.R
+56829ce5340d26f9488b8b9894ff5cd0 *R/distn.R
18daa9fee800c7d8d5827405738531e2 *R/hidden-hmodels.R
0f904468b8ce8fc5c3869f87d11da905 *R/hidden.R
f799711e5d3fb7140cbc0995360339b3 *R/procrust.R
2e00b2b1593d22b5d1c579e92b294106 *R/scythe.R
+4de1106e38b8ca3b12396ec0d84bba4c *R/testpanelGroupBreak.R
+cb8aeff19ebf8c85748bad1a57484e7a *R/testpanelSubjectBreak.R
c67ac5238de55cdc8c49d017e0226ed9 *R/tomog.R
6377da3453982f578df4a886f673122c *R/utility.R
-a4b896c744c34c6cece2779d541b538d *R/zzz.R
+7e2021def8a88edb41fe78d0a0ae642c *R/zzz.R
e42f762c27e26e0ace372c03c47ac4e5 *README
46e955baeb55718a380a11fb34d2ea67 *cleanup
e09d1c5ac6250f832af90566a9e565fb *configure
21e1a006fb388ed39daf16e6d2ec9549 *configure.ac
e17202d22ce5ab5ab7b01e3b248a4008 *data/Nethvote.rda
-bfb9323ec59f28654e261f570bd1fb9a *data/PErisk.rda
+d03a27617f7933fbbb47732af740b129 *data/PErisk.rda
93e636ac58b966f9a74315f51be07d25 *data/Rehnquist.rda
accb7b5b46207a02e49de2702a6faff4 *data/Senate.rda
-f5d5cc1cd9f9f1699e249946f8767179 *data/SupremeCourt.rda
+2b95cdf5122fe0f86191ba668da542dc *data/SupremeCourt.rda
a5cfcf73e21c03eeaf4f3baa5c492c14 *inst/CITATION
2b42855343b5f0c51c7e9411aef57509 *man/BayesFactor.Rd
+872365fa451a1717a950a8102453eb8e *man/HMMpanelFE.Rd
+2e322e903df0d9c20513670710537995 *man/HMMpanelRE.Rd
293ff17f852d0e60dd1aac1f419353d2 *man/MCMCSVDreg.Rd
-d479bf1112276aef579b481f05eb3a3f *man/MCMCbinaryChange.Rd
+9af751257fcacc3a57543e4ae6a5e6cc *man/MCMCbinaryChange.Rd
5078ab3b9d86c7165909791acf1fc886 *man/MCMCdynamicEI.Rd
102a94ced183cf023da7696f948bb82f *man/MCMCdynamicIRT1d.Rd
29e89e27c6b11794c5dc6f50618375e0 *man/MCMCfactanal.Rd
@@ -76,14 +83,15 @@ b894aa07a0e68f9ffdeb3b09bab0aee2 *man/MCMChpoisson.Rd
751e56edabc12da2d8246f7c34cd7613 *man/MCMCmixfactanal.Rd
b7dfcf274126008379b0ee903f134b25 *man/MCMCmnl.Rd
cf07e2c9f307215c8fe841f5381b62f8 *man/MCMCoprobit.Rd
-ac737a1e0cc404c1ac8dc91e2e870725 *man/MCMCoprobitChange.Rd
+1556620313c5a4ca175fbe0b68f84161 *man/MCMCoprobitChange.Rd
a491294e97dd0b7aa5d0de525542ea8f *man/MCMCordfactanal.Rd
03f2e93e6876a4f0e87b3469bfe38d2f *man/MCMCpoisson.Rd
-3d38651d02ff76f5d7430bdccf56c1d2 *man/MCMCpoissonChange.Rd
+c80d0b7b63dbaf0df8e23d78dd9e65cc *man/MCMCpoissonChange.Rd
1917e296040b80c92525e5cb8ad02e71 *man/MCMCprobit.Rd
-0b79ed083594613481da272097d0d359 *man/MCMCprobitChange.Rd
+c3631c612ccc2d8975a7e164c8506331 *man/MCMCprobitChange.Rd
2f9caf8983e405f912eb5b25f29784eb *man/MCMCquantreg.Rd
-df339d8d9dce904214b8d3d120e88182 *man/MCMCregress.Rd
+e1852779dbd3598219f162860d60e8f5 *man/MCMCregress.Rd
+02144cc6cca72312f4eafb013cfe23c3 *man/MCMCresidualBreakAnalysis.Rd
c314e36afdd60449cf4f8ec08e00f80d *man/MCMCtobit.Rd
a96eacae6b64827e40505661fb636cd4 *man/MCbinomialbeta.Rd
193401fb83df78594bb965c7bbac968f *man/MCmultinomdirichlet.Rd
@@ -103,18 +111,24 @@ e6849fc07eb2dd45dcc07300e7cd149d *man/Senate.Rd
5203f2019ccb8b2dbf032c1689335ec3 *man/dtomog.Rd
b6f23a0e12fa9bb16f571261b265921f *man/invgamma.Rd
0051db062affc35be76c3ef5591dac82 *man/iwishart.Rd
+a5378f49cfb5f7ae58439b828e10fb84 *man/make.breaklist.Rd
629023d17bc9299449eec0b74a73340d *man/mptable.Rd
96ab1db8e6a5fb7c6640f06dd95a162c *man/noncenhypergeom.Rd
170bc2168de62e08f7157a8cbc635e06 *man/plotChangepoint.Rd
c64abc923f2c6f56a167a815c881ffb1 *man/plotState.Rd
6a88d877e88b5abfd7a66d659c9b4e6a *man/procrust.Rd
41feb17dda3b00cffa2f6176ba584b74 *man/readscythe.Rd
+f151473ede68e2cff21c06afd3060373 *man/testpanelGroupBreak.Rd
+8a3c03825fcdcdf22359faa43b330574 *man/testpanelSubjectBreak.Rd
254168bb5258fd00a7b766fb7d3bbfd9 *man/tomog.Rd
f5ba8eb700c67f188505a4585dd39f6b *man/topmodels.Rd
a4c560360ba912bb342a429e434bcc96 *man/vech.Rd
8699162f8e39a03d96dd633866f8c4ee *man/wishart.Rd
9d56915bf90694106df6586cec1bada0 *man/writescythe.Rd
ee8f580a79a520bbdf458410b49aae2a *man/xpnd.Rd
+103dadac598fcb2a6e20af0cfa7b3b69 *src/HMMmultivariateGaussian.cc
+16ec2efa868e3f0339d32bfc4917b0cd *src/HMMpanelFE.cc
+37113b38bf4693823d6a05cc8e3d2e99 *src/HMMpanelRE.cc
1c68060796310df44a049e66f10b14f7 *src/MCMCSVDreg.cc
6474ba3567cea0fc10d61bbdc35c649a *src/MCMCbinaryChange.cc
c6353855b63e402d7e60d248b22e6f50 *src/MCMCdynamicEI.cc
@@ -130,7 +144,7 @@ ea2d9f60541c419765dde67824ca9139 *src/MCMChregress.cc
fa7404ccc4ce053d13b2e0b9a86819b8 *src/MCMCirt1d.cc
abdf98660af4bed11d89e04074d9f1d1 *src/MCMCirtHier1d.cc
ac9b6d3f588b46cf009816bf028ee5ed *src/MCMCirtKdHet.cc
-01541fda7e59e9653da81371de726eb6 *src/MCMCirtKdRob.cc
+b6c0ff1884472c0e297b7d6770d61ca0 *src/MCMCirtKdRob.cc
82710c136632704b70c71cc8ee5dca5c *src/MCMClogit.cc
8e8ab807995eedcd108f11ad639475a9 *src/MCMClogituserprior.cc
d886c33fcacd18ba4435f4171159824f *src/MCMCmetrop1R.cc
@@ -147,29 +161,30 @@ a78375184abf7c97bf7fbf3f7614cad3 *src/MCMCoprobit.cc
8a2f4a79c6bc951310074f5757041424 *src/MCMCprobitChange.cc
a6a01c3906a8b92b9f1a12cfa92a1a63 *src/MCMCprobitres.cc
806af936660056856d4b908f8fb0cfa4 *src/MCMCquantreg.cc
-32d7032c2fda6bdd5a8177fd9858ffe4 *src/MCMCregress.cc
+aedab719c419ce9d80f534e967d2648d *src/MCMCregress.cc
+1cab26aaeea8a8b701325f41aadff24c *src/MCMCresidualBreakAnalysis.cc
ef5566f824887d9bbef5b76f1b64d839 *src/MCMCrng.h
f28dc81e54a6ca54a2202a0d6f52ed40 *src/MCMCtobit.cc
2add70b59d5e52089b08311d62f5011e *src/Makevars
30978b84a8789fd5a22a0c606b94c02f *src/Makevars.in
6ba1ab7daac8f7149663686fd7abfd53 *src/SSVSquantreg.cc
-514000a986da5e497a7474e17379c626 *src/algorithm.h
-da3e44a3ae5138017ebd112e8e7af5c4 *src/datablock.h
+beed5b98c04a5ffe412bafc994ebbc47 *src/algorithm.h
+70554a1bc2942623afb4f305d5b590ca *src/datablock.h
d489c387b86b46548f3f149fa6778734 *src/defs.h
2acfadfc5d2ac2a532dc37e9c1f694b1 *src/distributions.h
-8e2a5e1e884ef5dd8486920ce3646f30 *src/error.h
+60f270724a096968cc222cfcd1e6d6ea *src/error.h
176c95d1fc3f3f7073c5c49df98a3abb *src/ide.h
d9bbcd86bcee5ba918f3d368bd270712 *src/la.h
ef7b879e29de6e8e11527700115686d6 *src/lapack.h
cfc7f8d03a981a3453f81f2bd67d02c5 *src/lecuyer.cc
12d602243cc945e693e56174f31874a9 *src/lecuyer.h
-9ccd363601ec9907b2b6226355445ada *src/matrix.h
+56dab3f0ca5ab44b5038cf8c271992fe *src/matrix.h
7ab16dbb6757a5de37dc7de9867a3c10 *src/matrix_bidirectional_iterator.h
f924b218cd6eaed80846ed0cba534b03 *src/matrix_forward_iterator.h
3dc8b8e303daa7ffdf3e2d29b2fa43d9 *src/matrix_random_access_iterator.h
b0f0c5a5c5197df32d18b77af2b85957 *src/mersenne.h
-563652e8eec5ba3bd78ad6951a3ca27a *src/optimize.h
-1df60939513617beba99fe82eb7374c9 *src/rng.h
+b49f8f2678a9436166f96b70e711c83b *src/optimize.h
+120974a25cd887e85fbc748b3225eff4 *src/rng.h
b48baf77b93da895a04154ef15c8683f *src/rtmvnorm.h
8aa95ee02342c8863d2d03b0e722aa4a *src/smath.h
-c3770569055925fec0945088672f5b61 *src/stat.h
+6f0c8e88a1aa10a5198551f9856ad516 *src/stat.h
diff --git a/NAMESPACE b/NAMESPACE
index 6ce78ee..2910a27 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -12,6 +12,8 @@ export(
dnoncenhypergeom,
dtomogplot,
dwish,
+ HMMpanelFE,
+ HMMpanelRE,
MCbinomialbeta,
MCpoissongamma,
MCnormalnormal,
@@ -20,7 +22,6 @@ export(
MCMCdynamicEI,
MCMCdynamicIRT1d,
MCMCfactanal,
-## MCMChierBetaBinom,
MCMChierEI,
MCMChlogit,
MCMChpoisson,
@@ -45,6 +46,7 @@ export(
MCMCregress,
MCMCSVDreg,
MCMCtobit,
+ make.breaklist,
mptable,
plotState,
plotChangepoint,
@@ -57,6 +59,8 @@ export(
rnoncenhypergeom,
rwish,
SSVSquantreg,
+ testpanelGroupBreak,
+ testpanelSubjectBreak,
tomogplot,
topmodels,
vech,
diff --git a/R/HMMpanelFE.R b/R/HMMpanelFE.R
new file mode 100644
index 0000000..f014f14
--- /dev/null
+++ b/R/HMMpanelFE.R
@@ -0,0 +1,168 @@
+####################################################################
+## HMM Gaussian Panel with Time Varying Intercepts
+##
+## first written by Jong Hee Park on 02/2009
+## revised for MCMCpack inclusion on 09/2011
+######################################################################
+"HMMpanelFE" <-
+ function(subject.id, y, X, m,
+ mcmc=1000, burnin=1000, thin=1, verbose=0,
+ b0=0, B0=0.001, c0 = 0.001, d0 = 0.001, delta0=0, Delta0=0.001,
+ a = NULL, b = NULL, seed = NA, ...){
+ ## m should be a vector with a number of breaks for each group
+ ## id is a numeric list
+ ## p is a lag order
+ ## offset is the first time period from which each group starts
+
+ ## seeds
+ seeds <- form.seeds(seed)
+ lecuyer <- seeds[[1]]
+ seed.array <- seeds[[2]]
+ lecuyer.stream <- seeds[[3]]
+
+ ## Data
+ Y <- as.matrix(y);
+ X <- as.matrix(cbind(1, X)) ## the intercept is not reported
+ N <- nrow(Y);
+ K <- ncol(X)
+
+ mvn.prior <- form.mvn.prior(b0, B0, K)
+ b0 <- mvn.prior[[1]]
+ B0 <- mvn.prior[[2]]
+
+ mvn.prior <- form.mvn.prior(delta0, Delta0, 1)
+ delta0 <- mvn.prior[[1]]
+ Delta0 <- mvn.prior[[2]]
+
+ nstore <- mcmc/thin
+ nsubj <- length(unique(subject.id))
+
+ ## groupinfo matrix
+ ## col1: subj ID, col2: offset (first time C indexing), col3: #time periods
+ if (unique(subject.id)[1] != 1){
+ stop("subject.id should start 1!")
+ }
+ subject.offset <- c(0, which(diff(sort(subject.id))==1)[-nsubj])
+ ## col1: subj ID, col2: offset (C indexing), col3: #time periods in each subject
+ nsubject.vec <- rep(NA, nsubj)
+ for (i in 1:nsubj){
+ nsubject.vec[i] <- sum(subject.id==unique(subject.id)[i])
+ }
+ subject.groupinfo <- cbind(unique(subject.id), subject.offset, nsubject.vec)
+
+ ## maximum time length
+ ntime <- max(nsubject.vec)
+ m.max <- max(m)
+ m.min <- min(m)
+
+ ## prior inputs
+ P0data <- NULL
+ Pstart <- NULL
+ for (i in 1:nsubj){
+ if(m[i] == 0){
+ P0current <- 1
+ Pstartcurrent <- 1
+ }
+ else{
+ P0current <- trans.mat.prior(m=m[i], n=nsubject.vec[i], a=a, b=b)
+ Pstartcurrent <- trans.mat.prior(m=m[i], n=nsubject.vec[i], a=.9, b=.1)
+ }
+ P0data <- c(P0data, P0current)
+ Pstart <- c(Pstart, Pstartcurrent)
+ }
+
+ ## starting values
+ ols <- lm(Y~X-1)
+ beta.start <- coef(ols)
+ sigma2.start <- summary(ols)$sigma^2
+ deltastart <- NULL
+ Ytilde <- Y - X%*%beta.start
+ deltaformula <- Ytilde ~ 1 ## without intercept
+ for (i in 1:nsubj){
+ deltacurrent <- rep(as.vector(coef(lm(Ytilde ~ 1))), m[i] + 1)
+ deltastart <- c(deltastart, deltacurrent)
+ }
+
+ ## Storage
+ totalstates0 <- sum(m+1)
+ betadraws0 <- matrix(0, nstore, K)
+ deltadraws0 <- matrix(data=0, nstore, totalstates0)
+ sigmadraws0 <- matrix(data=0, nstore, totalstates0)
+ statedraws0 <- matrix(data=0, nstore, totalstates0)
+
+ ## call C++ code to draw sample
+ posterior <- .C("HMMpanelFE",
+ deltadraws = as.double(deltadraws0),
+ sigmadraws = as.double(sigmadraws0),
+ statedraws = as.double(statedraws0),
+
+ betadraws = as.double(betadraws0),
+ betarow = as.integer(nrow(betadraws0)),
+ betacol = as.integer(ncol(betadraws0)),
+ totalstates = as.integer(totalstates0),
+
+ nsubj = as.integer(nsubj),
+ ntime = as.integer(ntime),
+ nobs = as.integer(N),
+ subjectid = as.integer(subject.id),
+
+ m = as.integer(m),
+ mmax = as.integer(m.max),
+ mmin = as.integer(m.min),
+
+ 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)),
+
+ 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),
+
+ betastartdata = as.double(beta.start),
+ sigma2start = as.double(sigma2.start),
+ deltastartdata = as.double(deltastart),
+ deltastartrow = as.integer(length(deltastart)),
+
+ b0data = as.double(b0),
+ B0data = as.double(B0),
+
+ delta0data = as.double(delta0),
+ Delta0data = as.double(Delta0),
+
+ c0 = as.double(c0),
+ d0 = as.double(d0),
+
+ P0data = as.double(P0data),
+ P0row = as.integer(length(P0data)),
+ Pstartdata = as.double(Pstart),
+
+ subject_groupinfodata = as.double(subject.groupinfo),
+ PACKAGE="MCMCpack")
+
+ ## pull together matrix and build MCMC object to return
+ betadraws <- matrix(posterior$betadraws,
+ posterior$betarow,
+ posterior$betacol)
+
+ sigma.samp <- as.mcmc(matrix(posterior$sigmadraws, nstore, totalstates0))
+ delta.samp <- as.mcmc(matrix(posterior$deltadraws, nstore, totalstates0))
+ state.samp <- as.mcmc(matrix(posterior$statedraws, nstore, totalstates0))
+
+ ## output <- mcmc(betadraws, start=burnin+1, end=burnin+mcmc, thin=thin)
+ output <- as.mcmc(betadraws[, -1])## drop the intercept
+ attr(output, "title") <- "HMMpanelFE Sample"
+ attr(output, "m") <- m
+ attr(output, "sigma") <- sigma.samp
+ attr(output, "state") <- state.samp
+ attr(output, "delta") <- delta.samp
+ return(output)
+ }
diff --git a/R/HMMpanelRE.R b/R/HMMpanelRE.R
new file mode 100644
index 0000000..0f4e3aa
--- /dev/null
+++ b/R/HMMpanelRE.R
@@ -0,0 +1,223 @@
+####################################################################
+## HMM Gaussian Panel Random Effects Model
+## y_it = x_it'*beta + w_it*bi + e_it,
+## e_it ~ N(0, sigma2)
+##
+## bi ~ N(0, D) : random effects coefficient
+## D ~ IW(r0, R0) : covariance matrix for multiple random effects
+## beta ~ N(b0, B0) : fixed effect coefficient
+## sigma2 ~ IG(c0/2, d0/2) : random error
+##
+## written by Jong Hee Park 03/2009
+## modified and integrated with other codes on 09/2011
+######################################################################
+
+"HMMpanelRE" <-
+ function(subject.id, time.id, y, X, W, m=1,
+ mcmc=1000, burnin=1000, thin=1, verbose=0,
+ b0=0, B0=0.001, c0 = 0.001, d0 = 0.001, r0, R0, a = NULL, b = NULL,
+ seed = NA, beta.start = NA, sigma2.start = NA, D.start= NA, P.start = NA,
+ marginal.likelihood = c("none", "Chib95"), ...){
+
+ cl <- match.call()
+ ## seeds
+ seeds <- form.seeds(seed)
+ lecuyer <- seeds[[1]]
+ seed.array <- seeds[[2]]
+ lecuyer.stream <- seeds[[3]]
+
+ ## Data
+ ns <- m + 1
+ Y <- as.matrix(y)
+ X <- as.matrix(X)
+ W <- as.matrix(W)
+ formula <- Y ~ X-1
+ ols <- lm(formula)
+
+ N <- nrow(Y)
+ K <- ncol(X)
+ Q <- ncol(W)
+
+ nobs <- nrow(X)
+
+ ## Sort Data based on time.id
+ oldTSCS <- cbind(time.id, subject.id, y, X, W)
+ newTSCS <- oldTSCS[order(oldTSCS[,1]),]
+ YT <- as.matrix(newTSCS[,3])
+ XT <- as.matrix(newTSCS[,4:(4+K-1)])
+ WT <- as.matrix(newTSCS[,(4+K):(4+K+Q-1)])
+
+ mvn.prior <- form.mvn.prior(b0, B0, K)
+ b0 <- mvn.prior[[1]]
+ B0 <- mvn.prior[[2]]
+
+ R0 <- as.matrix(R0)
+
+ nstore <- mcmc/thin
+ nsubj <- length(unique(subject.id))
+ if (unique(subject.id)[1] != 1){
+ stop("subject.id should start 1!")
+ }
+ ## subject.offset is the obs number from which a new subject unit starts
+ subject.offset <- c(0, which(diff(sort(subject.id))==1)[-nsubj])
+ ## col1: subj ID, col2: offset (C indexing), col3: #time periods in each subject
+ nsubject.vec <- rep(NA, nsubj)
+ for (i in 1:nsubj){
+ nsubject.vec[i] <- sum(subject.id==unique(subject.id)[i])
+ }
+ subject.groupinfo <- cbind(unique(subject.id), subject.offset, nsubject.vec)
+
+
+ ## time.groupinfo
+ ## col1: time ID, col2: offset (C indexing), col3: # subjects in each time
+ if(unique(time.id)[1] != 1){
+ time.id <- time.id - unique(time.id)[1] + 1
+ cat("time.id does not start from 1. So it is modified by subtracting the first unit of time.")
+ }
+ ntime <- max(nsubject.vec)## maximum time length
+ ntime.vec <- rep(NA, ntime)
+ for (i in 1:ntime){
+ ntime.vec[i] <- sum(time.id==unique(time.id)[i])
+ }
+ ## time.offset is the obs number from which a new time unit starts when we stack data by time.id
+ time.offset <- c(0, which(diff(sort(time.id))==1)[-ntime])
+ time.groupinfo <- cbind(unique(time.id), time.offset, ntime.vec)
+
+
+ ## prior inputs
+ if (m > 0){
+ P0 <- trans.mat.prior(m=m, n=ntime, a=a, b=b)
+ ## initial values
+ Pstart <- check.P(P.start, m, a=a, b=b)
+ }
+ else {
+ Pstart <- P0 <- matrix(1, 1, 1)
+ }
+ if (is.na(beta.start[1])) {
+ betastart <- coef(ols)
+ }
+ else{
+ betastart <- beta.start
+ }
+ if (is.na(sigma2.start[1])) {
+ sigma2start <- summary(ols)$sigma^2
+ }
+ else{
+ sigma2start <- sigma2.start
+ }
+
+ betadraws <- matrix(data=0, nstore, ns*K)
+ sigmadraws <- matrix(data=0, nstore, ns)
+ Ddraws <- matrix(data=0, nstore, ns*Q*Q)
+ psdraws <- matrix(data=0, ntime, ns)
+ sdraws <- matrix(data=0, nstore, ntime*ns)
+
+ ## 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
+ }
+
+ ## call C++ code to draw sample
+ posterior <- .C("HMMpanelRE",
+ betadata = as.double(betadraws),
+ betarow = as.integer(nrow(betadraws)),
+ betacol = as.integer(ncol(betadraws)),
+ sigmadata = as.double(sigmadraws),
+ Ddata = as.double(Ddraws),
+ psout = as.double(psdraws),
+ sout = as.double(sdraws),
+
+ nsubj = as.integer(nsubj), ntime = as.integer(ntime), m = as.integer(m),
+ nobs = as.integer(nobs),
+ subjectid = as.integer(subject.id),
+ timeid = as.integer(time.id),
+
+ 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)),
+ Wdata = as.double(W), Wrow = as.integer(nrow(W)), Wcol = as.integer(ncol(W)),
+ YTdata = as.double(YT), XTdata = as.double(XT), WTdata = as.double(WT),
+ 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),
+
+ betastartdata = as.double(betastart),
+ sigma2start = as.double(sigma2start),
+ Pstart = as.double(Pstart),
+
+ b0data = as.double(b0), B0data = as.double(B0),
+ c0 = as.double(c0), d0 = as.double(d0),
+ r0 = as.integer(r0), R0data = as.double(R0),
+
+ subject_groupinfodata = as.double(subject.groupinfo),
+ time_groupinfodata = as.double(time.groupinfo),
+
+ logmarglikeholder = as.double(0),
+ loglikeholder = as.double(0),
+ chib = as.integer(chib),
+ PACKAGE="MCMCpack"
+ )
+
+ ## pull together matrix and build MCMC object to return
+ beta.samp <- matrix(posterior$betadata,
+ posterior$betarow,
+ posterior$betacol)
+ ## stored by the order of (11, 12, 13, 21, 22, 23)
+ sigma.samp <- matrix(posterior$sigmadata,
+ posterior$betarow,
+ ns)
+ D.samp <- matrix(posterior$Ddata,
+ posterior$betarow,
+ Q*Q*ns)
+ xnames <- sapply(c(1:K), function(i){paste("beta", i, sep = "")})
+ Dnames <- sapply(c(1:(Q*Q)), function(i){paste("D", i, sep = "")})
+ if (m == 0){
+ output <- as.mcmc(cbind(beta.samp, sigma.samp, D.samp))
+ names <- c(xnames, "sigma2", Dnames)
+ varnames(output) <- as.list(names)
+ }
+ else{
+ output1 <- mcmc(data=beta.samp, start=burnin+1, end=burnin + mcmc, thin=thin)
+ output2 <- mcmc(data=sigma.samp, start=burnin+1, end=burnin + mcmc, thin=thin)
+ output3 <- mcmc(data=D.samp, start=burnin+1, end=burnin + mcmc, thin=thin)
+ varnames(output1) <- sapply(c(1:ns),
+ function(i){
+ paste(xnames, "_regime", i, sep = "")
+ })
+ varnames(output2) <- sapply(c(1:ns),
+ function(i){
+ paste("sigma2_regime", i, sep = "")
+ })
+ varnames(output3) <- sapply(c(1:ns),
+ function(i){
+ paste(Dnames, "_regime", i, sep = "")
+ })
+
+ output <- as.mcmc(cbind(output1, output2, output3))
+ ps.holder <- matrix(posterior$psout, ntime, ns)
+ s.holder <- matrix(posterior$sout, nstore, )
+ }
+
+ attr(output, "title") <- "HMMpanelRE Posterior Sample"
+ attr(output, "call") <- cl
+ attr(output, "y") <- y[1:ntime]
+ attr(output, "X") <- X[1:ntime, ]
+ attr(output, "m") <- m
+ attr(output, "nsubj") <- nsubj
+ attr(output, "ntime") <- ntime
+ if (m > 0){
+ attr(output, "s.store") <- s.holder
+ attr(output, "prob.state") <- ps.holder/nstore
+ }
+ attr(output, "logmarglike") <- posterior$logmarglikeholder
+ attr(output, "loglike") <- posterior$loglikeholder
+
+ return(output)
+ }
diff --git a/R/MCMCbinaryChange.R b/R/MCMCbinaryChange.R
old mode 100755
new mode 100644
diff --git a/R/MCMCresidualBreakAnalysis.R b/R/MCMCresidualBreakAnalysis.R
new file mode 100644
index 0000000..cfae683
--- /dev/null
+++ b/R/MCMCresidualBreakAnalysis.R
@@ -0,0 +1,128 @@
+#########################################################
+## residual break analysis
+## JHP 07/01/2007
+## JHP 03/03/2009
+#########################################################
+
+"MCMCresidualBreakAnalysis"<-
+ function(resid, m = 1,
+ b0 = 0, B0 = 0.001, c0 = 0.1, d0 = 0.1, 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.matrix(resid, ,1)
+ n <- nrow(y)
+ ns <- m + 1 # number of states
+
+ ## check iteration parameters
+ check.mcmc.parameters(burnin, mcmc, thin)
+ totiter <- mcmc + burnin
+ cl <- match.call()
+ nstore <- mcmc/thin
+
+ ## seeds
+ seeds <- form.seeds(seed)
+ lecuyer <- seeds[[1]]
+ seed.array <- seeds[[2]]
+ lecuyer.stream <- seeds[[3]]
+ if(!is.na(seed)) set.seed(seed)
+
+ ## get marginal likelihood argument
+ marginal.likelihood <- match.arg(marginal.likelihood)
+
+ ## following MCMCregress, set chib as binary
+ logmarglike <- loglike <- NULL
+ chib <- 0
+ if (marginal.likelihood == "Chib95"){
+ chib <- 1
+ }
+
+ ## initial values
+ if(m == 0){
+ output <- MCMCregress(y~1, mcmc=mcmc, burnin=burnin, verbose=verbose, thin=thin,
+ b0=b0, B0=B0, c0=c0, d0=d0, seed = seed,
+ beta.start=beta.start, marginal.likelihood = marginal.likelihood)
+ attr(output, "y") <- y
+ }
+ else{
+ ## prior
+ A0 <- trans.mat.prior(m=m, n=n, a=a, b=b)
+ Pstart <- check.P(P.start, m, a=a, b=b)
+ betastart <- rep(mean(y), ns)
+ Sigmastart <- rep(var(y), ns)
+ statestart <- sort(sample(1:ns, n, replace=T))
+
+ ## call C++ code to draw sample
+ posterior <- .C("MCMCresidualBreakAnalysis",
+ betaout = as.double(rep(0.0, nstore*ns)),
+ Sigmaout = as.double(rep(0.0, nstore*ns)),
+ psout = as.double(rep(0.0, n*ns)),
+
+ Ydata = as.double(y),
+ Yrow = as.integer(nrow(y)),
+ Ycol = as.integer(ncol(y)),
+
+ m = as.integer(m),
+ burnin = as.integer(burnin),
+ mcmc = as.integer(mcmc),
+ thin = as.integer(thin),
+ verbose = as.integer(verbose),
+
+ lecuyer=as.integer(lecuyer),
+ seedarray=as.integer(seed.array),
+ lecuyerstream=as.integer(lecuyer.stream),
+
+ 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),
+ 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, ns)
+ Sigma.holder <- matrix(posterior$Sigmaout, nstore, ns)
+ ps.holder <- matrix(posterior$psout, n, )
+
+ output1 <- mcmc(data=beta.holder, start=burnin+1, end=burnin + mcmc, thin=thin)
+ varnames(output1) <- sapply(c(1:ns),
+ function(i){
+ paste("mu_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("sigma2_regime", i, sep = "")
+ })
+ output <- as.mcmc(cbind(output1, output2))
+
+ attr(output, "title") <- "MCMCresidualBreakAnalysis Posterior Sample"
+ attr(output, "formula") <- formula
+ attr(output, "y") <- y
+ attr(output, "m") <- m
+ attr(output, "call") <- cl
+ attr(output, "prob.state") <- ps.holder/nstore
+ attr(output, "logmarglike") <- logmarglike
+ attr(output, "loglike") <- loglike
+ }
+ return(output)
+
+}## end of MCMC function
+
diff --git a/R/btsutil.R b/R/btsutil.R
index d04a58a..872857a 100644
--- a/R/btsutil.R
+++ b/R/btsutil.R
@@ -216,6 +216,7 @@
return(tau)
}
+
## beta starting values in MCMCpoissonChange()
"beta.change.start"<- function (beta.start, ns, k, formula, family, data){
## if a user does not specify beta.start, use a coefficient vector from mle
@@ -245,3 +246,32 @@
}
return(beta.start)
}
+
+
+
+## find the most reasonable model
+## by choosing a model with a larger break number
+## If the difference of natural log marginal likelihood is smaller than the threshold,
+## choose the model with a smaller break number
+## to avoid the over-identification problem
+"make.breaklist" <- function(BF, threshold=3){
+ eBF <- exp(BF)
+ N <- nrow(BF)
+ out <- rep(NA, N)
+ for (i in 1:N){
+ order.i <- order(eBF[i,], decreasing=TRUE)
+ if(eBF[i, order.i[1]] / eBF[i, order.i[2]] > threshold){
+ out[i] <- order.i[1]
+ }
+ else if (eBF[i, order.i[1]] / eBF[i, order.i[2]] < threshold & order.i[1]<=order.i[2]){
+ out[i] <- order.i[1]
+ }
+ else if (eBF[i, order.i[1]] / eBF[i, order.i[2]] < threshold & order.i[1]>order.i[2]){
+ out[i] <- order.i[2]
+ }
+ else{
+ cat("\n Error occurs at i = ", i)
+ }
+ }
+ return(out-1)
+}
diff --git a/R/distn.R b/R/distn.R
index 6539524..d268460 100644
--- a/R/distn.R
+++ b/R/distn.R
@@ -87,7 +87,7 @@
stop(message="W not square in dwish()\n\n")
}
if (!is.matrix(W))
- S <- matrix(W)
+ W <- matrix(W)
if (nrow(W) != ncol(W)){
stop(message="W not square in dwish()\n\n")
}
diff --git a/R/testpanelGroupBreak.R b/R/testpanelGroupBreak.R
new file mode 100644
index 0000000..2d06f9f
--- /dev/null
+++ b/R/testpanelGroupBreak.R
@@ -0,0 +1,173 @@
+####################################################################
+## test group-level breaks from panel residuals
+##
+## written by Jong Hee Park 03/2009
+## modified and integrated with other codes by JHP 07/2011
+######################################################################
+
+"testpanelGroupBreak" <-
+ function(subject.id, time.id, resid, m=1,
+ mcmc=1000, burnin=1000, thin=1, verbose=0,
+ b0, B0, c0, d0, a = NULL, b = NULL,
+ seed = NA, marginal.likelihood = c("none", "Chib95"), ...){
+ ## beta.start and sigma2.start are not arguments. OLS estimates will be used!
+ ## subject.id is a numeric list indicating the group number. It should start from 1.
+ ## time.id is a numeric list indicating the time, starting from 1.
+ ## subject.offset is the obs number from which a new subject unit starts
+ ## time.offset is the obs number from which a new time unit starts when we stack data by time.id
+
+ cl <- match.call()
+ ## seeds
+ seeds <- form.seeds(seed)
+ lecuyer <- seeds[[1]]
+ seed.array <- seeds[[2]]
+ lecuyer.stream <- seeds[[3]]
+
+ ## Data
+ ns <- m + 1
+ nobs <- length(subject.id)
+ newY <- matrix(resid, nobs, 1)
+ newX <- matrix(1, nobs, 1)
+ K <- 1
+
+ ## Sort Data based on time.id
+ oldTSCS <- cbind(time.id, subject.id, newY, newX)
+ newTSCS <- oldTSCS[order(oldTSCS[,1]),]
+ newYT <- as.matrix(newTSCS[,3])
+ newXT <- as.matrix(newTSCS[,4])
+ b0 <- as.matrix(b0)
+ B0 <- as.matrix(B0)
+
+ nstore <- mcmc/thin
+ nsubj <- length(unique(subject.id))
+ ## subject.groupinfo matrix
+ if (unique(subject.id)[1] != 1){
+ stop("subject.id should start 1!")
+ }
+ subject.offset <- c(0, which(diff(sort(subject.id))==1)[-nsubj])
+ nsubject.vec <- rep(NA, nsubj)
+ for (i in 1:nsubj){
+ nsubject.vec[i] <- sum(subject.id==unique(subject.id)[i])
+ }
+ subject.groupinfo <- cbind(unique(subject.id), subject.offset, nsubject.vec)
+
+ ## time.groupinfo
+ if(unique(time.id)[1] != 1){
+ time.id <- time.id - unique(time.id)[1] + 1
+ cat("time.id does not start from 1. So it is modified by subtracting the first unit of time.")
+ }
+ ntime <- max(nsubject.vec)## maximum time length
+ ntime.vec <- rep(NA, ntime)
+ for (i in 1:ntime){
+ ntime.vec[i] <- sum(time.id==unique(time.id)[i])
+ }
+ time.offset <- c(0, which(diff(sort(time.id))==1)[-ntime])
+ time.groupinfo <- cbind(unique(time.id), time.offset, ntime.vec)
+
+ ## prior inputs
+ if (m > 0){
+ P0 <- trans.mat.prior(m=m, n=ntime, a=a, b=b)
+ }
+ else {
+ P0 <- matrix(1, 1, 1)
+ }
+ betadraws <- matrix(data=0, nstore, ns*K)
+ sigmadraws <- matrix(data=0, nstore, ns)
+ psdraws <- matrix(data=0, ntime, ns)
+ ols <- lm(newY ~ newX - 1)
+ beta.start <- rep(coef(ols)[1], ns)
+ sigma2.start <- summary(ols)$sigma^2
+
+ ## 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
+ }
+
+ ## call C++ code to draw sample
+ posterior <- .C("HMMmultivariateGaussian",
+ betadata = as.double(betadraws),
+ betarow = as.integer(nrow(betadraws)),
+ betacol = as.integer(ncol(betadraws)),
+
+ sigmadata = as.double(sigmadraws),
+ psout = as.double(psdraws),
+ nsubj = as.integer(nsubj), ntime = as.integer(ntime),
+ m = as.integer(m),
+ nobs = as.integer(nobs),
+ subjectid = as.integer(subject.id),
+ timeid = as.integer(time.id),
+
+ Ydata = as.double(newY), Yrow = as.integer(nrow(newY)), Ycol = as.integer(ncol(newY)),
+ Xdata = as.double(newX), Xrow = as.integer(nrow(newX)), Xcol = as.integer(ncol(newX)),
+ YTdata = as.double(newYT), XTdata = as.double(newXT),
+ 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),
+
+ betastartdata = as.double(beta.start),
+ sigma2start = as.double(sigma2.start),
+
+ b0data = as.double(b0), B0data = as.double(B0),
+ c0 = as.double(c0), d0 = as.double(d0),
+ P0data = as.double(P0),
+ P0row = as.integer(nrow(P0)),
+ P0col = as.integer(ncol(P0)),
+
+ subject_groupinfodata = as.double(subject.groupinfo),
+ time_groupinfodata = as.double(time.groupinfo),
+
+ logmarglikeholder = as.double(0),
+ loglikeholder = as.double(0),
+ chib = as.integer(chib),
+ PACKAGE="MCMCpack"
+ )
+
+ ## pull together matrix and build MCMC object to return
+ beta.samp <- matrix(posterior$betadata,
+ posterior$betarow,
+ posterior$betacol)
+ ## stored by the order of (11, 12, 13, 21, 22, 23)
+ sigma.samp <- matrix(posterior$sigmadata,
+ posterior$betarow,
+ ns)
+ xnames <- sapply(c(1:K), function(i){paste("beta", i, sep = "")})
+ output1 <- mcmc(data=beta.samp, start=burnin+1, end=burnin + mcmc, thin=thin)
+ output2 <- mcmc(data=sigma.samp, start=burnin+1, end=burnin + mcmc, thin=thin)
+ if(m>1){
+ varnames(output1) <- sapply(c(1:ns),
+ function(i){
+ paste(xnames, "_regime", i, sep = "")
+ })
+ varnames(output2) <- sapply(c(1:ns),
+ function(i){
+ paste("sigma2_regime", i, sep = "")
+ })
+ }
+ output <- as.mcmc(cbind(output1, output2))
+
+ attr(output, "title") <- "testpanelGroupBreak Posterior Sample"
+ attr(output, "call") <- cl
+ attr(output, "y") <- resid[1:ntime]
+ attr(output, "m") <- m
+ attr(output, "nsubj") <- nsubj
+ attr(output, "ntime") <- ntime
+ if(m>0){
+ ps.holder <- matrix(posterior$psout, ntime, ns)
+ attr(output, "prob.state") <- ps.holder/nstore
+ }
+ attr(output, "logmarglike") <- posterior$logmarglikeholder
+ attr(output, "loglike") <- posterior$loglikeholder
+
+ ## report the results in a simple manner
+
+
+ return(output)
+ }
diff --git a/R/testpanelSubjectBreak.R b/R/testpanelSubjectBreak.R
new file mode 100644
index 0000000..7128cf4
--- /dev/null
+++ b/R/testpanelSubjectBreak.R
@@ -0,0 +1,103 @@
+####################################################################
+## test subject-level breaks from panel residuals
+##
+## written by Jong Hee Park 03/2009
+## modified and integrated with other codes by JHP 07/2011
+######################################################################
+
+"testpanelSubjectBreak" <-
+ function(subject.id, time.id, resid, max.break=2, minimum = 10,
+ mcmc=1000, burnin=1000, thin=1, verbose=0,
+ b0, B0, c0, d0, a = NULL, b = NULL, seed = NA,
+ Time = NULL, start.id = NULL, end.id = NULL, ps.out = FALSE){
+ ## seeds
+ seeds <- form.seeds(seed)
+ lecuyer <- seeds[[1]]
+ seed.array <- seeds[[2]]
+ lecuyer.stream <- seeds[[3]]
+
+ ## Data
+ N <- length(subject.id)
+
+ ## groupinfo matrix
+ ## col1: subj ID, col2: offset (first time C indexing), col3: #time periods
+ if (min(subject.id) != 1){
+ stop("subject.id should start 1!")
+ }
+ if (min(time.id) != 1){
+ stop("time.id should start 1!")
+ }
+ if (is.null(Time)){
+ Time <- rep(N, 1)
+ }
+ if (is.null(start.id)){
+ start.id <- min(subject.id)
+ }
+ if (is.null(end.id)){
+ end.id <- max(subject.id)
+ }
+ NC <- length(unique(subject.id))
+ time.list <- as.numeric(table(subject.id))
+
+ ## Make a residula list
+ resid.list <- as.list(rep(NA, NC))
+ start <- 1; end <- 0
+ for (i in 1:NC){
+ end <- start + time.list[i] - 1
+ resid.list[[i]] <- ts(resid[start:end], start=Time[start])
+ start <- end + 1
+ }
+ ## Do the break analysis
+ BFout <- matrix(NA, NC, max.break + 1)
+ if (ps.out ==TRUE){
+ psout <- NULL
+ }
+ else {
+ psout <- array(NA, c(max(time.list), sum(2:(max.break+1)), NC))
+ }
+ for (i in start.id:end.id){
+ residual <- resid.list[[i]]
+ nk <- length(residual)
+ out <- as.list(rep(NA, max.break))
+ if(nk > minimum){
+ for (k in 0:max.break){
+ out[[k+1]] <- MCMCpack:::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")
+ if (ps.out ==TRUE&k>0){
+ if(k==1){
+ start <- 1
+ }
+ else{
+ start <- sum(2:k)+1
+ }
+ probstate <- attr(out[[k+1]], "prob.state")
+ psout[1:length(probstate[,1]), start:(start+k), i] <- probstate
+ }
+ ## if no convergence diagnostic
+ BFout[i, k+1] <- attr(out[[k+1]], "logmarglike")
+ }
+ }
+ if (verbose > 0){
+ cat("\n ------------------------------------------------------------- ")
+ cat("\n Break analysis for subject=", i, "is just finished! \n")
+ }
+ }
+ if (ps.out ==TRUE){
+ attr(BFout, "psout") <- psout
+ }
+ model.prob.mat <- matrix(NA, NC, max.break + 1)
+ for (i in 1:NC){
+ model.prob <- exp(BFout[i, ])/sum(exp(BFout[i, ]))
+ winner <- which.max(model.prob)
+ if (verbose > 0){
+ cat("\nPr(no residual break) for subject", i, "=",
+ model.prob[1])
+ }
+ model.prob.mat[i,] <- model.prob
+ }
+ attr(BFout, "model.prob") <- model.prob.mat
+ return(BFout)
+ }
+
diff --git a/R/zzz.R b/R/zzz.R
index 8af0a9b..ea71f22 100644
--- a/R/zzz.R
+++ b/R/zzz.R
@@ -18,14 +18,14 @@
this.year <- substr(date, x[1], x[1] + attr(x, "match.length") - 1)
# echo output to screen
- cat("##\n## Markov Chain Monte Carlo Package (MCMCpack)\n")
- cat("## Copyright (C) 2003-", this.year,
- " Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park\n", sep="")
- cat("##\n## Support provided by the U.S. National Science Foundation\n")
- cat("## (Grants SES-0350646 and SES-0350613)\n##\n")
- require(coda, quietly=TRUE)
- require(MASS, quietly=TRUE)
- require(stats, quietly=TRUE)
+ packageStartupMessage("##\n## Markov Chain Monte Carlo Package (MCMCpack)")
+ packageStartupMessage("## Copyright (C) 2003-", this.year,
+ " Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park", sep="")
+ packageStartupMessage("##\n## Support provided by the U.S. National Science Foundation")
+ packageStartupMessage("## (Grants SES-0350646 and SES-0350613)\n##")
+ ## require(coda, quietly=TRUE)
+ ## require(MASS, quietly=TRUE)
+ ## require(stats, quietly=TRUE)
}
.onUnload <- function(libpath) {
diff --git a/data/PErisk.rda b/data/PErisk.rda
index 7688d70..aaf01b7 100644
Binary files a/data/PErisk.rda and b/data/PErisk.rda differ
diff --git a/data/SupremeCourt.rda b/data/SupremeCourt.rda
index 1f1d1e7..9ffadab 100644
Binary files a/data/SupremeCourt.rda and b/data/SupremeCourt.rda differ
diff --git a/man/HMMpanelFE.Rd b/man/HMMpanelFE.Rd
new file mode 100644
index 0000000..c7e4b3a
--- /dev/null
+++ b/man/HMMpanelFE.Rd
@@ -0,0 +1,205 @@
+\name{HMMpanelFE}
+
+\alias{HMMpanelFE}
+
+\title{Markov Chain Monte Carlo for the Hidden Markov Fixed-effects Model}
+
+\description{HMMpanelFE generates a sample from the posterior
+ distribution of the fixed-effects model with varying individual effects model discussed in Park (2011).
+ The code works for both balanced and unbalanced panel data as long as there is no missing data in the middle of each group.
+ This model uses a multivariate
+ Normal prior for the fixed effects parameters and varying individual effects, an Inverse-Gamma
+ prior on the residual error variance, and Beta prior for transition probabilities. 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{HMMpanelFE(subject.id, y, X, m,
+ mcmc=1000, burnin=1000, thin=1, verbose=0,
+ b0=0, B0=0.001, c0 = 0.001, d0 = 0.001, delta0=0, Delta0=0.001,
+ a = NULL, b = NULL, seed = NA, ...)}
+
+\arguments{
+
+ \item{subject.id}{A numeric vector indicating the group number. It should start from 1.}
+
+ \item{y}{The response variable.}
+
+ \item{X}{The model matrix excluding the constant.}
+
+ \item{m}{A vector of break numbers for each subject in the panel.}
+
+
+ \item{mcmc}{The number of MCMC iterations after burn-in.}
+
+ \item{burnin}{The number of burn-in iterations for the sampler.}
+
+ \item{thin}{The thinning interval used in the simulation. The number of
+ MCMC iterations must be divisible by this value.}
+
+ \item{verbose}{A switch which determines whether or not the progress of
+ the sampler is printed to the screen. If \code{verbose} is greater
+ than 0, the iteration number and the posterior density samples are printed to the screen every \code{verbose}th iteration.}
+
+ \item{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{delta0}{The prior mean of \eqn{\alpha}{alpha}.}
+
+ \item{Delta0}{The prior precision of \eqn{\alpha}{alpha}.}
+
+ \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{seed}{The seed for the random number generator. If NA, current R system seed is used.}
+
+ \item{...}{further arguments to be passed}
+
+}
+
+\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:
+ \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:
+ \deqn{\beta \sim \mathcal{N}(b_0,B_0^{-1})}{beta ~ N(b0,B0^(-1))}
+ And:
+ \deqn{\sigma^{-2} \sim \mathcal{G}amma(c_0/2, d_0/2)}{sigma^(-2) ~
+ Gamma(c0/2, d0/2)}
+ And:
+ \deqn{\alpha \sim \mathcal{N}(delta_0,Delta_0^{-1})}{alpha ~ N(delta0, Delta0^(-1))}
+ \eqn{\beta}{beta}, \eqn{\alpha}{alpha} and \eqn{\sigma^{-2}}{sigma^(-2)} are assumed
+ \emph{a priori} independent.
+
+ And:
+ \deqn{p_{mm} \sim \mathcal{B}eta(a, b),\;\; m = 1, \ldots, M}{p_mm ~ Beta(a, b), m = 1,...,M.}
+ Where \eqn{M}{M} is the number of states.
+
+ OLS estimates are used for starting values.
+
+ }
+
+\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.
+ The object contains an attribute \code{sigma} storage matrix that contains time-varying residual variance, an attribute \code{state} storage
+ matrix that contains posterior samples of hidden states, and an attribute \code{delta} storage matrix containing time-varying intercepts.
+}
+
+\references{
+ Jong Hee Park, 2011. ``A Unified Method for Dynamic and Cross-Sectional Heterogeneity:
+ Introducing Hidden Markov Panel Models." Working Paper.
+
+ Siddhartha Chib. 1998. ``Estimation and comparison of multiple change-point models.''
+ \emph{Journal of Econometrics}. 86: 221-241.
+
+}
+
+\examples{
+\dontrun{
+ ## data generating
+ set.seed(1974)
+ N <- 30
+ T <- 80
+ NT <- N*T
+
+ ## true parameter values
+ true.beta <- c(1, 1)
+ true.sigma <- 3
+ x1 <- rnorm(NT)
+ x2 <- runif(NT, 2, 4)
+
+ ## group-specific breaks
+ break.point = rep(T/2, N); break.sigma=c(rep(1, N));
+ break.list <- rep(1, N)
+
+ X <- as.matrix(cbind(x1, x2), NT, );
+ y <- rep(NA, NT)
+ id <- rep(1:N, each=NT/N)
+ K <- ncol(X);
+ true.beta <- as.matrix(true.beta, K, 1)
+
+ ## compute the break probability
+ ruler <- c(1:T)
+ W.mat <- matrix(NA, T, N)
+ for (i in 1:N){
+ W.mat[, i] <- pnorm((ruler-break.point[i])/break.sigma[i])
+ }
+ Weight <- as.vector(W.mat)
+
+ ## draw time-varying individual effects and sample y
+ j = 1
+ true.sigma.alpha <- 30
+ true.alpha <- true.alpha1 <- true.alpha2 <- rep(NA, N)
+ for (i in 1:N){
+ Xi <- X[j:(j+T-1), ]
+ true.mean <- Xi%*%true.beta
+ 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])
+ j <- j + T
+ }
+
+ ## extract the standardized residuals from the OLS with fixed-effects
+ FEols <- lm(y ~ X + as.factor(id) -1 )
+ resid.all <- rstandard(FEols)
+ time.id <- rep(1:80, N)
+
+ ## model fitting
+ G <- 100
+ BF <- testpanelSubjectBreak(subject.id=id, time.id=time.id,
+ resid= resid.all, max.break=3, minimum = 10,
+ mcmc=G, burnin = G, thin=1, verbose=G,
+ b0=0, B0=1/100, c0=2, d0=2, Time = time.id)
+
+ ## get the estimated break numbers
+ estimated.breaks <- make.breaklist(BF, thresh=3)
+
+ ## model fitting
+ out <- HMMpanelFE(subject.id = id, y, X=X, m = estimated.breaks,
+ mcmc=G, burnin=G, thin=1, verbose=G,
+ b0=0, B0=1/1000, c0=2, d0=2, delta0=0, Delta0=1/1000)
+
+ ## print out the slope estimate
+ ## true values are 1 and 1
+ summary(out)
+
+ ## compare them with the result from the constant fixed-effects
+ summary(FEols)
+}
+}
+\keyword{models}
diff --git a/man/HMMpanelRE.Rd b/man/HMMpanelRE.Rd
new file mode 100644
index 0000000..50f0a5c
--- /dev/null
+++ b/man/HMMpanelRE.Rd
@@ -0,0 +1,226 @@
+\name{HMMpanelRE}
+
+\alias{HMMpanelRE}
+
+\title{Markov Chain Monte Carlo for the Hidden Markov Random-effects Model}
+
+\description{HMMpanelRE generates a sample from the posterior
+ distribution of the hidden Markov random-effects model discussed in Park (2011).
+ The code works for panel data with the same starting point.
+ The sampling of panel parameters is based on Algorithm 2 of Chib and Carlin (1999). This model uses a multivariate
+ Normal prior for the fixed effects parameters and varying individual effects, an Inverse-Wishart prior on the
+ random-effects parameters, an Inverse-Gamma prior on the residual error variance, and Beta prior for transition probabilities.
+ 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{
+HMMpanelRE(subject.id, time.id, y, X, W, m=1,
+ mcmc=1000, burnin=1000, thin=1, verbose=0,
+ b0=0, B0=0.001, c0 = 0.001, d0 = 0.001, r0, R0, a = NULL, b = NULL,
+ seed = NA, beta.start = NA, sigma2.start = NA, D.start= NA, P.start = NA,
+ marginal.likelihood = c("none", "Chib95"), ...)}
+
+\arguments{
+
+ \item{subject.id}{A numeric vector indicating the group number. It should start from 1.}
+
+ \item{time.id}{A numeric vector indicating the time unit. It should start from 1.}
+
+ \item{y}{The dependent variable}
+
+ \item{X}{The model matrix of the fixed-effects}
+
+ \item{W}{The model matrix of the random-effects. W should be a subset of X.}
+
+ \item{m}{The number of changepoints.}
+
+ \item{mcmc}{The number of MCMC iterations after burn-in.}
+
+ \item{burnin}{The number of burn-in iterations for the sampler.}
+
+ \item{thin}{The thinning interval used in the simulation. The number of
+ MCMC iterations must be divisible by this value.}
+
+ \item{verbose}{A switch which determines whether or not the progress of
+ the sampler is printed to the screen. If \code{verbose} is greater
+ than 0, the iteration number and the posterior density samples are printed to the screen every \code{verbose}th iteration.}
+
+ \item{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{r0}{The shape parameter for the Inverse-Wishart prior on variance
+ matrix for the random effects. Set r=q for an uninformative prior where q is the number of random effects}
+
+ \item{R0}{The scale matrix for the Inverse-Wishart prior on variance matrix
+ for the random effects. This must be a square q-dimension
+ matrix. Use plausible variance regarding random effects for the
+ diagonal of R.}
+
+ \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{seed}{The seed for the random number generator. If NA, current R system seed is used.}
+
+ \item{beta.start}{The starting values for the beta vector. This can either be a scalar or a column vector with dimension equal to the number of betas. The default value of NA will use draws from the Uniform distribution with the same boundary with the data as the starting value. If this is a scalar, that value will serve as the starting value mean for all of the betas. When there is no covariate, the log value of means should be used.}
+
+ \item{sigma2.start}{The starting values for \eqn{\sigma^2}{sigma^2}. This can either be a scalar or a column vector with dimension equal to the number of states.}
+
+ \item{D.start}{The starting values for the beta vector. This can either be a scalar or a column vector with dimension equal to the number of betas. The default value of NA will use draws from the Uniform distribution with the same boundary with the data as the starting value. If this is a scalar, that value will serve as the starting value mean for all of the betas. When there is no covariate, the log value of means should be used.}
+
+ \item{P.start}{The starting values for the transition matrix. A user should provide a square matrix with dimension equal to the number of states. By default, draws from the \code{Beta(0.9, 0.1)} are used to construct a proper transition matrix for each raw except the last raw.}
+
+ \item{marginal.likelihood}{How should the marginal likelihood be
+ calculated? Options are: \code{none} in which case the marginal
+ likelihood will not be calculated and \code{Chib95} in which case the method of Chib (1995) is used.}
+
+ \item{...}{further arguments to be passed}
+
+}
+
+\details{
+
+ \code{HMMpanelRE} simulates from the random-effect hidden Markov panel model introduced by Park (2011).
+
+ The model takes the following form:
+ \deqn{y_i = X_i \beta_m + W_i b_i + \varepsilon_i\;\; m = 1, \ldots, M}{y_i = X_i * beta_m + W_i *
+ b_i + epsilon_i, m = 1,..., M.}
+ Where each group \eqn{i} have \eqn{k_i} observations.
+ Random-effects parameters are assumed to be time-varying at the system level:
+ \deqn{b_i \sim \mathcal{N}_q(0, D_m)}{b_i ~ N_q(0, D_m)}
+ \deqn{\varepsilon_i \sim \mathcal{N}(0, \sigma^2_m I_{k_i})}{epsilon_i ~ N(0, sigma^2_m I_{k_i})}
+
+ And the errors:
+ We assume standard, conjugate priors:
+ \deqn{\beta \sim \mathcal{N}_p(b0, B0)}{beta ~ N_p(b0, B0)}
+ And:
+ \deqn{\sigma^{2} \sim \mathcal{IG}amma(c0/2, d0/2)}{sigma^2 ~ IGamma(c0/2, d0/2)}
+ And:
+ \deqn{D \sim \mathcal{IW}ishart(r0, R0)}{D ~ IWishart(r0, R0)}
+ See Chib and Carlin (1999) for more details.
+
+ And:
+ \deqn{p_{mm} \sim \mathcal{B}eta(a, b),\;\; m = 1, \ldots, M}{p_mm ~ Beta(a, b), m = 1,...,M.}
+ Where \eqn{M}{M} is the number of states.
+
+ \emph{NOTE:} We do not provide default parameters for the priors on
+ the precision matrix for the random effects. When fitting one of
+ these models, it is of utmost importance to choose a prior that
+ reflects your prior beliefs about the random effects. Using the
+ \code{dwish} and \code{rwish} functions might be useful in choosing
+ these values.
+
+ }
+
+\value{
+ An mcmc object that contains the posterior sample. This
+ object can be summarized by functions provided by the coda package.
+ The object contains an attribute \code{prob.state} storage matrix that contains the probability of \eqn{state_i}{state_i}
+ for each period, and
+ the log-marginal likelihood of the model (\code{logmarglike}).
+}
+
+\references{
+ Jong Hee Park, 2011. ``A Unified Method for Dynamic and Cross-Sectional Heterogeneity:
+ Introducing Hidden Markov Panel Models." Working Paper.
+
+ 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
+ set.seed(1977)
+ Q <- 3
+ true.beta1 <- c(1, 1, 1) ; true.beta2 <- c(-1, -1, -1)
+ true.sigma2 <- c(2, 5); true.D1 <- diag(.5, Q); true.D2 <- diag(2.5, Q)
+ N=30; T=100;
+ NT <- N*T
+ x1 <- runif(NT, 1, 2)
+ x2 <- runif(NT, 1, 2)
+ X <- cbind(1, x1, x2); W <- X; y <- rep(NA, NT)
+
+ ## true break numbers are one and at the center
+ break.point = rep(T/2, N); break.sigma=c(rep(1, N));
+ break.list <- rep(1, N)
+ id <- rep(1:N, each=NT/N)
+ K <- ncol(X);
+ ruler <- c(1:T)
+
+ ## compute the weight for the break
+ W.mat <- matrix(NA, T, N)
+ for (i in 1:N){
+ W.mat[, i] <- pnorm((ruler-break.point[i])/break.sigma[i])
+ }
+ Weight <- as.vector(W.mat)
+
+ ## data generating by weighting two means and variances
+ j = 1
+ for (i in 1:N){
+ Xi <- X[j:(j+T-1), ]
+ Wi <- W[j:(j+T-1), ]
+ true.V1 <- true.sigma2[1]*diag(T) + Wi\%*\%true.D1\%*\%t(Wi)
+ true.V2 <- true.sigma2[2]*diag(T) + Wi\%*\%true.D2\%*\%t(Wi)
+ true.mean1 <- Xi\%*\%true.beta1
+ true.mean2 <- Xi\%*\%true.beta2
+ weight <- Weight[j:(j+T-1)]
+ y[j:(j+T-1)] <- (1-weight)*true.mean1 + (1-weight)*chol(true.V1)\%*\%rnorm(T) +
+ weight*true.mean2 + weight*chol(true.V2)\%*\%rnorm(T)
+ j <- j + T
+ }
+ ## model fitting
+ subject.id <- c(rep(1:N, each=T))
+ time.id <- c(rep(1:T, N))
+
+ ## model fitting
+ G <- 100
+ b0 <- rep(0, K) ; B0 <- solve(diag(100, K))
+ c0 <- 2; d0 <- 2
+ r0 <- 5; R0 <- diag(c(1, 0.1, 0.1))
+ subject.id <- c(rep(1:N, each=T))
+ time.id <- c(rep(1:T, N))
+ out1 <- HMMpanelRE(subject.id, time.id, y, X, W, m=1,
+ mcmc=G, burnin=G, thin=1, verbose=G,
+ b0=b0, B0=B0, c0=c0, d0=d0, r0=r0, R0=R0)
+
+ ## latent state changes
+ plotState(out1)
+
+ ## print mcmc output
+ summary(out1)
+
+
+
+}
+}
+\keyword{models}
diff --git a/man/MCMCbinaryChange.Rd b/man/MCMCbinaryChange.Rd
index e9fcbe0..3678c76 100755
--- a/man/MCMCbinaryChange.Rd
+++ b/man/MCMCbinaryChange.Rd
@@ -87,12 +87,8 @@
\author{Jong Hee Park, \email{jhp at uchicago.edu}, \url{http://home.uchicago.edu/~jhp/}.}
\references{
- Siddhartha Chib. 1995. ``Marginal Likelihood from the Gibbs Output.''
- \emph{Journal of the American Statistical Association}. 90:
- 1313-1321.
-
- Siddhartha Chib. 1998. ``Estimation and comparison of multiple change-point models.''
- \emph{Journal of Econometrics}. 86: 221-241.
+ 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.
Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011.
``MCMCpack: Markov Chain Monte Carlo in R.'',
@@ -103,6 +99,10 @@
\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.
+
}
\examples{
diff --git a/man/MCMCoprobitChange.Rd b/man/MCMCoprobitChange.Rd
index fe4a735..55bc1ae 100644
--- a/man/MCMCoprobitChange.Rd
+++ b/man/MCMCoprobitChange.Rd
@@ -124,16 +124,23 @@
Where \eqn{M}{M} is the number of states, and \eqn{\gamma_{c, m}}{gamma_(c, m)} and \eqn{\beta_m}{beta_m}
are paramters when a state is \eqn{m}{m} at \eqn{t}{t}.
+ We assume Gaussian distribution for prior of \eqn{\beta}{beta}:
+ \deqn{\beta_m \sim \mathcal{N}(b_0,B_0^{-1}),\;\; m = 1, \ldots, M}{beta_m ~ N(b0,B0^(-1)), m = 1,...,M.}
+
+ And:
+ \deqn{p_{mm} \sim \mathcal{B}eta(a, b),\;\; m = 1, \ldots, M}{p_mm ~ Beta(a, b), m = 1,...,M.}
+ Where \eqn{M}{M} is the number of states.
+
+
Note that when the fitted changepoint model has very few observations in any of states,
the marginal likelihood outcome can be ``nan," which indicates that too many breaks are assumed
given the model and data.
}
-\references{
-
+\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}.
+ 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.'',
diff --git a/man/MCMCpoissonChange.Rd b/man/MCMCpoissonChange.Rd
index ff94937..fea7a03 100644
--- a/man/MCMCpoissonChange.Rd
+++ b/man/MCMCpoissonChange.Rd
@@ -97,8 +97,9 @@
We assume Gaussian distribution for prior of \eqn{\beta}{beta}:
\deqn{\beta_m \sim \mathcal{N}(b_0,B_0^{-1}),\;\; m = 1, \ldots, M}{beta_m ~ N(b0,B0^(-1)), m = 1,...,M.}
+
And:
- \deqn{p_{mm} \sim \mathcal{B}eta{a}{b},\;\; m = 1, \ldots, k}{p_mm ~ Beta(a, b), m = 1,...,M.}
+ \deqn{p_{mm} \sim \mathcal{B}eta(a, b),\;\; m = 1, \ldots, M}{p_mm ~ Beta(a, b), m = 1,...,M.}
Where \eqn{M}{M} is the number of states.
}
diff --git a/man/MCMCprobitChange.Rd b/man/MCMCprobitChange.Rd
index c0963eb..090b937 100644
--- a/man/MCMCprobitChange.Rd
+++ b/man/MCMCprobitChange.Rd
@@ -107,13 +107,22 @@
Pr(y_t = 1) = Phi(x_i'beta_m)}
Where \eqn{M}{M} is the number of states, and \eqn{\beta_m}{beta_m}
is a parameter when a state is \eqn{m}{m} at \eqn{t}{t}.
+
+ We assume Gaussian distribution for prior of \eqn{\beta}{beta}:
+ \deqn{\beta_m \sim \mathcal{N}(b_0,B_0^{-1}),\;\; m = 1, \ldots, M}{beta_m ~ N(b0,B0^(-1)), m = 1,...,M.}
+
+ And:
+ \deqn{p_{mm} \sim \mathcal{B}eta(a, b),\;\; m = 1, \ldots, M}{p_mm ~ Beta(a, b), m = 1,...,M.}
+ Where \eqn{M}{M} is the number of states.
+
+
}
\author{Jong Hee Park, \email{jhp at uchicago.edu}, \url{http://home.uchicago.edu/~jhp/}.}
\references{
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}.
+ 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.'',
@@ -128,6 +137,7 @@ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011.
}
\examples{
+\dontrun{
set.seed(1973)
x1 <- rnorm(300, 0, 1)
true.beta <- c(-.5, .2, 1)
@@ -166,6 +176,7 @@ BayesFactor(out0, out1, out2, out3)
plotState(out2)
plotChangepoint(out2)
}
+}
\keyword{models}
diff --git a/man/MCMCregress.Rd b/man/MCMCregress.Rd
index 43847b8..022c69c 100644
--- a/man/MCMCregress.Rd
+++ b/man/MCMCregress.Rd
@@ -118,7 +118,7 @@ MCMCregress(formula, data = NULL, burnin = 1000, mcmc = 10000,
}
\references{
-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/MCMCregress.Rd b/man/MCMCresidualBreakAnalysis.Rd
similarity index 60%
copy from man/MCMCregress.Rd
copy to man/MCMCresidualBreakAnalysis.Rd
index 43847b8..39e849b 100644
--- a/man/MCMCregress.Rd
+++ b/man/MCMCresidualBreakAnalysis.Rd
@@ -1,27 +1,54 @@
-\name{MCMCregress}
-\alias{MCMCregress}
-\title{Markov Chain Monte Carlo for Gaussian Linear Regression}
+\name{MCMCresidualBreakAnalysis}
+\alias{MCMCresidualBreakAnalysis}
+\title{Break Analysis of Univariate Time Series using Markov Chain Monte Carlo}
\description{
- This function generates a sample from the posterior distribution
- of a linear regression model with Gaussian errors using
- Gibbs sampling (with a multivariate Gaussian prior on the
- beta vector, and an inverse Gamma prior on the conditional
- error variance). The 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.
- }
+ This function performs a break analysis for univariate time series data using
+ a linear Gaussian changepoint model. The code is written mainly for an internal
+ use in \code{testpanelSubjectBreak}.}
\usage{
-MCMCregress(formula, data = NULL, burnin = 1000, mcmc = 10000,
- thin = 1, verbose = 0, seed = NA, beta.start = NA,
- b0 = 0, B0 = 0, c0 = 0.001, d0 = 0.001,
- marginal.likelihood = c("none", "Laplace", "Chib95"), ...) }
+MCMCresidualBreakAnalysis(resid, m = 1,
+ b0 = 0, B0 = 0.001, c0 = 0.1, d0 = 0.1, 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{formula}{Model formula.}
+ \item{resid}{Univariate time series}
- \item{data}{Data frame.}
+ \item{m}{The number of 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{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.}
@@ -54,34 +81,14 @@ MCMCregress(formula, data = NULL, burnin = 1000, mcmc = 10000,
scalar, that value will serve as the starting value
mean for all of the betas.}
- \item{b0}{The prior mean of \eqn{\beta}{beta}. This can either be a
- scalar or a
- column vector with dimension equal to the number of betas. If this
- takes a scalar value, then that value will serve as the prior
- mean for all of the betas.}
-
- \item{B0}{The prior precision of \eqn{\beta}{beta}. This can either be a
- scalar or a square matrix with dimensions equal to the number of betas.
- If this
- takes a scalar value, then that value times an identity matrix serves
- as the prior precision of beta. Default value of 0 is equivalent to
- an improper uniform prior for beta.}
-
- \item{c0}{\eqn{c_0/2}{c0/2} is the shape parameter for the inverse
- Gamma prior on \eqn{\sigma^2}{sigma^2} (the variance of the
- disturbances). The amount of information in the inverse Gamma prior
- is something like that from \eqn{c_0}{c0} pseudo-observations.}
-
- \item{d0}{\eqn{d_0/2}{d0/2} is the scale parameter for the
- inverse Gamma prior on \eqn{\sigma^2}{sigma^2} (the variance of the
- disturbances). In constructing the inverse Gamma prior,
- \eqn{d_0}{d0} acts like the sum of squared errors from the
- \eqn{c_0}{c0} pseudo-observations.}
+ \item{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
+ \item{marginal.likelihood}{How should the marginal likelihood be
calculated? Options are: \code{none} in which case the marginal
- likelihood will not be calculated, \code{Laplace} in which case the
- Laplace approximation (see Kass and Raftery, 1995) is used, and
+ likelihood will not be calculated, and
\code{Chib95} in which case the method of Chib (1995) is used.}
\item{...}{further arguments to be passed}
@@ -93,7 +100,7 @@ MCMCregress(formula, data = NULL, burnin = 1000, mcmc = 10000,
}
\details{
- \code{MCMCregress} simulates from the posterior distribution using
+ \code{MCMCresidualBreakAnalysis} simulates from the posterior distribution using
standard Gibbs sampling (a multivariate Normal draw for the betas, and an
inverse Gamma draw for the conditional error variance). The simulation
proper is done in compiled C++ code to maximize efficiency. Please consult
@@ -101,50 +108,37 @@ MCMCregress(formula, data = NULL, burnin = 1000, mcmc = 10000,
used to analyze the posterior sample.
The model takes the following form:
- \deqn{y_i = x_i ' \beta + \varepsilon_{i}}{y_i = x_i'beta + epsilon_i}
- Where the errors are assumed to be Gaussian:
- \deqn{\varepsilon_{i} \sim \mathcal{N}(0, \sigma^2)}{epsilon_i ~ N(0,
- sigma^2)}
- We assume standard, semi-conjugate priors:
+ \deqn{y_{i} \sim \mathcal{N}(\beta_{m}, \sigma^2_{m}) \;\; m = 1, \ldots, M}{y_i ~ N(beta_m, sigma^2_m), m = 1,...,M.}
+
+ We assume standard, semi-conjugate priors:
\deqn{\beta \sim \mathcal{N}(b_0,B_0^{-1})}{beta ~ N(b0,B0^(-1))}
And:
\deqn{\sigma^{-2} \sim \mathcal{G}amma(c_0/2, d_0/2)}{sigma^(-2) ~
Gamma(c0/2, d0/2)}
Where \eqn{\beta}{beta} and \eqn{\sigma^{-2}}{sigma^(-2)} are assumed
- \emph{a priori} independent. Note that only starting values for
- \eqn{\beta}{beta} are allowed because simulation is done using
- Gibbs sampling with the conditional error variance
- as the first block in the sampler.
+ \emph{a priori} independent.
+
+ And:
+ \deqn{p_{mm} \sim \mathcal{B}eta(a, b),\;\; m = 1, \ldots, M}{p_mm ~ Beta(a, b), m = 1,...,M.}
+ Where \eqn{M}{M} is the number of states.
}
\references{
-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/}.
-
- Siddhartha Chib. 1995. ``Marginal Likelihood from the Gibbs Output.''
- \emph{Journal of the American Statistical Association}. 90:
- 1313-1321.
-
- Robert E. Kass and Adrian E. Raftery. 1995. ``Bayes Factors.''
- \emph{Journal of the American Statistical Association}. 90: 773-795.
-
- Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007.
- \emph{Scythe Statistical Library 1.0.} \url{http://scythe.wustl.edu}.
+ Jong Hee Park, 2011. ``A Unified Method for Dynamic and Cross-Sectional Heterogeneity:
+ Introducing Hidden Markov Panel Models." Working Paper.
- 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. 1998. ``Estimation and comparison of multiple change-point models.''
+ \emph{Journal of Econometrics}. 86: 221-241.
}
\examples{
\dontrun{
line <- list(X = c(-2,-1,0,1,2), Y = c(1,3,3,3,5))
-posterior <- MCMCregress(Y~X, data=line, verbose=1000)
-plot(posterior)
-raftery.diag(posterior)
+ols <- lm(Y~X)
+residual <- rstandard(ols)
+posterior <- MCMCresidualBreakAnalysis(residual, m = 1, data=line, mcmc=1000, verbose=200)
+plotState(posterior)
summary(posterior)
}
}
diff --git a/man/make.breaklist.Rd b/man/make.breaklist.Rd
new file mode 100644
index 0000000..14aefa6
--- /dev/null
+++ b/man/make.breaklist.Rd
@@ -0,0 +1,30 @@
+\name{make.breaklist}
+\alias{make.breaklist}
+\title{Vector of break numbers}
+\description{This function generates a vector of break numbers using the output of \code{testpanelSubjectBreak}.
+The function performs a pairwise comparison of models using Bayes Factors. }
+
+\usage{
+ make.breaklist(BF, threshold=3)
+}
+
+\arguments{
+
+\item{BF}{output of \code{testpanelSubjectBreak}.}
+
+\item{threshold}{The Bayes Factor threshold to pick the best model.
+ If a Bayes factor of two models is smaller than \code{threshold},
+ the model with a smaller number of break is chosen to avoid the over-identification problem.
+ Users can change threshold into any positive number.
+ The default value of 3 is chosen as it indicates the existence of "substantial evidence" in favor of the model in the numerator
+ according to Jeffreys' scale.}
+
+}
+\references{
+ Jong Hee Park, 2011. ``A Unified Method for Dynamic and Cross-Sectional Heterogeneity:
+ Introducing Hidden Markov Panel Models." Working Paper.
+
+ Harold Jeffreys, 1961. The Theory of Probability. Oxford University Press.
+}
+
+\seealso{\code{\link{testpanelSubjectBreak}}}
diff --git a/man/testpanelGroupBreak.Rd b/man/testpanelGroupBreak.Rd
new file mode 100644
index 0000000..02650dd
--- /dev/null
+++ b/man/testpanelGroupBreak.Rd
@@ -0,0 +1,179 @@
+\name{testpanelGroupBreak}
+
+\alias{testpanelGroupBreak}
+
+\title{A Test for the Group-level Break using a Multivariate Linear Regression Model with Breaks}
+
+\description{testpanelGroupBreak fits a multivariate linear regression model with parametric breaks using
+ panel residuals to test the existence of group-level breaks in panel residuals. The details are discussed in Park (2011).}
+
+\usage{
+testpanelGroupBreak(subject.id, time.id, resid, m=1,
+ mcmc=1000, burnin=1000, thin=1, verbose=0,
+ b0, B0, c0, d0, a = NULL, b = NULL,
+ seed = NA, marginal.likelihood = c("none", "Chib95"), ...)}
+
+\arguments{
+
+ \item{subject.id}{A numeric vector indicating the group number. It should start from 1.}
+
+ \item{time.id}{A numeric vector indicating the time unit. It should start from 1.}
+
+ \item{resid}{A vector of panel residuals}
+
+ \item{m}{The number of changepoints.}
+
+ \item{mcmc}{The number of MCMC iterations after burn-in.}
+
+ \item{burnin}{The number of burn-in iterations for the sampler.}
+
+ \item{thin}{The thinning interval used in the simulation. The number of
+ MCMC iterations must be divisible by this value.}
+
+ \item{verbose}{A switch which determines whether or not the progress of
+ the sampler is printed to the screen. If \code{verbose} is greater
+ than 0, the iteration number and the posterior density samples are printed to the screen every \code{verbose}th iteration.}
+
+ \item{b0}{The prior mean of the residual mean.}
+
+ \item{B0}{The prior precision of the residual variance}
+
+ \item{c0}{\eqn{c_0/2}{c0/2} is the shape parameter for the inverse
+ Gamma prior on \eqn{\sigma^2}{sigma^2}. 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}.}
+
+ \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{seed}{The seed for the random number generator. If NA, current R system seed is used.}
+
+ \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}
+
+}
+
+\details{
+
+ \code{testpanelGroupBreak} fits a multivariate linear regression model with parametric breaks using panel residuals
+ to detect the existence of system-level breaks in unobserved factors as discussed in Park (2011).
+
+ The model takes the following form:
+ \deqn{e_{i} \sim \mathcal{N}(\beta_{m}, \sigma^2_m I)\;\; m = 1, \ldots, M}{epsilon_i ~ N(beta_m, sigma^2_m I_{k_i}), m = 1,..., M.}
+
+ We assume standard, semi-conjugate priors:
+ \deqn{\beta \sim \mathcal{N}(b0, B0)}{beta ~ N(b0, B0)}
+ And:
+ \deqn{\sigma^{-2} \sim \mathcal{G}amma(c_0/2, d_0/2)}{sigma^(-2) ~
+ Gamma(c0/2, d0/2)}
+ Where \eqn{\beta}{beta} and \eqn{\sigma^{-2}}{sigma^(-2)} are assumed
+ \emph{a priori} independent.
+
+ And:
+ \deqn{p_{mm} \sim \mathcal{B}eta(a, b),\;\; m = 1, \ldots, M}{p_mm ~ Beta(a, b), m = 1,...,M.}
+ Where \eqn{M}{M} is the number of states.
+
+ }
+
+\value{
+ An mcmc object that contains the posterior sample. This
+ object can be summarized by functions provided by the coda package.
+ The object contains an attribute \code{prob.state} storage matrix that contains the probability of \eqn{state_i}{state_i}
+ for each period, and
+ the log-marginal likelihood of the model (\code{logmarglike}).
+}
+
+\references{
+ Jong Hee Park, 2011. ``A Unified Method for Dynamic and Cross-Sectional Heterogeneity:
+ Introducing Hidden Markov Panel Models." Working Paper.
+
+ 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
+ set.seed(1977)
+ Q <- 3
+ true.beta1 <- c(1, 1, 1) ; true.beta2 <- c(1, -1, -1)
+ true.sigma2 <- c(1, 3); true.D1 <- diag(.5, Q); true.D2 <- diag(2.5, Q)
+ N=20; T=100;
+ NT <- N*T
+ x1 <- rnorm(NT)
+ x2 <- runif(NT, 5, 10)
+ X <- cbind(1, x1, x2); W <- X; y <- rep(NA, NT)
+
+ ## true break numbers are one and at the center
+ break.point = rep(T/2, N); break.sigma=c(rep(1, N));
+ break.list <- rep(1, N)
+ id <- rep(1:N, each=NT/N)
+ K <- ncol(X);
+ ruler <- c(1:T)
+
+ ## compute the weight for the break
+ W.mat <- matrix(NA, T, N)
+ for (i in 1:N){
+ W.mat[, i] <- pnorm((ruler-break.point[i])/break.sigma[i])
+ }
+ Weight <- as.vector(W.mat)
+
+ ## data generating by weighting two means and variances
+ j = 1
+ for (i in 1:N){
+ Xi <- X[j:(j+T-1), ]
+ Wi <- W[j:(j+T-1), ]
+ true.V1 <- true.sigma2[1]*diag(T) + Wi\%*\%true.D1\%*\%t(Wi)
+ true.V2 <- true.sigma2[2]*diag(T) + Wi\%*\%true.D2\%*\%t(Wi)
+ true.mean1 <- Xi\%*\%true.beta1
+ true.mean2 <- Xi\%*\%true.beta2
+ weight <- Weight[j:(j+T-1)]
+ y[j:(j+T-1)] <- (1-weight)*true.mean1 + (1-weight)*chol(true.V1)\%*\%rnorm(T) +
+ weight*true.mean2 + weight*chol(true.V2)\%*\%rnorm(T)
+ j <- j + T
+ }
+ ## model fitting
+ subject.id <- c(rep(1:N, each=T))
+ time.id <- c(rep(1:T, N))
+
+
+ resid <- rstandard(lm(y ~X-1 + as.factor(subject.id)))
+ G <- 100
+ out0 <- testpanelGroupBreak(subject.id, time.id, resid, m=0,
+ mcmc=G, burnin=G, thin=1, verbose=G,
+ b0=0, B0=1/100, c0=2, d0=2, marginal.likelihood = "Chib95")
+ out1 <- testpanelGroupBreak(subject.id, time.id, resid, m=1,
+ mcmc=G, burnin=G, thin=1, verbose=G,
+ b0=0, B0=1/100, c0=2, d0=2, marginal.likelihood = "Chib95")
+ out2 <- testpanelGroupBreak(subject.id, time.id, resid, m=2,
+ mcmc=G, burnin=G, thin=1, verbose=G,
+ b0=0, B0=1/100, c0=2, d0=2, marginal.likelihood = "Chib95")
+ out3 <- testpanelGroupBreak(subject.id, time.id, resid, m=3,
+ mcmc=G, burnin=G, thin=1, verbose=G,
+ b0=0, B0=1/100, c0=2, d0=2, marginal.likelihood = "Chib95")
+
+ ## Note that the code is for a hypothesis test of no break in panel residuals.
+ ## When breaks exist, the estimated number of break in the mean and variance of panel residuals
+ ## tends to be larger than the number of break in the data generating process.
+ ## This is due to the difference in parameter space, not an error of the code.
+ BayesFactor(out0, out1, out2, out3)
+
+ ## In order to identify the number of breaks in panel parameters,
+ ## use HMMpanelRE() instead.
+
+}
+}
+\keyword{models}
diff --git a/man/testpanelSubjectBreak.Rd b/man/testpanelSubjectBreak.Rd
new file mode 100644
index 0000000..1365fd0
--- /dev/null
+++ b/man/testpanelSubjectBreak.Rd
@@ -0,0 +1,179 @@
+\name{testpanelSubjectBreak}
+
+\alias{testpanelSubjectBreak}
+
+\title{A Test for the Subject-level Break using a Unitivariate Linear Regression Model with Breaks}
+
+\description{testpanelSubjectBreak fits a unitivariate linear regression model with parametric breaks using
+ panel residuals to test the existence of subject-level breaks in panel residuals. The details are discussed in Park (2011).}
+
+\usage{testpanelSubjectBreak(subject.id, time.id, resid, max.break=2,
+ minimum = 10, mcmc=1000, burnin=1000, thin=1, verbose=0,
+ b0, B0, c0, d0, a = NULL, b = NULL, seed = NA,
+ Time = NULL, start.id = NULL, end.id = NULL, ps.out = FALSE)}
+
+\arguments{
+ \item{subject.id}{A numeric vector indicating the group number. It should start from 1.}
+
+ \item{time.id}{A numeric vector indicating the time unit. It should start from 1.}
+
+ \item{resid}{A vector of panel residuals.}
+
+ \item{max.break}{An upper bound of break numbers for the test.}
+
+ \item{minimum}{A minimum length of time series for the test. The test will skip a subject with a time series shorter than this.}
+
+ \item{mcmc}{The number of MCMC iterations after burn-in.}
+
+ \item{burnin}{The number of burn-in iterations for the sampler.}
+
+ \item{thin}{The thinning interval used in the simulation. The number of
+ MCMC iterations must be divisible by this value.}
+
+ \item{verbose}{A switch which determines whether or not the progress of
+ the sampler is printed to the screen. If \code{verbose} is greater
+ than 0, the iteration number and the posterior density samples are printed to the screen every \code{verbose}th iteration.}
+
+ \item{b0}{The prior mean of the residual mean.}
+
+ \item{B0}{The prior precision of the residual variance}
+
+ \item{c0}{\eqn{c_0/2}{c0/2} is the shape parameter for the inverse
+ Gamma prior on \eqn{\sigma^2}{sigma^2}. 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}.}
+
+ \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{seed}{The seed for the random number generator. If NA, current R system seed is used.}
+
+ \item{Time}{Times of the observations. This will be used to find the time of the first observations in panel residuals. }
+
+ \item{start.id}{If users want to check the residual breaks within a subset of panel subjects, users can specify the starting subject id and
+ the ending subject id.}
+
+ \item{end.id}{See start.id.}
+
+ \item{ps.out}{If ps.out == TRUE, state probabilities are exported. If the number of panel subjects is huge, users can turn it off to save memory.}
+
+ \item{...}{further arguments to be passed}
+
+}
+
+\details{
+
+ \code{testpanelSubjectBreak} fits a univariate linear regression model for subject-level residuals from a panel model.
+ The details are discussed in Park (2011).
+
+ The model takes the following form:
+ \deqn{e_{it} = \alpha_{im} + \varepsilon_{it}\;\; m = 1, \ldots, M}{y_it = alpha_im + epsilon_it, m = 1,...,M.}
+ The errors are assumed to be time-varying at the subject level:
+ \deqn{\varepsilon_{it} \sim \mathcal{N}(0, \sigma^2_{im})}{epsilon_it ~ N(0, sigma^2_im)}
+
+ We assume standard, semi-conjugate priors:
+ \deqn{\beta \sim \mathcal{N}(b_0,B_0^{-1})}{beta ~ N(b0,B0^(-1))}
+ And:
+ \deqn{\sigma^{-2} \sim \mathcal{G}amma(c_0/2, d_0/2)}{sigma^(-2) ~
+ Gamma(c0/2, d0/2)}
+ Where \eqn{\beta}{beta} and \eqn{\sigma^{-2}}{sigma^(-2)} are assumed
+ \emph{a priori} independent.
+
+ And:
+ \deqn{p_{mm} \sim \mathcal{B}eta(a, b),\;\; m = 1, \ldots, M}{p_mm ~ Beta(a, b), m = 1,...,M.}
+ Where \eqn{M}{M} is the number of states.
+
+ OLS estimates are used for starting values.
+
+ }
+
+\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.
+ If psout == TRUE, the returned object has an array attribute \code{psout} containing state probabilities for all HMMs.
+}
+
+\references{
+ Jong Hee Park, 2011. ``A Unified Method for Dynamic and Cross-Sectional Heterogeneity:
+ Introducing Hidden Markov Panel Models." Working Paper.
+
+ Siddhartha Chib. 1998. ``Estimation and comparison of multiple change-point models.''
+ \emph{Journal of Econometrics}. 86: 221-241.
+
+}
+
+\examples{
+\dontrun{
+ set.seed(1974)
+ N <- 30
+ T <- 80
+ NT <- N*T
+
+ ## true parameter values
+ true.beta <- c(1, 1)
+ true.sigma <- 3
+ x1 <- rnorm(NT)
+ x2 <- runif(NT, 2, 4)
+
+ ## group-specific breaks
+ break.point = rep(T/2, N); break.sigma=c(rep(1, N));
+ break.list <- rep(1, N)
+
+ X <- as.matrix(cbind(x1, x2), NT, );
+ y <- rep(NA, NT)
+ id <- rep(1:N, each=NT/N)
+ K <- ncol(X);
+ true.beta <- as.matrix(true.beta, K, 1)
+
+ ## compute the break probability
+ ruler <- c(1:T)
+ W.mat <- matrix(NA, T, N)
+ for (i in 1:N){
+ W.mat[, i] <- pnorm((ruler-break.point[i])/break.sigma[i])
+ }
+ Weight <- as.vector(W.mat)
+
+ ## draw time-varying individual effects and sample y
+ j = 1
+ true.sigma.alpha <- 30
+ true.alpha <- true.alpha1 <- true.alpha2 <- rep(NA, N)
+ for (i in 1:N){
+ Xi <- X[j:(j+T-1), ]
+ true.mean <- Xi%*%true.beta
+ 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])
+ j <- j + T
+ }
+
+ ## extract the standardized residuals from the OLS with fixed-effects
+ FEols <- lm(y ~ X + as.factor(id) -1 )
+ resid.all <- rstandard(FEols)
+ time.id <- rep(1:80, N)
+
+ ## model fitting
+ G <- 1000
+ BF <- testpanelSubjectBreak(subject.id=id, time.id=time.id,
+ resid= resid.all, max.break=3, minimum = 10,
+ mcmc=G, burnin = G, thin=1, verbose=G,
+ b0=0, B0=1/100, c0=2, d0=2, Time = time.id)
+
+ ## estimated break numbers
+ estimated.breaks <- make.breaklist(BF, thresh=3)
+
+ ## print all posterior model probabilities
+ print(attr(BF, "model.prob"))
+}
+}
+\keyword{models}
diff --git a/src/HMMmultivariateGaussian.cc b/src/HMMmultivariateGaussian.cc
new file mode 100644
index 0000000..dc0e49f
--- /dev/null
+++ b/src/HMMmultivariateGaussian.cc
@@ -0,0 +1,996 @@
+//////////////////////////////////////////////////////////////////////////
+//HMMmultivariateGaussian.cc is C++ code to estimate a Gaussian panel model with a structural break
+// y_{it} = \x'_{it}\b + \varepsilon_{it}
+// \varepsilon_{it} \sim \normdist{0}{\sigma^2}
+// Parameters = {beta, sigma2, P_i}// static + changing
+//
+// Written by Jong Hee Park
+// 11/19/2008
+// Modified by JHP
+// 07/04/2011
+//
+// Copyright (C) 2008-present Jong Hee Park
+//////////////////////////////////////////////////////////////////////////
+#ifndef HMMGAUSSIANPANELRE_CC
+#define HMMGAUSSIANPANELRE_CC
+
+#include<vector>
+#include<algorithm>
+
+#include "MCMCrng.h"
+#include "MCMCfcds.h"
+#include "matrix.h"
+#include "distributions.h"
+#include "stat.h"
+#include "la.h"
+#include "ide.h"
+#include "smath.h"
+#include "rng.h"
+#include "mersenne.h"
+#include "lecuyer.h"
+#include "lapack.h"
+
+#include <R.h> // needed to use Rprintf()
+#include <R_ext/Utils.h> // needed to allow user interrupts
+
+using namespace std;
+using namespace scythe;
+
+static double ln_invgamma(double theta, double a, double b) {
+ double logf = a * log(b) - lngammafn(a) + -(a+1) * log(theta) +
+ -b/theta;
+ return (logf);
+ //pow(b, a) / gammafn(a) * pow(theta, -(a+1)) * exp(-b/theta);
+}
+
+
+// used to access 1d arrays holding R matrices like a 2d array
+#define M(ROW,COL,NROWS) (COL*NROWS+ROW)
+
+template <typename RNGTYPE>
+void MultivariateGaussian_impl(rng<RNGTYPE>& stream,
+ unsigned int nsubj, unsigned int ntime, unsigned int nobs,
+ const Matrix<int>& subjectid, const Matrix<int>& timeid,
+ const Matrix<>& Y, const Matrix<>& X,
+ const Matrix<>& YT, const Matrix<>& XT,
+ unsigned int burnin, unsigned int mcmc,
+ unsigned int thin, unsigned int verbose,
+ double sigma2start,
+ const Matrix<>& b0, const Matrix<>& B0,
+ const double c0, const double d0,
+ const Matrix<>& time_groupinfo,
+ const Matrix<>& subject_groupinfo,
+ Matrix<>& beta_store, Matrix<>& sigma_store,
+ double& logmarglike, double& loglike,
+ unsigned int chib)
+{ // redefine constants
+ const unsigned int K = X.cols();; // ncol(X)
+ const int NOBS = nobs;
+ const Matrix<> B0inv = invpd(B0);
+
+ const int tot_iter = burnin + mcmc;
+ const int nstore = mcmc / thin; // number of draws to store
+
+ // generate posk_arr and post_arr
+ // Number of observations by group k
+ int *nobsk = new int[nsubj];
+ for (int k=0; k<nsubj; k++) {
+ nobsk[k]=0;
+ for (int n=0; n<NOBS; n++) {
+ if (subjectid[n]==k+1) {
+ nobsk[k]+=1;
+ }
+ }
+ }
+
+ // Position of each group in the data-set
+ int **posk_arr = new int*[nsubj];
+ for (int k=0; k<nsubj; k++) {
+ posk_arr[k] = new int[nobsk[k]];
+ int repk=0;
+ for (int n=0; n<NOBS; n++) {
+ if (subjectid[n]==k+1) {
+ posk_arr[k][repk]=n;
+ repk++;
+ }
+ }
+ }
+
+ // Small fixed matrices indexed on k for data access
+ Matrix<double> *Yk_arr = new Matrix<double>[nsubj];
+ Matrix<double> *Xk_arr = new Matrix<double>[nsubj];
+ for(int k=0; k<nsubj; k++) {
+ Xk_arr[k] = Matrix<double>(nobsk[k], K);
+ Yk_arr[k] = Matrix<double>(nobsk[k], 1);
+ for (int l=0; l<nobsk[k]; l++) {
+ for (int p=0; p< K; p++) {
+ Xk_arr[k](l, p) = X[p*NOBS + posk_arr[k][l]];
+ }
+ Yk_arr[k](l,0) = Y[posk_arr[k][l]];
+ }
+ }
+
+ Matrix<double> *cpXk_arr = new Matrix<double>[nsubj];
+ Matrix<double> *tXYk_arr = new Matrix<double>[nsubj];
+ for(int k=0; k<nsubj; k++) {
+ cpXk_arr[k] = crossprod(Xk_arr[k]);
+ tXYk_arr[k] = t(Xk_arr[k])*Yk_arr[k];
+ }
+
+ // starting values
+ Matrix<> beta(K, 1);
+ double sigma2 = sigma2start;
+
+ // MCMC iterations start here
+ int sampcount = 0;
+
+ Rprintf("\n ///////////////////////////////////////////////// \n");
+ Rprintf("\n MCMC for Multivariate Gaussian loop starts! \n");
+ Rprintf("\n ///////////////////////////////////////////////// \n");
+
+ /////////////////////////////////////////////////
+ // initialize Yhat for the first loop only
+ /////////////////////////////////////////////////
+ for (int iter=0; iter < tot_iter; ++iter){
+ //////////////////////
+ // Step 1. Sample beta
+ //////////////////////
+ Matrix<> XVX(K, K);
+ Matrix<> XVY(K, 1);
+ for(int s = 0; s<nsubj; ++s) {
+ XVX = XVX + cpXk_arr[s];
+ XVY = XVY + tXYk_arr[s];
+ }
+ Matrix<> beta_var = invpd(B0 + XVX/sigma2);
+ Matrix<> beta_mean = beta_var*(B0*b0 + XVY/sigma2);
+ beta = stream.rmvnorm(beta_mean, beta_var);
+
+
+ /////////////////////////////////////////////////
+ // Step 2. Sample sigma2
+ /////////////////////////////////////////////////
+ double SSE = 0;
+ int counter = 0;
+ for(int s=0;s<nsubj; ++s){
+ int ntime_s = subject_groupinfo(s, 2);
+ Matrix<> e = t(Yk_arr[s]-Xk_arr[s]*beta)*(Yk_arr[s] - Xk_arr[s]*beta);
+ SSE = SSE + e[0];
+ counter = counter + ntime_s;
+ }
+ double nu = (c0 + NOBS)/2;
+ double delta = (d0 + SSE)/2;
+ sigma2 = stream.rigamma(nu, delta);
+
+ /////////////////////////////////////////////////
+ // STORE
+ /////////////////////////////////////////////////
+ if (iter >= burnin && ((iter % thin) == 0)) {
+ for (int i=0; i< K; ++i){
+ beta_store(sampcount, i) = beta(i);// stored by the order of (11, 12, 13, 21, 22, 23)
+ }
+ sigma_store(sampcount) = sigma2;
+ ++sampcount;
+ }
+
+
+ /////////////////////////////////////////////////
+ // REPORT
+ /////////////////////////////////////////////////
+ if(verbose > 0 && iter % verbose == 0){
+ Rprintf("\n ----------------------------------------------------------------------- \n");
+ Rprintf("Multivaraite Gaussian iteration %i of %i \n", (iter+1), tot_iter);
+ Rprintf("\n beta is ");
+ for(int i=0;i<K; ++i) {
+ Rprintf("%10.5f", beta(i));
+ }
+ Rprintf("\n sigma2 = %10.5f\n", sigma2);
+ }
+ }// END of MCMC loop
+
+ //////////////////////////////////////
+ if(chib == 1){
+ //////////////////////////////////////
+ Matrix<> beta_st(K, 1);
+ beta_st(_, 0) = meanc(beta_store);
+ const double sigma2_st = mean(sigma_store);
+
+ Matrix<> density_beta(nstore, 1);
+ //////////////////////////////////////////////////////////////////
+ // 1. pdf.beta | sigma2_g
+ //////////////////////////////////////////////////////////////////
+ for (int iter = 0; iter<nstore; ++iter){
+ Matrix<> XVX(K, K);
+ Matrix<> XVY(K, 1);
+ for(int s = 0; s<nsubj; ++s) {
+ XVX = XVX + cpXk_arr[s];
+ XVY = XVY + tXYk_arr[s];
+ }
+ Matrix<> beta_var = invpd(B0 + XVX/sigma_store(iter));
+ Matrix<> beta_mean = beta_var*(B0*b0 + XVY/sigma_store(iter));
+ beta = stream.rmvnorm(beta_mean, beta_var);
+
+ if (K == 1){
+ density_beta(iter) = dnorm(beta_st(0), beta_mean[0], sqrt(beta_var[0]));
+ }
+ else{
+ density_beta(iter) = ::exp(lndmvn(beta_st, beta_mean, beta_var));
+ }
+ }
+ double pdf_beta = log(mean(density_beta));
+
+ //////////////////////////////////////////////////////////////////
+ // 2. pdf.sigma2
+ //////////////////////////////////////////////////////////////////
+ Matrix<> density_sigma2(nstore, 1);
+ for (int iter = 0; iter<nstore; ++iter){
+ double SSE = 0;
+ for(int s=0;s<nsubj; ++s){
+ int ntime_s = subject_groupinfo(s, 2);
+ Matrix<> e = t(Yk_arr[s]-Xk_arr[s]*beta_st)*(Yk_arr[s] - Xk_arr[s]*beta_st);
+ SSE = SSE + e[0];
+ }
+ double nu = (c0 + NOBS)/2;
+ double delta = (d0 + SSE)/2;
+ sigma2 = stream.rigamma(nu, delta);
+ density_sigma2(iter) = ::exp(ln_invgamma(sigma2_st, nu, delta));
+ }
+ double pdf_sigma2 = log(mean(density_sigma2));
+
+ //////////////////////////////////////////////////////////////////
+ // likelihood f(y|beta_st, D_st, sigma2_st, P_st)
+ //////////////////////////////////////////////////////////////////
+ loglike = 0;
+ for(int s = 0; s<nsubj; ++s) {
+ int ntime_s = subject_groupinfo(s, 2);
+ Matrix<> Sigma = sigma2_st*eye(ntime_s);
+ Matrix<> Mu = Xk_arr[s]*beta_st;
+ loglike += lndmvn(Yk_arr[s], Mu, Sigma);
+ }
+
+ //////////////////////
+ // log prior ordinate
+ //////////////////////
+ double density_beta_prior = 0;
+ 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);
+ }
+ double density_sigma2_prior = ln_invgamma(sigma2_st, c0/2, d0/2);
+
+ // compute marginal likelihood
+ double logprior = density_beta_prior + density_sigma2_prior;
+ logmarglike = (loglike + logprior) - (pdf_beta + pdf_sigma2);
+ if(verbose > 0){
+ Rprintf("\nlogmarglike = %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_sigma2);
+ // Rprintf("pdf_P = %10.5f\n", pdf_P);
+ }
+ }
+
+ delete[] nobsk;
+ for(int k=0; k<nsubj; k++) {
+ delete[] posk_arr[k];
+ }
+ delete [] posk_arr;
+ delete [] Xk_arr;
+ delete [] Yk_arr;
+ delete [] tXYk_arr;
+ delete [] cpXk_arr;
+
+}// end of MultivariateGaussian_impl
+
+
+template <typename RNGTYPE>
+void HMMmultivariateGaussian_impl(rng<RNGTYPE>& stream,
+ unsigned int nsubj, unsigned int ntime,
+ unsigned int m, unsigned int nobs,
+ const Matrix<int>& subjectid, const Matrix<int>& timeid,
+ const Matrix<>& Y, const Matrix<>& X,
+ const Matrix<>& YT, const Matrix<>& XT,
+ unsigned int burnin, unsigned int mcmc,
+ unsigned int thin, unsigned int verbose,
+ Matrix<>& betastart, double sigma2start,
+ const Matrix<>& b0, const Matrix<>& B0,
+ const double c0, const double d0,
+ const Matrix<>& P0,
+ const Matrix<>& time_groupinfo,
+ const Matrix<>& subject_groupinfo,
+ Matrix<>& beta_store, Matrix<>& sigma_store,
+ Matrix<>& P_store, Matrix<>& ps_store, Matrix<>& s_store,
+ double& logmarglike, double& loglike,
+ unsigned int chib)
+{ // redefine constants
+ const unsigned int K = X.cols();; // ncol(X)
+ const int NOBS = nobs;
+ const Matrix<> B0inv = invpd(B0);
+
+ const int tot_iter = burnin + mcmc;
+ const int nstore = mcmc / thin; // number of draws to store
+ const int ns = m + 1;
+
+ // generate posk_arr and post_arr
+ // Number of observations by group k
+ int *nobsk = new int[nsubj];
+ for (int k=0; k<nsubj; k++) {
+ nobsk[k]=0;
+ for (int n=0; n<NOBS; n++) {
+ if (subjectid[n]==k+1) {
+ nobsk[k]+=1;
+ }
+ }
+ }
+
+ // Position of each group in the data-set
+ int **posk_arr = new int*[nsubj];
+ for (int k=0; k<nsubj; k++) {
+ posk_arr[k] = new int[nobsk[k]];
+ int repk=0;
+ for (int n=0; n<NOBS; n++) {
+ if (subjectid[n]==k+1) {
+ posk_arr[k][repk]=n;
+ repk++;
+ }
+ }
+ }
+
+ // Small fixed matrices indexed on k for data access
+ Matrix<double> *Yk_arr = new Matrix<double>[nsubj];
+ Matrix<double> *Xk_arr = new Matrix<double>[nsubj];
+ for(int k=0; k<nsubj; k++) {
+ Xk_arr[k] = Matrix<double>(nobsk[k], K);
+ Yk_arr[k] = Matrix<double>(nobsk[k], 1);
+ for (int l=0; l<nobsk[k]; l++) {
+ for (int p=0; p< K; p++) {
+ Xk_arr[k](l, p) = X[p*NOBS + posk_arr[k][l]];
+ }
+ Yk_arr[k](l,0) = Y[posk_arr[k][l]];
+ }
+ }
+ // number of observations by time t
+ int *nobst = new int[ntime];
+ for (int k=0; k<ntime; k++) {
+ nobst[k]=0;
+ for (int n=0; n<NOBS; n++) {
+ if (timeid[n]==k + 1) {
+ nobst[k]+=1;
+ }
+ }
+ }
+ // Position of each group in the data-set
+ int **post_arr = new int*[ntime];
+ for (int k=0; k<ntime; k++) {
+ post_arr[k] = new int[nobst[k]];
+ int repk=0;
+ for (int n=0; n<NOBS; n++) {
+ if (timeid[n]==k+1) {
+ post_arr[k][repk]=n;
+ repk++;
+ }
+ }
+ }
+
+ // XTarr is data transformed for multivariate TS analsysi
+ Matrix<>* Xt_arr = new Matrix<>[ntime];
+ Matrix<>* Yt_arr = new Matrix<>[ntime];
+ for(int k=0; k<ntime; k++) {
+ if (nobst[k] > 0){
+ Xt_arr[k] = Matrix<double>(nobst[k], K);
+ Yt_arr[k] = Matrix<double>(nobst[k], 1);
+ for (int l = 0; l<nobst[k]; l++) {
+ for (int p = 0; p < K; p++) {
+ Xt_arr[k](l, p) = X[p*NOBS+post_arr[k][l]];
+ }
+ Yt_arr[k](l,0) = Y[post_arr[k][l]];
+ }
+ }
+ }
+
+ Matrix<double> *cpXt_arr = new Matrix<double>[ntime];
+ Matrix<double> *tXt_arr = new Matrix<double>[ntime];
+ Matrix<double> *tXYt_arr = new Matrix<double>[ntime];
+ for(int k=0; k<ntime; k++) {
+ cpXt_arr[k] = crossprod(Xt_arr[k]);
+ tXt_arr[k] = t(Xt_arr[k]);
+ tXYt_arr[k] = t(Xt_arr[k])*Yt_arr[k];
+ }
+
+ // starting values
+ Matrix<> beta(ns, K);
+ Matrix<> sigma2(ns, 1);
+ for (int j=0; j<ns; ++j){
+ beta(j,_) = betastart;
+ sigma2(j) = sigma2start;
+ }
+ Matrix<> P = P0;
+
+ // MCMC iterations start here
+ int sampcount = 0;
+
+ /////////////////////////////////////////////////
+ // initialize Yhat for the first loop only
+ /////////////////////////////////////////////////
+ for (int iter=0; iter < tot_iter; ++iter){
+ //////////////////////
+ // Step 1. Sample state
+ //////////////////////
+ Matrix<> F(ntime, ns);
+ Matrix<> pr1(ns, 1);
+ pr1[0] = 1;
+ Matrix<> py(ns, 1);
+ Matrix<> pstyt1(ns, 1);
+ for (int tt=0; tt<ntime ; ++tt) {
+ if(nobst[tt]>0){
+ int nsubj_s = time_groupinfo(tt, 2);
+ for (int j=0; j<ns;++j){
+ Matrix<> Sigma = eye(nsubj_s)*sigma2[j];
+ Matrix<> Mu = Xt_arr[tt]*::t(beta(j,_));
+ py[j] = ::exp(lndmvn(Yt_arr[tt], Mu, Sigma));
+ }
+ if (tt==0) pstyt1 = pr1;
+ else {
+ pstyt1 = ::t(F(tt-1,_)*P);
+ }
+ Matrix<> unnorm_pstyt = pstyt1%py;
+ Matrix<> pstyt = unnorm_pstyt/sum(unnorm_pstyt);
+ for (int j=0; j<ns ; ++j) {
+ F(tt,j) = pstyt(j);
+ }
+ }
+ }// end of F matrix filtering
+ Matrix<int> state(ntime, 1);
+ Matrix<> ps = Matrix<>(ntime, ns);
+ ps(ntime-1,_) = F(ntime-1,_);
+ state(ntime-1) = ns;
+
+ Matrix<> pstyn = Matrix<>(ns, 1);
+ double pone = 0.0;
+ int tt = ntime-2;
+
+ while (tt >= 0){
+ if(nobst[tt]>0){
+ int st = state(tt+1);
+ Matrix<> Pst_1 = ::t(P(_, st-1));
+ Matrix<> unnorm_pstyn = F(tt,_)%Pst_1;
+ pstyn = unnorm_pstyn/sum(unnorm_pstyn);
+
+ if (st==1){
+ state(tt) = 1;
+ }
+ else{
+ pone = pstyn(st-2);
+ if(stream.runif() < pone) state(tt) = st-1;
+ else state(tt) = st;
+ }
+ ps(tt,_) = pstyn;
+ --tt;
+ }
+ }// end of while loop
+
+ //////////////////////
+ // Step 2. Sample beta
+ //////////////////////
+ int beta_count = 0;
+ Matrix<int> nstate(ns, 1);
+
+ for (int j = 0; j<ns; ++j){
+ for (int i = 0; i<ntime; ++i){
+ if (state[i] == j + 1) {
+ nstate[j] = nstate[j] + 1;
+ }// end of if
+ }// end of int i<n
+ beta_count = beta_count + nstate[j];
+
+ // BETA UPDATE
+ Matrix<> XVX(K, K);
+ Matrix<> XVY(K, 1);
+ for(int tt = (beta_count - nstate[j]); tt<beta_count; ++tt) {
+ if (nobst[tt] > 0){
+ XVX = XVX + cpXt_arr[tt];
+ XVY = XVY + tXYt_arr[tt];
+ }
+ }
+ Matrix<> beta_var = invpd(B0 + XVX/sigma2[j]);
+ Matrix<> beta_mean = beta_var*(B0*b0 + XVY/sigma2[j]);
+ beta(j,_) = stream.rmvnorm(beta_mean, beta_var);
+ }
+
+ /////////////////////////////////////////////////
+ // Step 3. Sample sigma2
+ /////////////////////////////////////////////////
+ beta_count = 0;
+ Matrix<int> YN(ns, 1);
+ Matrix<> SSE(ns, 1);
+ for (int j = 0; j<ns; ++j){
+ beta_count = beta_count + nstate[j];
+ for(int s = 0; s<nsubj; ++ s) {
+ Matrix<> yj = Yk_arr[s]((beta_count - nstate[j]), 0, (beta_count - 1), 0);
+ Matrix<> Xj = Xk_arr[s]((beta_count - nstate[j]), 0, (beta_count - 1), K-1);
+ Matrix<> e = ::t(yj - Xj*::t(beta(j,_)))*(yj - Xj*::t(beta(j,_)));
+ YN[j] = YN[j] + yj.rows();
+ SSE[j] = SSE[j] + e[0];
+ }
+ double nu = c0 + (double)YN[j]*0.5;
+ double scale = d0 + SSE[j]*0.5;
+ sigma2[j] = stream.rigamma(nu, scale);
+ }
+
+ //////////////////////
+ // Step 4. Sample P
+ //////////////////////
+ double shape1 = 0;
+ double shape2 = 0;
+ for (int j =0; j<m; ++j){
+ shape1 = std::abs(P0(j,j) + nstate[j] - 1);
+ shape2 = P0(j,j+1) + 1; //
+ P(j,j) = stream.rbeta(shape1, shape2);
+ P(j,j+1) = 1 - P(j,j);
+ }
+ P(m, m) = 1; //no jump at the last state
+
+ /////////////////////////////////////////////////
+ // STORE
+ /////////////////////////////////////////////////
+ if (iter >= burnin && ((iter % thin) == 0)) {
+ Matrix<> tbeta = ::t(beta);
+ for (int i=0; i<(ns*K); ++i){
+ beta_store(sampcount, i) = tbeta[i];
+ }
+ for (int i=0; i<ns; ++i){
+ sigma_store(sampcount, i) = sigma2[i];
+ }
+ for (int j=0; j<ns*ns; ++j){
+ P_store(sampcount,j)= P[j];
+ }
+ s_store(sampcount,_) = state;
+ for (int l=0; l<ntime ; ++l){
+ ps_store(l,_) = ps_store(l,_) + ps(l,_);
+ }
+ ++sampcount;
+ }
+
+
+ /////////////////////////////////////////////////
+ // REPORT
+ /////////////////////////////////////////////////
+ if(verbose > 0 && iter % verbose == 0){
+ Rprintf("\n ----------------------------------------------------------------------- \n");
+ Rprintf("HMM Multivaraite Gaussian iteration %i of %i \n", (iter+1), tot_iter);
+ for(int j=0;j<ns; ++j) {
+ Rprintf("\n number of obs from state %i is %i", j, nstate[j]);
+ }
+ for(int j=0;j<ns; ++j) {
+ Rprintf("\n beta for state %i = ", j);
+ for(int i=0;i<K; ++i) {
+ Rprintf("%10.5f", beta(j, i));
+ }
+ }
+ for(int j=0;j<ns; ++j) {
+ Rprintf("\n sigma2 for state %i is %10.5f", j, sigma2[j]);
+ }
+ }
+
+ }// END of MCMC loop
+
+ //////////////////////////////////////
+ if(chib == 1){
+ //////////////////////////////////////
+
+ // posterior ordinate
+ 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> sigma2st = meanc(sigma_store);
+ Matrix<double, Row> sigma2_st(ns, 1);
+ for (int j = 0; j< ns; ++j){
+ sigma2_st[j] = sigma2st[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];
+ }
+
+ Matrix<> density_beta(nstore, ns);
+ //////////////////////////////////////////////////////////////////
+ // 1. pdf.beta | sigma2_g, P_g
+ //////////////////////////////////////////////////////////////////
+ for (int iter = 0; iter<nstore; ++iter){
+ int beta_count = 0;
+ Matrix<int> nstate(ns, 1);
+ for (int j = 0; j<ns; ++j){
+ for (int i = 0; i<ntime; ++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];
+
+ Matrix<> XVX(K, K);
+ Matrix<> XVY(K, 1);
+ for(int tt = (beta_count - nstate[j]); tt<beta_count; ++tt) {
+ XVX = XVX + cpXt_arr[tt];
+ XVY = XVY + tXYt_arr[tt];
+ }
+ Matrix<> beta_var = invpd(B0 + XVX/sigma_store(iter, j));
+ Matrix<> beta_mean = beta_var*(B0*b0 + XVY/sigma_store(iter, j));
+ beta(j,_) = stream.rmvnorm(beta_mean, beta_var);
+
+ if (K == 1){
+ density_beta(iter, j) = ::exp(log(dnorm(beta_st(j), beta_mean[0], sqrt(beta_var[0]))));
+ }
+ else{
+ density_beta(iter, j) = ::exp(lndmvn(::t(beta_st(j,_)), beta_mean, beta_var));
+ }
+ }
+ }
+ double pdf_beta = log(prod(meanc(density_beta)));
+
+ //////////////////////////////////////////////////////////////////
+ // 2. pdf.sigma2
+ //////////////////////////////////////////////////////////////////
+ Matrix<> density_sigma2(nstore, ns);
+ for (int iter = 0; iter<nstore; ++iter){
+ Matrix<> F(ntime, ns);
+ Matrix<> pr1(ns, 1);
+ pr1[0] = 1;
+ Matrix<> py(ns, 1);
+ Matrix<> pstyt1(ns, 1);
+
+ for (int tt=0; tt<ntime ; ++tt) {
+ int nsubj_s = time_groupinfo(tt, 2);
+ for (int j=0; j<ns;++j){
+ Matrix<> Sigma = eye(nsubj_s)*sigma2[j];
+ Matrix<> Mu = Xt_arr[tt]*::t(beta_st(j,_));
+ py[j] = ::exp(lndmvn(Yt_arr[tt], Mu, Sigma));
+ }
+ if (tt==0) pstyt1 = pr1;
+ else {
+ pstyt1 = ::t(F(tt-1,_)*P);
+ }
+
+ Matrix<> unnorm_pstyt = pstyt1%py;
+ Matrix<> pstyt = unnorm_pstyt/sum(unnorm_pstyt);
+ for (int j=0; j<ns ; ++j) {
+ F(tt,j) = pstyt(j);
+ }
+ }// end of F matrix filtering
+ Matrix<int> state(ntime, 1);
+ Matrix<> ps = Matrix<>(ntime, ns);
+ ps(ntime-1,_) = F(ntime-1,_);
+ state(ntime-1) = ns;
+ Matrix<> pstyn = Matrix<>(ns, 1);
+ double pone = 0.0;
+ int tt = ntime-2;
+ while (tt >= 0){
+ int st = state(tt+1);
+ Matrix<> Pst_1 = ::t(P(_, st-1));
+ Matrix<> unnorm_pstyn = F(tt,_)%Pst_1;
+ pstyn = unnorm_pstyn/sum(unnorm_pstyn);
+ if (st==1){
+ state(tt) = 1;
+ }
+ else{
+ pone = pstyn(st-2);
+ if(stream.runif() < pone) state(tt) = st-1;
+ else state(tt) = st;
+ }
+ ps(tt,_) = pstyn;
+ --tt;
+
+ }// end of while loop
+
+ /////////////////////////////////////////////////
+ // 2.2. sigma2| beta_st, s, P
+ /////////////////////////////////////////////////
+ Matrix<int> nstate(ns, 1);
+ for (int j = 0; j<ns; ++j){
+ for (int i = 0; i<ntime; ++i){
+ if (state(i) == j + 1) {
+ nstate[j] = nstate[j] + 1;
+ }// end of if
+ }// end of int i<n
+ }
+
+ Matrix<> SSE(ns, 1);
+ int beta_count = 0;
+ Matrix<int> YN(ns, 1);
+ for (int j = 0; j<ns; ++j){
+ beta_count = beta_count + nstate[j];
+ for(int s = 0; s<nsubj; ++ s) {
+ Matrix<> yj = Yk_arr[s]((beta_count - nstate[j]), 0, (beta_count - 1), 0);
+ Matrix<> Xj = Xk_arr[s]((beta_count - nstate[j]), 0, (beta_count - 1), K-1);
+ YN[j] = YN[j] + yj.rows();
+ Matrix<> e = ::t(yj - Xj*::t(beta_st(j,_)))*(yj - Xj*::t(beta_st(j,_)));
+ SSE[j] = SSE[j] + e[0];
+ }
+ double nu = (c0 + (double)YN[j])/2;
+ double scale = (d0 + SSE[j])/2;
+ sigma2[j] = stream.rigamma(nu, scale);
+
+ density_sigma2(iter, j) = ::exp(ln_invgamma(sigma2_st[j], nu, scale));
+ }
+
+ /////////////////////////////////////////////////
+ // 2.3. P| beta_st, sigma2, s
+ /////////////////////////////////////////////////
+ double shape1 = 0;
+ double shape2 = 0;
+ for (int j =0; j<m; ++j){
+ shape1 = std::abs(P0(j,j) + nstate[j] - 1);
+ shape2 = P0(j,j+1) + 1; //
+ P(j,j) = stream.rbeta(shape1, shape2);
+ P(j,j+1) = 1 - P(j,j);
+ }
+ P(m, m) = 1;
+ }
+ double pdf_sigma2 = log(prod(meanc(density_sigma2)));
+
+ //////////////////////
+ // 3. pdf.P:
+ //////////////////////
+ Matrix<> density_P(nstore, ns);
+ // Matrix<> density_local_P(ns, 1);
+
+ for (int iter = 0; iter < nstore; ++iter){
+ Matrix<> F(ntime, ns);
+ Matrix<> pr1(ns, 1);
+ pr1[0] = 1;
+ Matrix<> py(ns, 1);
+ Matrix<> pstyt1(ns, 1);
+ for (int tt=0; tt<ntime ; ++tt) {
+ if(nobst[tt]>0){
+ int nsubj_s = time_groupinfo(tt, 2);
+ for (int j=0; j<ns;++j){
+ Matrix<> Sigma = eye(nsubj_s)*sigma2_st[j];
+ Matrix<> Mu = Xt_arr[tt]*::t(beta_st(j,_));
+ py[j] = ::exp(lndmvn(Yt_arr[tt], Mu, Sigma));
+ }
+ if (tt==0) pstyt1 = pr1;
+ else {
+ pstyt1 = ::t(F(tt-1,_)*P);
+ }
+ Matrix<> unnorm_pstyt = pstyt1%py;
+ Matrix<> pstyt = unnorm_pstyt/sum(unnorm_pstyt);
+ for (int j=0; j<ns ; ++j) {
+ F(tt,j) = pstyt(j);
+ }
+ }
+ }// end of F matrix filtering
+ Matrix<int> state(ntime, 1);
+ Matrix<> ps = Matrix<>(ntime, ns);
+ ps(ntime-1,_) = F(ntime-1,_);
+ state(ntime-1) = ns;
+ Matrix<> pstyn = Matrix<>(ns, 1);
+ double pone = 0.0;
+ int tt = ntime-2;
+
+ while (tt >= 0){
+ if(nobst[tt]>0){
+ int st = state(tt+1);
+ Matrix<> Pst_1 = ::t(P(_, st-1));
+ Matrix<> unnorm_pstyn = F(tt,_)%Pst_1;
+ pstyn = unnorm_pstyn/sum(unnorm_pstyn);
+ if (st==1){
+ state(tt) = 1;
+ }
+ else{
+ pone = pstyn(st-2);
+ if(stream.runif() < pone) state(tt) = st-1;
+ else state(tt) = st;
+ }
+ ps(tt,_) = pstyn;
+ --tt;
+ }
+ }// end of while loop
+
+ Matrix<int> nstate(ns, 1);
+ for (int j = 0; j<ns; ++j){
+ for (int i = 0; i<ntime; ++i){
+ if (state(i) == j + 1) {
+ nstate[j] = nstate[j] + 1;
+ }// end of if
+ }// end of int i<n
+ }
+ double shape1 = 0;
+ double shape2 = 0;
+ for (int j =0; j<m; ++j){
+ shape1 = std::abs(P0(j,j) + nstate[j] - 1);
+ shape2 = P0(j,j+1) + 1; //
+ P(j,j) = stream.rbeta(shape1, shape2);
+ P(j,j+1) = 1 - P(j,j);
+ density_P(iter, j) = dbeta(P_st(j,j), shape1, shape2);
+ }
+ density_P(iter, ns-1) = 1;
+ }
+
+ double pdf_P = log(prod(meanc(density_P)));
+
+
+ //////////////////////////////////////////////////////////////////
+ // likelihood f(y|beta_st, D_st, sigma2_st, P_st)
+ //////////////////////////////////////////////////////////////////
+ Matrix<> F(ntime, ns);
+ Matrix<> pr1(ns, 1);
+ Matrix<> like(ntime, 1);
+ pr1[0] = 1;
+ Matrix<> py(ns, 1);
+ Matrix<> pstyt1(ns, 1);
+ for (int tt=0; tt<ntime ; ++tt) {
+ for (int j=0; j<ns;++j){
+ int nsubj_s = time_groupinfo(tt, 2);
+ Matrix<> Sigma = eye(nsubj)*sigma2_st[j];
+ Matrix<> Mu = Xt_arr[tt]*::t(beta_st(j,_));
+ py[j] = ::exp(lndmvn(Yt_arr[tt], Mu, Sigma));
+ }
+ if (tt==0) pstyt1 = pr1;
+ else {
+ pstyt1 = ::t(F(tt-1,_)*P_st);
+ }
+ Matrix<> unnorm_pstyt = pstyt1%py;
+ Matrix<> pstyt = unnorm_pstyt/sum(unnorm_pstyt);
+ for (int j=0; j<ns ; ++j) {
+ F(tt,j) = pstyt(j);
+ }
+ like[tt] = sum(unnorm_pstyt);
+ }// end of F matrix filtering
+
+ loglike = sum(log(like));
+
+ //////////////////////
+ // log prior ordinate
+ //////////////////////
+ Matrix<double> density_beta_prior(ns, 1);
+ Matrix<double> density_sigma2_prior(ns, 1);
+ Matrix<double> density_P_prior(ns, 1);
+ density_P_prior[ns-1] = 0; //
+
+ for (int j=0; j<ns ; ++j){
+ if (K == 1){
+ density_beta_prior[j] = log(dnorm(beta_st(j), b0[0], sqrt(B0inv[0])));
+ }
+ else{
+ density_beta_prior[j] = lndmvn(::t(beta_st(j,_)), b0, B0inv);
+ }
+ density_sigma2_prior[j] = ln_invgamma(sigma2_st[j], c0/2, d0/2);
+ }
+ for (int j=0; j<m ; ++j){
+ density_P_prior[j] = log(dbeta(P_st(j,j), P0(j,j), P0(j,j+1)));
+ }
+
+ // compute marginal likelihood
+ double logprior = sum(density_beta_prior) + sum(density_sigma2_prior) +
+ sum(density_P_prior);
+ logmarglike = (loglike + logprior) - (pdf_beta + pdf_P + pdf_sigma2);
+ if(verbose > 0){
+ Rprintf("\nlogmarglike = %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_sigma2);
+ Rprintf("pdf_P = %10.5f\n", pdf_P);
+ }
+
+ }// end of marginal likelihood computation
+
+ delete[] nobst;
+ for(int k=0; k<ntime; k++) {
+ delete[] post_arr[k];
+ }
+ delete[] nobsk;
+ for(int k=0; k<nsubj; k++) {
+ delete[] posk_arr[k];
+ }
+ delete [] posk_arr;
+ delete [] Xk_arr;
+ delete [] Yk_arr;
+ delete [] Xt_arr;
+ delete [] tXt_arr;
+ delete [] Yt_arr;
+ delete [] cpXt_arr;
+ delete [] tXYt_arr;
+
+
+}// end of HMMmultivariateGaussian_impl
+
+extern "C" {
+ void HMMmultivariateGaussian(double* betadata, const int* betarow, const int* betacol,
+ double* sigmadata, double *psout,
+ const int* nsubj, const int* ntime,
+ const int* m, const int* nobs, const int* subjectid, const int* timeid,
+ const double* Ydata, const int* Yrow, const int* Ycol,
+ const double* Xdata, const int* Xrow, const int* Xcol,
+ const double* YTdata, const double* XTdata,
+ const int* burnin, const int* mcmc, const int* thin, const int* verbose,
+ const int *uselecuyer, const int *seedarray, const int *lecuyerstream,
+ const double* betastartdata, const double* sigma2start,
+ const double* b0data, const double* B0data,
+ const double* c0, const double* d0,
+ const double* P0data, const int* P0row, const int* P0col,
+ const double* subject_groupinfodata, const double* time_groupinfodata,
+ double *logmarglikeholder, double *loglikeholder,
+ const int *chib){
+ // pull together Matrix objects
+ Matrix<> Y(*Yrow, *Ycol, Ydata);
+ Matrix<> X(*Xrow, *Xcol, Xdata);
+ Matrix<> YT(*Yrow, *Ycol, YTdata);
+ Matrix<> XT(*Xrow, *Xcol, XTdata);
+ Matrix<> betastart(*Xcol, 1, betastartdata);
+ Matrix<> b0(*Xcol, 1, b0data);
+ Matrix<> B0(*Xcol, *Xcol, B0data);
+ Matrix<> Pstart(*P0row, *P0col, P0data);
+ Matrix<int> subjectid_mat(*nobs, 1, subjectid);
+ Matrix<int> timeid_mat(*nobs, 1, timeid);
+ Matrix<> subject_groupinfo(*nsubj, 3, subject_groupinfodata);
+ Matrix<> time_groupinfo(*ntime, 3, time_groupinfodata);
+ const int ns = *m + 1;
+ double logmarglike;
+ double loglike;
+
+ if (*m == 0) {
+ Matrix<> beta_store(*betarow, *betacol);
+ Matrix<> sigma_store(*betarow, 1);
+ MCMCPACK_PASSRNG2MODEL(MultivariateGaussian_impl,
+ *nsubj, *ntime, *nobs,
+ subjectid_mat, timeid_mat,
+ Y, X, YT, XT,
+ *burnin, *mcmc, *thin, *verbose,
+ *sigma2start,
+ b0, B0, *c0, *d0,
+ time_groupinfo, subject_groupinfo,
+ beta_store, sigma_store,
+ logmarglike, loglike, *chib);
+ // store marginal likelihood
+ logmarglikeholder[0] = logmarglike;
+ loglikeholder[0] = loglike;
+
+ for (int i=0; i < (*betarow* *betacol); ++i){
+ betadata[i] = beta_store(i);
+ }
+ for (int i=0; i < (*betarow*ns); ++i){
+ sigmadata[i] = sigma_store(i);
+ }
+ }// end of if m == 0
+ else {
+ Matrix<> beta_store(*betarow, *betacol);
+ Matrix<> sigma_store(*betarow, ns);
+ Matrix<> P_store(*betarow, ns*ns);
+ Matrix<> s_store(*betarow, *ntime);
+ Matrix<> ps_store(*ntime, ns);
+ MCMCPACK_PASSRNG2MODEL(HMMmultivariateGaussian_impl,
+ *nsubj, *ntime, *m, *nobs,
+ subjectid_mat, timeid_mat,
+ Y, X, YT, XT,
+ *burnin, *mcmc, *thin, *verbose,
+ betastart, *sigma2start,
+ b0, B0, *c0, *d0,
+ Pstart, time_groupinfo, subject_groupinfo,
+ beta_store, sigma_store, P_store, ps_store, s_store,
+ logmarglike, loglike, *chib);
+ // store marginal likelihood
+ logmarglikeholder[0] = logmarglike;
+ loglikeholder[0] = loglike;
+
+ for (int i=0; i < (*betarow* *betacol); ++i){
+ betadata[i] = beta_store(i);
+ }
+ for (int i=0; i < (*betarow*ns); ++i){
+ sigmadata[i] = sigma_store(i);
+ }
+ for (int i = 0; i<(*ntime *ns); ++i){
+ psout[i] = ps_store[i];
+ }
+
+ }// end of m>0
+ }// end of HMMmultivariateGaussian_CC
+}// end of extern C
+
+#endif /*HMMmultivariateGaussian_CC */
diff --git a/src/HMMpanelFE.cc b/src/HMMpanelFE.cc
new file mode 100644
index 0000000..c48f889
--- /dev/null
+++ b/src/HMMpanelFE.cc
@@ -0,0 +1,434 @@
+//////////////////////////////////////////////////////////////////////////
+// HMMpanelFE.cc is C++ code
+//
+// Written by Jong Hee Park
+// 11/19/2008
+// 09/20/2009
+// Copyright (C) 2008-present Jong Hee Park
+//////////////////////////////////////////////////////////////////////////
+#ifndef HMMPANELFE_CC
+#define HMMPANELFE_CC
+
+#include<vector>
+#include<algorithm>
+
+#include "MCMCrng.h"
+#include "MCMCfcds.h"
+#include "matrix.h"
+#include "distributions.h"
+#include "stat.h"
+#include "la.h"
+#include "ide.h"
+#include "smath.h"
+#include "rng.h"
+#include "mersenne.h"
+#include "lecuyer.h"
+#include "lapack.h"
+
+#include <R.h> // needed to use Rprintf()
+#include <R_ext/Utils.h> // needed to allow user interrupts
+
+using namespace std;
+using namespace scythe;
+
+
+// used to access 1d arrays holding R matrices like a 2d array
+#define M(ROW,COL,NROWS) (COL*NROWS+ROW)
+
+template <typename RNGTYPE>
+// HMM_Gaussian_Fixed_state_sampler(stream, mvector[s], ntime_s, Yres[s], Zarr[s], delta, P);
+// For better identification, first two and last two states are constrained to be 1 and ns
+// 10/31/2011 JHP
+Matrix<int> hetero_state_sampler(rng<RNGTYPE>& stream,
+ const int m,
+ const int ntime_s,
+ const Matrix<>& Y,
+ const Matrix<>& delta,
+ const Matrix<>& Sigma,
+ const Matrix<>& P){
+ const int ns = m + 1;
+ const int ntime = ntime_s;
+
+ Matrix<> F(ntime, ns);
+ Matrix<> pr1(ns, 1);
+ pr1[0] = 1;
+ Matrix<> py(ns, 1);
+ Matrix<> pstyt1(ns, 1);
+
+ for (int tt=0; tt<ntime ; ++tt){
+ for (int j = 0; j< ns; ++j){
+ py[j] = dnorm(Y[tt], delta[j], sqrt(Sigma[j]));
+ }
+ if (tt==0) pstyt1 = pr1;
+ else {
+ pstyt1 = ::t(F(tt-1,_)*P); // make it an ns by 1 matrix
+ }
+
+ Matrix<> unnorm_pstyt = pstyt1%py;
+ const Matrix<> 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
+ Matrix<int> state(ntime, 1);
+ Matrix<> ps = Matrix<>(ntime, ns);
+ ps(ntime-1,_) = F(ntime-1,_);
+ state(ntime-1) = ns;
+
+ Matrix<> pstyn = Matrix<>(ns, 1);
+ double pone = 0.0;
+ int tt = ntime-2;
+
+ while (tt >= 0){
+ int st = state(tt+1);
+ Matrix<> Pst_1 = ::t(P(_,st-1));
+ Matrix<> unnorm_pstyn = F(tt,_)%Pst_1;
+ pstyn = unnorm_pstyn/sum(unnorm_pstyn);
+
+ if (st==1){
+ state(tt) = 1;
+ }
+ else{
+ pone = pstyn(st-2);
+ if(stream.runif() < pone) state(tt) = st-1;
+ else state(tt) = st;
+ }
+ ps(tt,_) = pstyn;
+ --tt;
+ }// end of while loop
+
+ return state;
+
+} // end of state sampler
+
+
+template <typename RNGTYPE>
+void HMMpanelFE_impl(rng<RNGTYPE>& stream,
+ unsigned int nsubj, unsigned int ntime,
+ unsigned int mmax, unsigned int mmin,
+ const Matrix<int>& mvector,
+ unsigned int totalstates,
+ const Matrix<>& Y, const Matrix<>& X,
+ const Matrix<int>& subjectid,
+ unsigned int burnin, unsigned int mcmc,
+ unsigned int thin, unsigned int verbose,
+ Matrix<>& beta, double sigma2, Matrix<>& deltastart,
+ const Matrix<>& b0, const Matrix<>& B0,
+ const double delta0, const double Delta0,
+ const double c0, const double d0,
+ const Matrix<>& P0data,
+ const Matrix<>& Pstart,
+ const Matrix<>& subject_groupinfo,
+ Matrix<>& betastorage,
+ Matrix<>& statestorage,
+ Matrix<>& deltastorage,
+ Matrix<>& sigmastorage){
+
+ // redefine constants
+ const unsigned int K = X.cols();; // ncol(X)
+ const int NOBS = Y.rows();
+
+ const int tot_iter = burnin + mcmc;
+ const int nstore = mcmc / thin; // number of draws to store
+
+ vector< vector <double> > P0_holder;
+ vector< vector <double> > P_holder;
+ int count = 0;
+ for (int s=0; s< nsubj; ++s){
+ const int nms = mvector[s] + 1;
+ vector <double> P0mat;
+ vector <double> Pmat;
+ for(int ii=0; ii<(nms*nms); ++ii){
+ P0mat.push_back(P0data(count + ii));
+ Pmat.push_back(Pstart(count + ii));
+ }
+ count = count + nms*nms;
+
+ P0_holder.push_back(P0mat);
+ P_holder.push_back(Pmat);
+ }
+
+ vector< vector <double> > delta_holder;
+ vector< vector <double> > sigma2_holder;
+ vector< vector <int> > nstate;
+ count = 0;
+ for (int s=0; s< nsubj; ++s){
+ int nms = mvector[s] + 1;
+ int ntime_s = subject_groupinfo(s, 2);
+ vector <double> deltamat;
+ vector <double> sigmamat;
+ vector <int> nstatemat;
+ for(int ii=0; ii<nms ;++ii){
+ deltamat.push_back(0);
+ nstatemat.push_back(ntime_s);
+ sigmamat.push_back(sigma2);
+ }
+ count = count + nms;
+ delta_holder.push_back(deltamat);
+ nstate.push_back(nstatemat);
+ sigma2_holder.push_back(sigmamat);
+ }
+
+ int *nobsk = new int[nsubj];
+ for (int k=0; k<nsubj; k++) {
+ nobsk[k]=0;
+ for (int n=0; n<NOBS; n++) {
+ if (subjectid[n]==k+1) {
+ nobsk[k]+=1;
+ }
+ }
+ }
+
+ int **posk_arr = new int*[nsubj];
+ for (int k=0; k<nsubj; k++) {
+ posk_arr[k] = new int[nobsk[k]];
+ int repk=0;
+ for (int n=0; n<NOBS; n++) {
+ if (subjectid[n]==k+1) {
+ posk_arr[k][repk]=n;
+ repk++;
+ }
+ }
+ }
+
+ Matrix<double> *Yk_arr = new Matrix<double>[nsubj];
+ Matrix<double> *Xk_arr = new Matrix<double>[nsubj];
+ for(int k=0; k<nsubj; k++) {
+ Xk_arr[k] = Matrix<double>(nobsk[k], K);
+ Yk_arr[k] = Matrix<double>(nobsk[k], 1);
+ for (int l=0; l<nobsk[k]; l++) {
+ for (int p=0; p< K; p++) {
+ Xk_arr[k](l, p) = X[p*NOBS + posk_arr[k][l]];
+ }
+ Yk_arr[k](l,0) = Y[posk_arr[k][l]];
+ }
+ }
+
+ Matrix<double> *tXk_arr = new Matrix<double>[nsubj];
+ for(int k=0; k<nsubj; k++) {
+ tXk_arr[k] = ::t(Xk_arr[k]);
+ }
+
+ /////////////////////////////////////////////////
+ // initialize newY for the first loop only
+ /////////////////////////////////////////////////
+ Matrix<>* newY = new Matrix<>[nsubj]; // newY = Y - delta
+ Matrix<>* Yres = new Matrix<>[nsubj]; // Yres = Y - Xbeta
+ for (int s=0; s<nsubj; ++s){
+ newY[s] = Yk_arr[s] - delta0;
+ Yres[s] = Yk_arr[s];
+ }
+
+ // MCMC iterations start here
+ int sampcount = 0;
+ for (int iter=0; iter < tot_iter; ++iter){
+ double delta_sum = 0;
+ for (int s=0; s< nsubj; ++s) {
+ int ntime_s = subject_groupinfo(s, 2);
+ Matrix<> Zk_arr(ntime_s, 1);
+ for (int tt=0; tt< ntime_s; ++tt) {
+ Zk_arr(tt) = 1;
+ }
+ Yres[s] = Yk_arr[s] - Xk_arr[s]*beta;
+
+ if(mvector[s]==0){
+ double Dn = 1/(Delta0 + (double)ntime_s/sigma2_holder[s][0]);
+ double dn = Dn*(Delta0*delta0 + sum(Yres[s])/sigma2_holder[s][0]);
+ delta_holder[s][0] = stream.rnorm(dn, sqrt(Dn));
+ delta_sum = delta_sum + delta_holder[s][0];
+ newY[s] = Yk_arr[s] - Zk_arr*delta_holder[s][0];
+
+ // Sample sigma
+ double shape = (c0 + (double)ntime_s)/2;
+ const Matrix<> SSE = crossprod (newY[s]);
+ double scale =(d0 + SSE[0])/2;
+ sigma2_holder[s][0] = 1/stream.rgamma(shape, scale);
+ }
+ else {
+ const int nscur = mvector[s] + 1;
+ Matrix<> P(nscur, nscur);
+ Matrix<> P0(nscur, nscur);
+ Matrix<> delta(nscur, 1);
+ Matrix<> Sigma(nscur, 1);
+
+ for (int i=0;i<(nscur*nscur); ++i){
+ P0[i] = P0_holder[s][i];
+ P[i] = P_holder[s][i];
+ }
+
+ for (int i=0;i<nscur; ++i){
+ delta[i] = delta_holder[s][i];
+ Sigma[i] = sigma2_holder[s][i];
+ }
+ // Sample s
+ Matrix<int> state_s = hetero_state_sampler(stream, mvector[s], ntime_s, Yres[s], delta, Sigma, P);
+ // Sample delta and Sigma
+ int delta_count = 0;
+ for (int j = 0; j <nscur ; ++j){
+ nstate[s][j] = 0;
+ for (int i = 0; i<ntime_s; ++i){
+ if (state_s[i] == (j+1)) {
+ nstate[s][j] = nstate[s][j] + 1;
+ }// end of if
+ }// end of int i<n
+
+ delta_count = delta_count + nstate[s][j];
+ // sample delta
+ Matrix<> yj = Yres[s]((delta_count - nstate[s][j]), 0, (delta_count - 1), 0);
+ Matrix<> Yj = Yk_arr[s]((delta_count - nstate[s][j]), 0, (delta_count - 1), 0);
+ double Dn = 1/(Delta0 + (double)nstate[s][j]/sigma2_holder[s][j]);
+ double dn = Dn*(Delta0*delta0 + sum(yj)/sigma2_holder[s][j]);
+ delta_holder[s][j] = stream.rnorm(dn, sqrt(Dn));
+ delta_sum = delta_sum + delta_holder[s][j];
+ newY[s]((delta_count - nstate[s][j]), 0, (delta_count - 1), 0) = Yj - delta_holder[s][j];
+
+ // Sample sigma
+ double shape = (c0 + (double)nstate[s][j])/2;
+ const Matrix<> SSE = crossprod (newY[s]((delta_count - nstate[s][j]), 0, (delta_count - 1), 0));
+ double scale =(d0 + SSE[0])/2;
+ sigma2_holder[s][j] = 1/stream.rgamma(shape, scale);
+
+ }
+ // assure that there is no label switching problem
+ // the code needs to be added here
+
+ // Sample P
+ double shape1 = 0;
+ double shape2 = 0;
+
+ for (int j =0; j<(nscur-1); ++j){
+ shape1 = std::abs(P0(j,j) + nstate[s][j] - 1);
+ shape2 = P0(j,j+1) + 1; //
+ P(j,j) = stream.rbeta(shape1, shape2);
+ P(j,j+1) = 1 - P(j,j);
+ }
+ P(mvector[s], mvector[s]) = 1; //no jump at the last state
+
+ for(int ii=0; ii<(nscur*nscur) ;++ii) {
+ P_holder[s][ii] = P[ii];
+ }
+ }//end of else (mvector!=0)
+
+ }// end of subject specific looping
+ // Sample beta
+ Matrix<> XVX(K, K);
+ Matrix<> XVY(K, 1);
+ for(int s = 0; s<nsubj; ++s) {
+ int ntime_s = subject_groupinfo(s, 2);
+ int delta_count = 0;
+ Matrix<> Vi = eye(ntime_s);
+ for (int j = 0; j <(mvector[s] + 1); ++j){
+ delta_count = delta_count + nstate[s][j];
+ for(int i = (delta_count - nstate[s][j]); i<delta_count; ++i) {
+ Vi(i,i) = 1/sigma2_holder[s][j];
+ }
+ }
+ XVX = XVX + tXk_arr[s] * Vi * Xk_arr[s];
+ XVY = XVY + tXk_arr[s] * Vi * newY[s];
+ }
+ Matrix<> beta_var = invpd(B0 + XVX);
+ Matrix<> beta_mean = beta_var*(B0*b0 + XVY);
+ beta = stream.rmvnorm(beta_mean, beta_var);
+ // STORE
+ if (iter >= burnin && ((iter % thin) == 0)) {
+ for(int j=0;j<K; ++j) {
+ betastorage(sampcount,j) = beta(j);
+ }
+ int count = 0;
+ for (int s=0; s<nsubj; ++s){
+ for (int j=0; j<(mvector[s] + 1); ++j){
+ sigmastorage(sampcount, count) = sigma2_holder[s][j];
+ deltastorage(sampcount, count) = delta_holder[s][j];
+ double nstat = nstate[s][j];
+ statestorage(sampcount, count) = nstat;
+ ++ count;
+ }
+ }
+ ++sampcount;
+ }
+
+ // REPORT
+ if(verbose > 0 && iter % verbose == 0){
+ Rprintf("\n ----------------------------------------------------------------------- ");
+ Rprintf("\n\n HMMpanelFE %i of %i \n", iter, tot_iter);
+ Rprintf("\n beta = \n");
+ for(int i=0;i<K; ++i) {
+ Rprintf("%10.5f\n", beta(i));
+ }
+ }
+ }// END of MCMC loop
+
+ delete [] Yres;
+ delete [] tXk_arr;
+ delete [] Xk_arr;
+ delete [] Yk_arr;
+ delete [] newY;
+}
+
+
+extern "C" {
+ void HMMpanelFE(double *deltadraws, double* sigmadraws,
+ double *statedraws, //const int* statecol,
+ double* betadraws, const int* betarow, const int* betacol,
+ const int* totalstates,
+ const int* nsubj, const int* ntime, const int* nobs,
+ const int* subjectid,
+ const int* m,
+ const int* mmax, const int* mmin,
+ const double* Ydata, const int* Yrow, const int* Ycol,
+ const double* Xdata, const int* Xrow, const int* Xcol,
+ const int* burnin, const int* mcmc, const int* thin, const int* verbose,
+ const int *uselecuyer, const int *seedarray, const int *lecuyerstream,
+ const double* betastartdata, const double* sigma2start,
+ const double* deltastartdata, const int* deltastartrow,
+ const double* b0data, const double* B0data,
+ const double* delta0, const double* Delta0,
+ const double* c0, const double* d0,
+ const double* P0data, const int* P0row,
+ const double* Pstartdata,
+ const double* subject_groupinfodata){
+
+
+ // pull together Matrix objects
+ Matrix<> Y(*Yrow, *Ycol, Ydata);
+ Matrix<> X(*Xrow, *Xcol, Xdata);
+ Matrix<> betastart(*Xcol, 1, betastartdata);
+ Matrix<> deltastart(*deltastartrow, 1, deltastartdata);
+ Matrix<> b0(*Xcol, 1, b0data);
+ Matrix<> B0(*Xcol, *Xcol, B0data);
+ Matrix<int> subjectid_mat(*nobs, 1, subjectid);
+ Matrix<> subject_groupinfo(*nsubj, 3, subject_groupinfodata);
+ Matrix<> P0(*P0row, 1, P0data);
+ Matrix<> Pstart(*P0row, 1, Pstartdata);
+ Matrix<int> mvector(*nsubj, 1, m);
+
+ Matrix<> betastorage(*betarow, *betacol);
+ Matrix<> sigmastorage(*betarow, *totalstates);
+ Matrix<> deltastorage(*betarow, *totalstates);
+ Matrix<> statestorage(*betarow, *totalstates);
+
+ MCMCPACK_PASSRNG2MODEL(HMMpanelFE_impl,
+ *nsubj, *ntime, *mmax, *mmin, mvector,
+ *totalstates, Y, X,
+ subjectid_mat,
+ *burnin, *mcmc, *thin, *verbose,
+ betastart, *sigma2start, deltastart,
+ b0, B0, *delta0, *Delta0, *c0, *d0,
+ P0, Pstart, subject_groupinfo,
+ betastorage, statestorage, deltastorage, sigmastorage);
+ unsigned int deltasize = *betarow * *totalstates;
+ for (int i=0; i < deltasize; ++i){
+ deltadraws[i] = deltastorage(i);
+ sigmadraws[i] = sigmastorage(i);
+ statedraws[i] = statestorage(i);
+ }
+ unsigned int betasize = *betarow * *betacol;
+ for (int i=0; i < betasize; ++i){
+ betadraws[i] = betastorage(i);
+ }
+ }// end of HMMpanelFE
+
+}// end of extern C
+
+#endif /* HMMPANELFE_CC */
diff --git a/src/HMMpanelRE.cc b/src/HMMpanelRE.cc
new file mode 100644
index 0000000..22872d6
--- /dev/null
+++ b/src/HMMpanelRE.cc
@@ -0,0 +1,1413 @@
+//////////////////////////////////////////////////////////////////////////
+// HMMpanelRE.cc is C++ code to estimate a Gaussian panel model with a structural break
+// y_{it} = \x'_{it}\b + \w'_{it}\bi_i + \z'_{it}\d_{s_{it}} +\varepsilon_{it}
+// \varepsilon_{it} \sim \normdist{0}{\sigma^2}
+// \bi_i \sim \normdist{0}{\D}
+// Parameters = {beta, bi, D, sigma2, || {s_it}, delta, P_i}// static + changing
+//
+// First written by Jong Hee Park on 11/19/2008
+// Modified by JHP (consistent with MCMChregress by Ghislain Vieilledent) on 07/04/2011
+// Thanks for letting know the Woodbury matrix identity!
+// Included in MCMCpack b JHP on 09/2011
+//////////////////////////////////////////////////////////////////////////
+#ifndef HMMPANELRE_CC
+#define HMMPANELRE_CC
+
+#include<vector>
+#include<algorithm>
+
+#include "MCMCrng.h"
+#include "MCMCfcds.h"
+#include "matrix.h"
+#include "distributions.h"
+#include "stat.h"
+#include "la.h"
+#include "ide.h"
+#include "smath.h"
+#include "rng.h"
+#include "mersenne.h"
+#include "lecuyer.h"
+#include "lapack.h"
+
+#include <R.h> // needed to use Rprintf()
+#include <R_ext/Utils.h> // needed to allow user interrupts
+
+using namespace std;
+using namespace scythe;
+
+// used to access 1d arrays holding R matrices like a 2d array
+#define M(ROW,COL,NROWS) (COL*NROWS+ROW)
+
+static double lndwish (const Matrix<>& W,
+ unsigned int v,
+ const Matrix<>& S) {
+ const int K = S.rows();
+ double gammapart = 1;
+ for (int i=0; i<K; ++i){
+ double gammain = (v - i)/2;
+ gammapart = gammapart * gammafn(gammain);// since i starts from 0!
+ }
+ double logdenom = log(gammapart) + (v * K/2)*log(2) + (K * (K - 1)/4)*log(M_PI);
+ double detS = 0;
+ double detW = 0;
+ if(K==1){
+ detS = S[0];
+ detW = W[0];
+ }
+ else{
+ detS = det(S);
+ detW = det(W);
+ }
+ Matrix<> hold = invpd(S) * W;
+ Matrix<> diaghold(K, 1);
+ diaghold(_,0) = diag(hold);
+ double tracehold = sum(diaghold);
+ double lognum = -1*(v *.5)*log(detS) + (v - K - 1)/2*log(detW) - 1/2 * tracehold;
+ return(lognum - logdenom);
+}
+
+static double lndinvgamma(double theta, double a, double b) {
+ double logf = a * log(b) - lngammafn(a) + -(a+1) * log(theta) +
+ -b/theta;
+ return (logf);
+ //pow(b, a) / gammafn(a) * pow(theta, -(a+1)) * exp(-b/theta);
+}
+
+template <typename RNGTYPE>
+void GaussianPanelRE_impl(rng<RNGTYPE>& stream,
+ unsigned int nsubj, unsigned int ntime,
+ unsigned int nobs,
+ const Matrix<int>& subjectid, const Matrix<int>& timeid,
+ const Matrix<>& Y, const Matrix<>& X, const Matrix<>& W,
+ unsigned int burnin, unsigned int mcmc,
+ unsigned int thin, unsigned int verbose,
+ const double sigma2start, Matrix<>& D,
+ const Matrix<>& b0, const Matrix<>& B0,
+ const double c0, const double d0,
+ unsigned int r0, const Matrix<>& R0,
+ const Matrix<>& time_groupinfo,
+ const Matrix<>& subject_groupinfo,
+ Matrix<>& beta_store, Matrix<>& sigma_store,
+ Matrix<>& D_store,
+ double& logmarglike, double& loglike,
+ unsigned int chib)
+{ // redefine constants
+ const unsigned int K = X.cols();; // ncol(X)
+ const unsigned int Q = W.cols(); // ncol(W)
+ const int NOBS = nobs;
+ double sigma2 = sigma2start;
+
+ const Matrix<> R0inv = invpd(R0);
+ const Matrix<> B0inv = invpd(B0);
+
+ const int tot_iter = burnin + mcmc;
+ const int nstore = mcmc / thin; // number of draws to store
+
+ Matrix<double> Dinv = invpd(D);
+
+ // generate posk_arr and post_arr
+ // Number of observations by group k
+ int *nobsk = new int[nsubj];
+ for (int k=0; k<nsubj; k++) {
+ nobsk[k]=0;
+ for (int n=0; n<NOBS; n++) {
+ if (subjectid[n]==k+1) {
+ nobsk[k]+=1;
+ }
+ }
+ }
+
+ // Position of each group in the data-set
+ int **posk_arr = new int*[nsubj];
+ for (int k=0; k<nsubj; k++) {
+ posk_arr[k] = new int[nobsk[k]];
+ int repk=0;
+ for (int n=0; n<NOBS; n++) {
+ if (subjectid[n]==k+1) {
+ posk_arr[k][repk]=n;
+ repk++;
+ }
+ }
+ }
+
+ // Small fixed matrices indexed on k for data access
+ Matrix<double> *Yk_arr = new Matrix<double>[nsubj];
+ Matrix<double> *Xk_arr = new Matrix<double>[nsubj];
+ Matrix<double> *Wk_arr = new Matrix<double>[nsubj];
+ for(int k=0; k<nsubj; k++) {
+ Xk_arr[k] = Matrix<double>(nobsk[k], K);
+ Wk_arr[k] = Matrix<double>(nobsk[k], Q);
+ Yk_arr[k] = Matrix<double>(nobsk[k], 1);
+ for (int l=0; l<nobsk[k]; l++) {
+ for (int p=0; p< K; p++) {
+ Xk_arr[k](l, p) = X[p*NOBS + posk_arr[k][l]];
+ // Rprintf("\n Xk_arr ( %i, %i) is %10.5f", l, p, Xk_arr[k](l, p));
+ }
+ for (int q=0; q < Q; q++) {
+ Wk_arr[k](l, q) = W[q*NOBS + posk_arr[k][l]];
+ }
+ Yk_arr[k](l,0) = Y[posk_arr[k][l]];
+ }
+ }
+
+ Matrix<double> *tXk_arr = new Matrix<double>[nsubj];
+ Matrix<double> *tWk_arr = new Matrix<double>[nsubj];
+ Matrix<double> *tXWk_arr = new Matrix<double>[nsubj];
+ Matrix<double> *tWXk_arr = new Matrix<double>[nsubj];
+ Matrix<double> *tWYk_arr = new Matrix<double>[nsubj];
+ Matrix<double> *tXYk_arr = new Matrix<double>[nsubj];
+ Matrix<double> *cpXk_arr = new Matrix<double>[nsubj];
+ Matrix<double> *cpWk_arr = new Matrix<double>[nsubj];
+ for(int k=0; k<nsubj; k++) {
+ tXk_arr[k] = t(Xk_arr[k]);
+ tXYk_arr[k] = t(Xk_arr[k])*Yk_arr[k];
+ tXWk_arr[k] = t(Xk_arr[k])*Wk_arr[k];
+ tWk_arr[k] = t(Wk_arr[k]);
+ tWXk_arr[k] = t(Wk_arr[k])*Xk_arr[k];
+ tWYk_arr[k] = t(Wk_arr[k])*Yk_arr[k];
+ cpWk_arr[k] = crossprod(Wk_arr[k]);
+ cpXk_arr[k] = crossprod(Xk_arr[k]);
+ }
+ // MCMC iterations start here
+ int sampcount = 0;
+
+ Rprintf("\n ///////////////////////////////////////////////// \n");
+ Rprintf("\n MCMC for HMM Gaussian Panel Randome Effects loop starts! \n");
+ Rprintf("\n ///////////////////////////////////////////////// \n");
+
+ Matrix<> beta(K, 1);
+
+ /////////////////////////////////////////////////
+ // initialize Yhat for the first loop only
+ /////////////////////////////////////////////////
+ for (int iter=0; iter < tot_iter; ++iter){
+ /////////////////////////////////////////////////
+ // Step 1. Sample beta (fixed-effects coef)
+ /////////////////////////////////////////////////
+ Matrix<double> XVX(K, K);
+ Matrix<double> XVY(K, 1);
+ for(int s = 0; s<nsubj; ++s) {
+ XVX += (1/sigma2)*cpXk_arr[s]-pow(1/sigma2,2)*tXWk_arr[s]*invpd(Dinv + cpWk_arr[s]*(1/sigma2))*tWXk_arr[s];
+ XVY += (1/sigma2)*tXYk_arr[s]-pow(1/sigma2,2)*tXWk_arr[s]*invpd(Dinv + cpWk_arr[s]*(1/sigma2))*tWYk_arr[s];
+ }
+ Matrix<> beta_var = invpd(B0 + XVX/sigma2);
+ Matrix<> beta_mean = beta_var*(B0*b0 + XVY/sigma2);
+ beta = stream.rmvnorm(beta_mean, beta_var);
+
+ /////////////////////////////////////////////////
+ // Step 2. Sample bi (random-effects coef)
+ /////////////////////////////////////////////////
+ Matrix<double> bi(Q, nsubj);
+ for(int s = 0; s<nsubj; ++ s) {
+ Matrix<double> b_var = invpd(Dinv + cpWk_arr[s]/sigma2);
+ Matrix<double> b_mean = b_var*tWk_arr[s]*(Yk_arr[s] - Xk_arr[s]*beta)/sigma2;
+ bi(_,s) = stream.rmvnorm(b_mean, b_var);
+ }
+
+ /////////////////////////////////////////////////
+ // Step 3. Sample sigma2
+ /////////////////////////////////////////////////
+ double SSE = 0;
+ for(int s=0;s<nsubj; ++s){
+ int ntime_s = subject_groupinfo(s, 2);
+ Matrix<> e = t(Yk_arr[s]-Xk_arr[s]*beta - Wk_arr[s]*bi(_,s))*
+ (Yk_arr[s] - Xk_arr[s]*beta - Wk_arr[s]*bi(_,s));
+ SSE = SSE + e[0];
+ }
+ double nu = (c0 + NOBS)/2;
+ double delta = (d0 + SSE)/2;
+ sigma2 = stream.rigamma(nu, delta);
+
+ /////////////////////////////////////////////////
+ // Step 4. Sample D
+ /////////////////////////////////////////////////
+ Matrix<double> SSB = bi*t(bi);
+ Matrix<double> D_scale = invpd(R0inv + SSB);
+ int D_dof = r0 + nsubj;
+ Dinv = stream.rwish(D_dof, D_scale);
+ D = inv(Dinv);
+
+ /////////////////////////////////////////////////
+ // STORE
+ /////////////////////////////////////////////////
+ if (iter >= burnin && ((iter % thin) == 0)) {
+ for(int j=0;j<K; ++j) {
+ beta_store(sampcount,j) = beta(j);
+ }
+ for(int j=0;j<(Q*Q); ++j) {
+ D_store(sampcount, j) = D(j);
+ }
+ sigma_store(sampcount) = sigma2;
+ ++sampcount;
+ }
+ /////////////////////////////////////////////////
+ // REPORT
+ /////////////////////////////////////////////////
+ if(verbose > 0 && iter % verbose == 0){
+ Rprintf("\n ----------------------------------------------------------------------- \n");
+ Rprintf("\n Gaussian Panel iteration %i of %i ", (iter+1), tot_iter);
+ // store global estimates
+ Rprintf("\n beta = ");
+ for(int j=0;j<K; ++j) {
+ Rprintf("%10.5f", beta(j));
+ }
+ Rprintf("\n D[,] = ");
+ for(int j=0;j<(Q*Q); ++j) {
+ Rprintf("%10.5f", D(j));
+ }
+ Rprintf("\n sigma2 = %10.5f\n", sigma2);
+ }
+ }// END of MCMC loop
+
+ //////////////////////////////////////
+ if(chib == 1){
+ //////////////////////////////////////
+ Matrix<> beta_st(K, 1);
+ beta_st(_, 0) = meanc(beta_store);
+ const double sigma2_st = mean(sigma_store);
+ Matrix<> Dst = meanc(D_store);
+ Matrix<> D_st(Q, Q);
+ for (int k = 0; k<(Q*Q); ++k){
+ D_st[k] = Dst[k];
+ }
+ const Matrix <> Dinv_st = invpd(D_st);
+ Matrix<> beta_density(nstore, 1);
+
+ //////////////////////////////////////////////////////////////////
+ // 1. pdf.beta | D_g, sigma2_g
+ //////////////////////////////////////////////////////////////////
+ for (int iter=0; iter < nstore; ++iter){
+ for(int j=0;j<(Q*Q); ++j) {
+ D(j) = D_store(iter,j);
+ }
+ Dinv = invpd(D);
+ Matrix<double> XVX(K, K);
+ Matrix<double> XVY(K, 1);
+ for(int s = 0; s<nsubj; ++s) {
+ XVX += (1/sigma_store(iter))*cpXk_arr[s]-pow(1/sigma_store(iter),2)*tXWk_arr[s]*invpd(Dinv + cpWk_arr[s]*(1/sigma_store(iter)))*tWXk_arr[s];
+ XVY += (1/sigma_store(iter))*tXYk_arr[s]-pow(1/sigma_store(iter),2)*tXWk_arr[s]*invpd(Dinv + cpWk_arr[s]*(1/sigma_store(iter)))*tWYk_arr[s];
+ }
+ Matrix<> beta_var = invpd(B0 + XVX/sigma_store(iter));
+ Matrix<> beta_mean = beta_var*(B0*b0 + XVY/sigma_store(iter));
+ if (K == 1){
+ beta_density(iter) = dnorm(beta_st(0), beta_mean[0], sqrt(beta_var[0]));
+ }
+ else{
+ beta_density(iter) = ::exp(lndmvn(beta_st, beta_mean, beta_var));
+ }
+ }
+ double pdf_beta = log(mean(beta_density));
+
+ //////////////////////////////////////////////////////////////////
+ // 2. pdf.D
+ //////////////////////////////////////////////////////////////////
+ Matrix <> D_density(nstore, 1);
+ for (int iter=0; iter < nstore; ++iter){
+
+ /////////////////////////////////////////////////
+ // Step 2.1 Sample bi (random-effects coef)
+ /////////////////////////////////////////////////
+ Matrix<double> bi(Q, nsubj);
+ for(int s = 0; s<nsubj; ++ s) {
+ Matrix<double> b_var = invpd(Dinv + cpWk_arr[s]/sigma2);
+ Matrix<double> b_mean = b_var*t(Wk_arr[s])*(Yk_arr[s] - Xk_arr[s]*beta_st)/sigma2;
+ bi(_,s) = stream.rmvnorm(b_mean, b_var);
+ }
+
+ /////////////////////////////////////////////////
+ // Step 2.2 Sample sigma2
+ /////////////////////////////////////////////////
+ double SSE = 0;
+ for(int s=0;s<nsubj; ++s){
+ int ntime_s = subject_groupinfo(s, 2);
+ Matrix<> e = t(Yk_arr[s]-Xk_arr[s]*beta_st - Wk_arr[s]*bi(_,s))*
+ (Yk_arr[s] - Xk_arr[s]*beta_st - Wk_arr[s]*bi(_,s));
+ SSE = SSE + e[0];
+ }
+ double nu = (c0 + NOBS)/2;
+ double delta = (d0 + SSE)/2;
+ sigma2 = stream.rigamma(nu, delta);
+
+ /////////////////////////////////////////////////
+ // Step 2.3 Sample D
+ /////////////////////////////////////////////////
+ Matrix<double> SSB = bi*t(bi);
+ Matrix<double> D_scale = invpd(R0inv + SSB);
+ int D_dof = r0 + nsubj;
+ Dinv = stream.rwish(D_dof, D_scale);
+ D_density(iter) = ::exp(lndwish(invpd(D_st), D_dof, D_scale));
+ D = inv(Dinv);
+ }
+ double pdf_D = log(mean(D_density));
+
+ //////////////////////////////////////////////////////////////////
+ // 3. pdf.sigma2
+ //////////////////////////////////////////////////////////////////
+ Matrix <> sigma_density(nstore, 1);
+ for (int iter=0; iter < nstore; ++iter){
+
+ /////////////////////////////////////////////////
+ // Step 3.1 Sample bi (random-effects coef)
+ /////////////////////////////////////////////////
+ Matrix<double> bi(Q, nsubj);
+ for(int s = 0; s<nsubj; ++ s) {
+ Matrix<double> b_var = invpd(Dinv_st + cpWk_arr[s]/sigma2);
+ Matrix<double> b_mean = b_var*tWk_arr[s]*(Yk_arr[s] - Xk_arr[s]*beta_st)/sigma2;
+ bi(_,s) = stream.rmvnorm(b_mean, b_var);
+ }
+
+ /////////////////////////////////////////////////
+ // Step 3.2 Sample sigma2
+ /////////////////////////////////////////////////
+ double SSE = 0;
+ for(int s=0;s<nsubj; ++s){
+ int ntime_s = subject_groupinfo(s, 2);
+ Matrix<> e = t(Yk_arr[s]-Xk_arr[s]*beta_st - Wk_arr[s]*bi(_,s))*
+ (Yk_arr[s] - Xk_arr[s]*beta_st - Wk_arr[s]*bi(_,s));
+ SSE = SSE + e[0];
+ }
+ double nu = (c0 + NOBS)/2;
+ double delta = (d0 + SSE)/2;
+ sigma2 = stream.rigamma(nu, delta);
+ sigma_density(iter) = ::exp(lndinvgamma(sigma2_st, nu, delta));
+ }
+ double pdf_sigma2 = log(mean(sigma_density));
+ //////////////////////////////////////////////////////////////////
+ // likelihood f(y|beta_st, D_st, sigma2_st)
+ //////////////////////////////////////////////////////////////////
+ loglike = 0;
+ for(int s = 0; s<nsubj; ++s) {
+ int ntime_s = subject_groupinfo(s, 2);
+ Matrix<> Sigma = sigma2_st*eye(ntime_s) + Wk_arr[s] *D_st* t(Wk_arr[s]);
+ Matrix<> Mu = Xk_arr[s]*beta_st;
+ loglike += lndmvn(Yk_arr[s], Mu, Sigma);
+ }
+
+ //////////////////////
+ // log prior ordinate
+ //////////////////////
+ double density_beta_prior = 0;
+ 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);
+ }
+ double density_D_prior = lndwish(Dinv_st, r0, R0);
+ double density_sigma2_prior = lndinvgamma(sigma2_st, c0/2, d0/2);
+
+ // compute marginal likelihood
+ double logprior = density_beta_prior + density_sigma2_prior + density_D_prior;
+ logmarglike = (loglike + logprior) - (pdf_beta + pdf_sigma2 + pdf_D);
+
+ if (verbose >0){
+ Rprintf("\n ----------------------------------------------------------------------- \n");
+ Rprintf("\n logmarglike %10.5f", logmarglike);
+ Rprintf("\n loglike %10.5f", loglike, "\n");
+ Rprintf("\n log_prior %10.5f", logprior, "\n");
+ Rprintf("\n pdf_beta is %10.5f", pdf_beta, "\n");
+ Rprintf("\n pdf_D is %10.5f", pdf_D, "\n");
+ Rprintf("\n pdf_sigma2 is %10.5f", pdf_sigma2, "\n");
+ }
+ }// end of marginal likelihood computation
+
+ delete [] Xk_arr;
+ delete [] Yk_arr;
+ delete [] Wk_arr;
+ delete [] cpWk_arr;
+ delete [] cpXk_arr;
+ delete [] tWk_arr;
+ delete [] tXk_arr;
+ delete [] tXYk_arr;
+ delete [] tXWk_arr;
+ delete [] tWXk_arr;
+ delete [] tWYk_arr;
+
+}
+
+
+template <typename RNGTYPE>
+void HMMpanelRE_impl(rng<RNGTYPE>& stream,
+ unsigned int nsubj, unsigned int ntime,
+ unsigned int m, unsigned int nobs,
+ const Matrix<int>& subjectid, const Matrix<int>& timeid,
+ const Matrix<>& Y, const Matrix<>& X, const Matrix<>& W,
+ const Matrix<>& YT, const Matrix<>& XT, const Matrix<>& WT,
+ unsigned int burnin, unsigned int mcmc,
+ unsigned int thin, unsigned int verbose,
+ Matrix<>& betastart, double sigma2start, Matrix<>& Dstart,
+ const Matrix<>& b0, const Matrix<>& B0,
+ const double c0, const double d0,
+ unsigned int r0, const Matrix<>& R0, const Matrix<>& P0,
+ const Matrix<>& time_groupinfo,
+ const Matrix<>& subject_groupinfo,
+ Matrix<>& beta_store, Matrix<>& sigma_store,
+ Matrix<>& D_store,
+ Matrix<>& P_store, Matrix<>& ps_store, Matrix<>& s_store,
+ double& logmarglike, double& loglike,
+ unsigned int chib)
+{ // redefine constants
+ const unsigned int K = X.cols();
+ const unsigned int Q = W.cols();
+ const int NOBS = nobs;
+
+ const Matrix<> R0inv = invpd(R0);
+ const Matrix<> B0inv = invpd(B0);
+
+ const int tot_iter = burnin + mcmc;
+ const int nstore = mcmc / thin; // number of draws to store
+ const int ns = m + 1;
+
+ Matrix<>* D = new Matrix<>[ns];
+ Matrix<>* Dinv = new Matrix<>[ns];
+ for (int j=0; j<ns; ++j){
+ D[j] = Dstart;
+ Dinv[j] = invpd(Dstart);
+ }
+
+ Matrix<> D_record(ns, Q*Q);// for record D in D_store
+ for (int j = 0; j<ns; ++j){
+ for (int i=0; i<(Q*Q); ++i){
+ D_record(j,i) = D[j](i);
+ }
+ }
+
+
+ // generate posk_arr and post_arr
+ // Number of observations by group k
+ int *nobsk = new int[nsubj];
+ for (int k=0; k<nsubj; k++) {
+ nobsk[k]=0;
+ for (int n=0; n<NOBS; n++) {
+ if (subjectid[n]==k+1) {
+ nobsk[k]+=1;
+ }
+ }
+ }
+
+ // Position of each group in the data-set
+ int **posk_arr = new int*[nsubj];
+ for (int k=0; k<nsubj; k++) {
+ posk_arr[k] = new int[nobsk[k]];
+ int repk=0;
+ for (int n=0; n<NOBS; n++) {
+ if (subjectid[n]==k+1) {
+ posk_arr[k][repk]=n;
+ repk++;
+ }
+ }
+ }
+
+ // Small fixed matrices indexed on k for data access
+ Matrix<double> *Yk_arr = new Matrix<double>[nsubj];
+ Matrix<double> *Xk_arr = new Matrix<double>[nsubj];
+ Matrix<double> *Wk_arr = new Matrix<double>[nsubj];
+ for(int k=0; k<nsubj; k++) {
+ Xk_arr[k] = Matrix<double>(nobsk[k], K);
+ Wk_arr[k] = Matrix<double>(nobsk[k], Q);
+ Yk_arr[k] = Matrix<double>(nobsk[k], 1);
+ for (int l=0; l<nobsk[k]; l++) {
+ for (int p=0; p< K; p++) {
+ Xk_arr[k](l, p) = X[p*NOBS + posk_arr[k][l]];
+ }
+ for (int q=0; q < Q; q++) {
+ Wk_arr[k](l, q) = W[q*NOBS + posk_arr[k][l]];
+ }
+ Yk_arr[k](l,0) = Y[posk_arr[k][l]];
+ }
+ }
+ // number of observations by time t
+ int *nobst = new int[ntime];
+ for (int k=0; k<ntime; k++) {
+ nobst[k]=0;
+ for (int n=0; n<NOBS; n++) {
+ if (timeid[n]==k + 1) {
+ nobst[k]+=1;
+ }
+ }
+ }
+ // Position of each group in the data-set
+ int **post_arr = new int*[ntime];
+ for (int k=0; k<ntime; k++) {
+ post_arr[k] = new int[nobst[k]];
+ int repk=0;
+ for (int n=0; n<NOBS; n++) {
+ if (timeid[n]==k+1) {
+ post_arr[k][repk]=n;
+ repk++;
+ }
+ }
+ }
+
+ // XTarr is data transformed for multivariate TS analsysi
+ Matrix<>* Xt_arr = new Matrix<>[ntime];
+ Matrix<>* Wt_arr = new Matrix<>[ntime];
+ Matrix<>* Yt_arr = new Matrix<>[ntime];
+ for(int k=0; k<ntime; k++) {
+ if (nobst[k] > 0){
+ Xt_arr[k] = Matrix<double>(nobst[k], K);
+ Wt_arr[k] = Matrix<double>(nobst[k], Q);
+ Yt_arr[k] = Matrix<double>(nobst[k], 1);
+ for (int l = 0; l<nobst[k]; l++) {
+ for (int p = 0; p < K; p++) {
+ Xt_arr[k](l, p) = X[p*NOBS+post_arr[k][l]];
+ }
+ for (int q = 0; q < Q; q++) {
+ Wt_arr[k](l, q) = W[q*NOBS+post_arr[k][l]];
+ }
+ Yt_arr[k](l,0) = Y[post_arr[k][l]];
+ }
+ }
+ }
+
+ Matrix<double> *cpWt_arr = new Matrix<double>[ntime];
+ Matrix<double> *cpXt_arr = new Matrix<double>[ntime];
+ Matrix<double> *tXWt_arr = new Matrix<double>[ntime];
+ Matrix<double> *tWXt_arr = new Matrix<double>[ntime];
+ Matrix<double> *tWYt_arr = new Matrix<double>[ntime];
+ Matrix<double> *tXYt_arr = new Matrix<double>[ntime];
+ Matrix<double> *tXt_arr = new Matrix<double>[ntime];
+ Matrix<double> *tWt_arr = new Matrix<double>[ntime];
+ for(int k=0; k<ntime; k++) {
+ cpXt_arr[k] = crossprod(Xt_arr[k]);
+ cpWt_arr[k] = crossprod(Wt_arr[k]);
+ tXt_arr[k] = t(Xt_arr[k]);
+ tWt_arr[k] = t(Wt_arr[k]);
+ tXWt_arr[k] = t(Xt_arr[k])*Wt_arr[k];
+ tWXt_arr[k] = t(Wt_arr[k])*Xt_arr[k];
+ tWYt_arr[k] = t(Wt_arr[k])*Yt_arr[k];
+ tXYt_arr[k] = t(Xt_arr[k])*Yt_arr[k];
+ }
+
+ // starting values
+ Matrix<double> *bk_run = new Matrix<double>[nsubj]; // Random effects
+ for (int k=0;k<nsubj;k++) {
+ bk_run[k] = Matrix<double>(Q,1);
+ }
+ Matrix<>* bi = new Matrix<>[ns];
+ Matrix<> beta(ns, K);
+ Matrix<> sigma2(ns, 1);
+ for (int j=0; j<ns; ++j){
+ bi[j] = Matrix<>(Q, nsubj);
+ beta(j,_) = betastart;
+ sigma2(j) = sigma2start;
+ }
+ Matrix<> P = P0;
+
+ // MCMC iterations start here
+ int sampcount = 0;
+ // initialize Yhat for the first loop only
+ for (int iter=0; iter < tot_iter; ++iter){
+ //////////////////////
+ // Step 1. Sample state
+ //////////////////////
+ Matrix<> F(ntime, ns);
+ Matrix<> pr1(ns, 1);
+ pr1[0] = 1;
+ Matrix<> py(ns, 1);
+ Matrix<> pstyt1(ns, 1);
+ for (int tt=0; tt<ntime ; ++tt) {
+ int nsubj_s = time_groupinfo(tt, 2);
+ // Matrix<> Wbi(nsubj_s, 1);
+ for (int j=0; j<ns;++j){
+ Matrix<> Sigma = eye(nsubj_s)*sigma2[j];
+ Matrix<> WDW = Wt_arr[tt]*D[j]*tWt_arr[tt];
+ Matrix<> Mu = Xt_arr[tt]*::t(beta(j,_));
+ py[j] = ::exp(lndmvn(Yt_arr[tt], Mu, WDW + Sigma));
+ }
+ if (tt==0) pstyt1 = pr1;
+ else {
+ pstyt1 = ::t(F(tt-1,_)*P);
+ }
+ Matrix<> unnorm_pstyt = pstyt1%py;
+ Matrix<> pstyt = unnorm_pstyt/sum(unnorm_pstyt);
+ for (int j=0; j<ns ; ++j) {
+ F(tt,j) = pstyt(j);
+ }
+ }
+ Matrix<int> state(ntime, 1);
+ Matrix<> ps = Matrix<>(ntime, ns);
+ ps(ntime-1,_) = F(ntime-1,_);
+ state(ntime-1) = ns;
+
+ Matrix<> pstyn = Matrix<>(ns, 1);
+ double pone = 0.0;
+ int tt = ntime-2;
+
+ while (tt >= 0){
+ int st = state(tt+1);
+ Matrix<> Pst_1 = ::t(P(_, st-1));
+ Matrix<> unnorm_pstyn = F(tt,_)%Pst_1;
+ pstyn = unnorm_pstyn/sum(unnorm_pstyn);
+
+ if (st==1){
+ state(tt) = 1;
+ }
+ else{
+ pone = pstyn(st-2);
+ if(stream.runif() < pone) state(tt) = st-1;
+ else state(tt) = st;
+ }
+ ps(tt,_) = pstyn;
+ --tt;
+ }// end of while loop
+
+ //////////////////////
+ // Step 2. Sample beta
+ //////////////////////
+ Matrix<int> nstate(ns, 1); // refresh the container for numbers of each state
+ int beta_count = 0;
+ for (int j = 0; j<ns; ++j){
+ for (int i = 0; i<ntime; ++i){
+ if (state[i] == j + 1) {
+ nstate[j] = nstate[j] + 1;
+ }
+ }
+ beta_count = beta_count + nstate[j];
+ Matrix<> XVX(K, K);
+ Matrix<> XVY(K, 1);
+ for(int tt = (beta_count - nstate[j]); tt<beta_count; ++tt) {
+ int nsubj_s = time_groupinfo(tt, 2);
+ XVX += (1/sigma2[j])*cpXt_arr[tt]-pow(1/sigma2[j],2)*tXWt_arr[tt]*invpd(Dinv[j] + cpWt_arr[tt]*(1/sigma2[j]))*tWXt_arr[tt];
+ XVY += (1/sigma2[j])*tXYt_arr[tt]-pow(1/sigma2[j],2)*tXWt_arr[tt]*invpd(Dinv[j] + cpWt_arr[tt]*(1/sigma2[j]))*tWYt_arr[tt];
+ }
+ Matrix<> beta_var = invpd(B0 + XVX/sigma2[j]);
+ Matrix<> beta_mean = beta_var*(B0*b0 + XVY/sigma2[j]);
+ beta(j,_) = stream.rmvnorm(beta_mean, beta_var);
+ }
+
+ /////////////////////////////////////////////////
+ // Step 3. Sample bi (random-effects coef)
+ /////////////////////////////////////////////////
+ beta_count = 0;
+ Matrix<int> YN(ns, 1);
+ Matrix<> SSE(ns, 1);
+ for (int j = 0; j<ns; ++j){
+ beta_count = beta_count + nstate[j];
+ for(int s = 0; s<nsubj; ++ s) {
+ Matrix<> yj = Yk_arr[s]((beta_count - nstate[j]), 0, (beta_count - 1), 0);
+ Matrix<> Xj = Xk_arr[s]((beta_count - nstate[j]), 0, (beta_count - 1), K-1);
+ Matrix<> Wj = Wk_arr[s]((beta_count - nstate[j]), 0, (beta_count - 1), Q-1);
+ Matrix<> b_var = invpd(Dinv[j]+ ::t(Wj)*Wj/sigma2[j]);
+ Matrix<> b_mean = b_var*::t(Wj)*(yj - Xj*::t(beta(j,_)))/sigma2[j];
+ bi[j](_,s) = stream.rmvnorm(b_mean, b_var);
+
+ // FOR SIGMA
+ YN[j] = YN[j] + yj.rows();
+ Matrix<> e = ::t(yj - Xj*::t(beta(j,_)) - Wj*bi[j](_,s))*(yj - Xj*::t(beta(j,_)) - Wj*bi[j](_,s));
+ SSE[j] = SSE[j] + e[0];
+ }
+ }
+
+ /////////////////////////////////////////////////
+ // Step 4. Sample sigma2
+ /////////////////////////////////////////////////
+ for (int j = 0; j<ns; ++j){
+ double nu = (c0 + (double)YN[j])*0.5;
+ double scale = (d0 + SSE[j])*0.5;
+ sigma2[j] = stream.rigamma(nu, scale);
+ }
+
+ /////////////////////////////////////////////////
+ // Step 5. Sample D
+ /////////////////////////////////////////////////
+ for (int j = 0; j<ns; ++j){
+ Matrix<> SSB = bi[j]*t(bi[j]);
+ Matrix<> D_scale = invpd(R0inv + SSB);
+ int D_dof = r0 + nsubj;
+ Dinv[j] = stream.rwish(D_dof, D_scale);
+ D[j] = invpd(Dinv[j]);
+
+ for (int i=0; i<(Q*Q); ++i){
+ D_record(j,i) = D[j](i);
+ }
+ }
+
+ //////////////////////
+ // Step 6. Sample P
+ //////////////////////
+ double shape1 = 0;
+ double shape2 = 0;
+ for (int j =0; j<m; ++j){
+ shape1 = std::abs(P0(j,j) + nstate[j] - 1);
+ shape2 = P0(j,j+1) + 1; //
+ P(j,j) = stream.rbeta(shape1, shape2);
+ P(j,j+1) = 1 - P(j,j);
+ }
+ P(m, m) = 1;
+
+ /////////////////////////////////////////////////
+ // STORE
+ /////////////////////////////////////////////////
+ if (iter >= burnin && ((iter % thin) == 0)) {
+ Matrix<> tbeta = ::t(beta); //transpose beta for R output
+ for (int i=0; i<(ns*K); ++i){
+ beta_store(sampcount, i) = tbeta[i];// stored by the order of (11, 12, 13, 21, 22, 23)
+ }
+ for (int i=0; i<ns; ++i){
+ sigma_store(sampcount, i) = sigma2[i];
+ }
+ Matrix<> tD_record = ::t(D_record);
+ for (int i=0; i<(ns*Q*Q); ++i){
+ D_store(sampcount,i) = tD_record(i);// stored by the order of (D11, D12, D13, 21, 22, 23)
+ }
+ for (int j=0; j<ns*ns; ++j){
+ P_store(sampcount,j)= P[j];
+ }
+ s_store(sampcount,_) = state;
+ for (int l=0; l<ntime ; ++l){
+ ps_store(l,_) = ps_store(l,_) + ps(l,_);
+ }
+ ++sampcount;
+ }
+
+
+ /////////////////////////////////////////////////
+ // REPORT
+ /////////////////////////////////////////////////
+ if(verbose > 0 && iter % verbose == 0){
+ Rprintf("\n ----------------------------------------------------------------------- \n");
+ Rprintf("HMMpanelRE iteration %i of %i \n", (iter+1), tot_iter);
+ for(int j=0;j<ns; ++j) {
+ Rprintf("\n number of obs from state %i is %i", j, nstate[j]);
+ }
+ for(int j=0;j<ns; ++j) {
+ Rprintf("\n beta for state %i = ", j);
+ for(int i=0;i<K; ++i) {
+ Rprintf("%10.5f", beta(j, i));
+ }
+ }
+ for(int j=0;j<ns; ++j) {
+ Rprintf("\n sigma2 for state %i is %10.5f", j, sigma2[j]);
+ }
+ }
+
+ }// END of MCMC loop
+
+ //////////////////////////////////////
+ if(chib == 1){
+ //////////////////////////////////////
+
+ // posterior ordinate
+ 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> sigma2st = meanc(sigma_store);
+ Matrix<double, Row> sigma2_st(ns, 1);
+ for (int j = 0; j< ns; ++j){
+ sigma2_st[j] = sigma2st[j];
+ }
+
+ Matrix<double> Dst = meanc(D_store);
+ Matrix<>* Dinv_st = new Matrix<>[ns];
+ Matrix<>* D_st = new Matrix<>[ns];
+ int Dcount = 0;
+ for (int j=0; j<ns; ++j){
+ Matrix<> Dtemp(Q, Q);
+ for (int k = 0; k<(Q*Q); ++k){
+ Dtemp[k] = Dst[Dcount + k];
+ }
+ Dcount = Dcount + Q*Q;
+ D_st[j] = Dtemp;
+ Dinv_st[j] = invpd(Dtemp);
+ }
+ Matrix<> Dst_record(ns, Q*Q);// for record D in D_store
+ for (int j = 0; j<ns; ++j){
+ for (int i=0; i<(Q*Q); ++i){
+ Dst_record(j,i) = D_st[j](i);
+ }
+ }
+
+ 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];
+ }
+
+ Matrix<> density_beta(nstore, ns);
+ Matrix<> density_local_beta(ns, 1);
+ //////////////////////////////////////////////////////////////////
+ // 1. pdf.beta | D_g, sigma2_g, P_g, bi_g
+ //////////////////////////////////////////////////////////////////
+ Matrix<> Dmcmc(Q, Q);
+
+ for (int iter = 0; iter<nstore; ++iter){
+ int Dcount = 0;
+ int beta_count = 0;
+ Matrix<int> nstate(ns, 1);
+ for (int j = 0; j<ns; ++j){
+ for (int i = 0; i<ntime; ++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];
+
+ for (int h = 0; h< (Q*Q); ++h){
+ Dmcmc[h] = D_store(iter, Dcount + h);
+ }
+ Dcount = Dcount + Q*Q;
+ Matrix<> Dinv_j = invpd(Dmcmc);
+ Matrix<> XVX(K, K);
+ Matrix<> XVY(K, 1);
+ for(int tt = (beta_count - nstate[j]); tt<beta_count; ++tt) {
+ int nsubj_s = time_groupinfo(tt, 2);
+ XVX += (1/sigma_store(iter, j))*cpXt_arr[tt]-pow(1/sigma_store(iter, j),2)*tXWt_arr[tt]*invpd(Dinv_j + cpWt_arr[tt]*(1/sigma_store(iter, j)))*tWXt_arr[tt];
+ XVY += (1/sigma_store(iter, j))*tXYt_arr[tt]-pow(1/sigma_store(iter, j),2)*tXWt_arr[tt]*invpd(Dinv_j + cpWt_arr[tt]*(1/sigma_store(iter, j)))*tWYt_arr[tt];
+ }
+ Matrix<> beta_var = invpd(B0 + XVX/sigma_store(iter, j));
+ Matrix<> beta_mean = beta_var*(B0*b0 + XVY/sigma_store(iter, j));
+ if (K == 1){
+ density_beta(iter, j) = dnorm(beta_st(j), beta_mean[0], sqrt(beta_var[0]));
+ }
+ else{
+ density_beta(iter, j) = ::exp(lndmvn(::t(beta_st(j,_)), beta_mean, beta_var));
+ }
+ }
+ }
+ double pdf_beta = log(prod(meanc(density_beta)));
+
+ //////////////////////////////////////////////////////////////////
+ // 2. pdf.D
+ //////////////////////////////////////////////////////////////////
+ Matrix<> density_D(nstore, ns);
+ for (int iter = 0; iter<nstore; ++iter){
+ Matrix<> SSE(ns, 1);
+ int beta_count = 0;
+ Matrix<int> YN(ns, 1);
+ Matrix<int> nstate(ns, 1);
+ /////////////////////////////////////////////////
+ // 2.1. s | beta_st, y, bi, sigma, D, P
+ /////////////////////////////////////////////////
+ Matrix<> F(ntime, ns);
+ Matrix<> pr1(ns, 1);
+ pr1[0] = 1;
+ Matrix<> py(ns, 1);
+ Matrix<> pstyt1(ns, 1);
+ for (int tt=0; tt<ntime ; ++tt) {
+ int nsubj_s = time_groupinfo(tt, 2);
+ Matrix<> Wbi(nsubj_s, 1);
+ for (int j=0; j<ns;++j){
+ Matrix<> Sigma = eye(nsubj_s)*sigma2[j];
+ Matrix<> WDW = Wt_arr[tt]*D[j]*tWt_arr[tt];
+ Matrix<> Mu = Xt_arr[tt]*::t(beta_st(j,_));
+ py[j] = ::exp(lndmvn(Yt_arr[tt], Mu, WDW+Sigma));
+ }
+ if (tt==0) pstyt1 = pr1;
+ else {
+ pstyt1 = ::t(F(tt-1,_)*P);
+ }
+ Matrix<> unnorm_pstyt = pstyt1%py;
+ Matrix<> pstyt = unnorm_pstyt/sum(unnorm_pstyt);
+ for (int j=0; j<ns ; ++j) {
+ F(tt,j) = pstyt(j);
+ }
+ }// end of F matrix filtering
+ Matrix<int> state(ntime, 1);
+ Matrix<> ps = Matrix<>(ntime, ns);
+ ps(ntime-1, _) = F(ntime-1, _);
+ state(ntime-1) = ns;
+
+ Matrix<> pstyn = Matrix<>(ns, 1);
+ double pone = 0.0;
+ int tt = ntime-2;
+
+ while (tt >= 0){
+ int st = state(tt+1);
+ Matrix<> Pst_1 = ::t(P(_, st-1));
+ Matrix<> unnorm_pstyn = F(tt,_)%Pst_1;
+ pstyn = unnorm_pstyn/sum(unnorm_pstyn);
+
+ if (st==1){
+ state(tt) = 1;
+ }
+ else{
+ pone = pstyn(st-2);
+ if(stream.runif() < pone) state(tt) = st-1;
+ else state(tt) = st;
+ }
+ ps(tt,_) = pstyn;
+ --tt;
+
+ }// end of while loop
+
+ /////////////////////////////////////////////////
+ // 2.2. P | beta_st, y, bi, sigma, s, D
+ /////////////////////////////////////////////////
+ for (int j = 0; j<ns; ++j){
+ for (int i = 0; i<ntime; ++i){
+ if (state(i) == j + 1) { // since j starts from 0
+ nstate[j] = nstate[j] + 1;
+ }// end of if
+ }// end of int i<n
+ }
+ double shape1 = 0;
+ double shape2 = 0;
+ for (int j =0; j<m; ++j){
+ shape1 = std::abs(P0(j,j) + nstate[j] - 1);
+ shape2 = P0(j,j+1) + 1; //
+ P(j,j) = stream.rbeta(shape1, shape2);
+ P(j,j+1) = 1 - P(j,j);
+ }
+ P(m, m) = 1;
+ /////////////////////////////////////////////////
+ // 2.3. bi | beta_st, y, D, sigma, s, P
+ /////////////////////////////////////////////////
+ beta_count = 0;
+ for (int j = 0; j<ns; ++j){
+ beta_count = beta_count + nstate[j];
+ for(int s = 0; s<nsubj; ++ s) {
+ Matrix<> yj = Yk_arr[s]((beta_count - nstate[j]), 0, (beta_count - 1), 0);
+ Matrix<> Xj = Xk_arr[s]((beta_count - nstate[j]), 0, (beta_count - 1), K-1);
+ Matrix<> Wj = Wk_arr[s]((beta_count - nstate[j]), 0, (beta_count - 1), Q-1);
+ Matrix<> b_var = invpd(Dinv[j]+ ::t(Wj)*Wj/sigma2[j]);
+ Matrix<> b_mean = b_var*::t(Wj)*(yj - Xj*::t(beta_st(j,_)))/sigma2[j];
+ bi[j](_,s) = stream.rmvnorm(b_mean, b_var);
+
+ // FOR SIGMA
+ YN[j] = YN[j] + yj.rows();
+ Matrix<> e = ::t(yj - Xj*::t(beta_st(j,_)) - Wj*bi[j](_,s))*(yj - Xj*::t(beta_st(j,_)) - Wj*bi[j](_,s));
+ SSE[j] = SSE[j] + e[0];
+ }
+ }
+ /////////////////////////////////////////////////
+ // 2.4. D | beta_st, y, bi, sigma, s, P
+ /////////////////////////////////////////////////
+ for (int j = 0; j<ns; ++j){
+ Matrix<> SSB = bi[j]*t(bi[j]);
+ Matrix<> D_scale = invpd(R0inv + SSB);
+ int D_dof = r0 + nsubj;
+ Dinv[j] = stream.rwish(D_dof, D_scale);
+ density_D(iter, j) = ::exp(lndwish(Dinv_st[j], D_dof, D_scale));
+ D[j] = inv(Dinv[j]);
+ }
+ /////////////////////////////////////////////////
+ // 2.5. sigma | beta_st, y, bi, D, s, P
+ /////////////////////////////////////////////////
+ for (int j = 0; j<ns; ++j){
+ double nu = (c0 + (double)YN[j])/2;
+ double scale = (d0 + SSE[j])/2;
+ sigma2[j] = stream.rigamma(nu, scale);
+ }
+ }// end of reduced run for step 2
+
+ double pdf_D = log(prod(meanc(density_D)));
+
+
+ //////////////////////////////////////////////////////////////////
+ // 3. pdf.sigma2 : sigma2| beta_st, D_st, y, bi, s, P
+ //////////////////////////////////////////////////////////////////
+ Matrix<> density_sigma2(nstore, ns);
+ for (int iter = 0; iter<nstore; ++iter){
+ Matrix<> SSE(ns, 1);
+ int beta_count = 0;
+ Matrix<int> YN(ns, 1);
+ Matrix<int> nstate(ns, 1);
+
+ /////////////////////////////////////////////////
+ // 3.1. s| beta_st, D_st, y, bi, sigma2, P
+ /////////////////////////////////////////////////
+ Matrix<> F(ntime, ns);
+ Matrix<> pr1(ns, 1);
+ pr1[0] = 1;
+ Matrix<> py(ns, 1);
+ Matrix<> pstyt1(ns, 1);
+
+ for (int tt=0; tt<ntime ; ++tt) {
+ int nsubj_s = time_groupinfo(tt, 2);
+ Matrix<> Wbi(nsubj_s, 1);
+ for (int j=0; j<ns;++j){
+ Matrix<> Sigma = eye(nsubj_s)*sigma2[j];
+ Matrix<> WDW = Wt_arr[tt]*D_st[j]*tWt_arr[tt];
+ Matrix<> Mu = Xt_arr[tt]*::t(beta_st(j,_));
+ py[j] = ::exp(lndmvn(Yt_arr[tt], Mu, WDW+Sigma));
+ }
+ if (tt==0) pstyt1 = pr1;
+ else {
+ pstyt1 = ::t(F(tt-1,_)*P);
+ }
+ Matrix<> unnorm_pstyt = pstyt1%py;
+ Matrix<> pstyt = unnorm_pstyt/sum(unnorm_pstyt);
+ for (int j=0; j<ns ; ++j) {
+ F(tt,j) = pstyt(j);
+ }
+
+ }// end of F matrix filtering
+ Matrix<int> state(ntime, 1);
+ Matrix<> ps = Matrix<>(ntime, ns);
+ ps(ntime-1,_) = F(ntime-1,_);
+ state(ntime-1) = ns;
+ Matrix<> pstyn = Matrix<>(ns, 1);
+ double pone = 0.0;
+ int tt = ntime-2;
+
+ while (tt >= 0){
+ int st = state(tt+1);
+ Matrix<> Pst_1 = ::t(P(_, st-1));
+ Matrix<> unnorm_pstyn = F(tt,_)%Pst_1;
+ pstyn = unnorm_pstyn/sum(unnorm_pstyn);
+
+ if (st==1){
+ state(tt) = 1;
+ }
+ else{
+ pone = pstyn(st-2);
+ if(stream.runif() < pone) state(tt) = st-1;
+ else state(tt) = st;
+ }
+ ps(tt,_) = pstyn;
+ --tt;
+
+ }// end of while loop
+
+ /////////////////////////////////////////////////
+ // 3.2. bi| beta_st, D_st, y, sigma2, s, P
+ /////////////////////////////////////////////////
+ for (int j = 0; j<ns; ++j){
+ for (int i = 0; i<ntime; ++i){
+ if (state(i) == j + 1) { // since j starts from 0
+ nstate[j] = nstate[j] + 1;
+ }// end of if
+ }// end of int i<n
+ }
+ for (int j = 0; j<ns; ++j){
+ beta_count = beta_count + nstate[j];
+ for(int s = 0; s<nsubj; ++ s) {
+ Matrix<> yj = Yk_arr[s]((beta_count - nstate[j]), 0, (beta_count - 1), 0);
+ Matrix<> Xj = Xk_arr[s]((beta_count - nstate[j]), 0, (beta_count - 1), K-1);
+ Matrix<> Wj = Wk_arr[s]((beta_count - nstate[j]), 0, (beta_count - 1), Q-1);
+ Matrix<> b_var = inv(Dinv_st[j]+ ::t(Wj)*Wj/sigma2[j]);
+ Matrix<> b_mean = b_var*::t(Wj)*(yj - Xj*::t(beta_st(j,_)))/sigma2[j];
+ bi[j](_,s) = stream.rmvnorm(b_mean, b_var);
+
+ // FOR SIGMA
+ YN[j] = YN[j] + yj.rows();
+ Matrix<> e = ::t(yj - Xj*::t(beta_st(j,_)) - Wj*bi[j](_,s))*(yj - Xj*::t(beta_st(j,_)) - Wj*bi[j](_,s));
+ SSE[j] = SSE[j] + e[0];
+ }
+
+ /////////////////////////////////////////////////
+ // 3.3. sigma2| beta_st, D_st, y, bi, s, P
+ /////////////////////////////////////////////////
+ double nu = (c0 + (double)YN[j])/2;
+ double scale = (d0 + SSE[j])/2;
+ sigma2[j] = stream.rigamma(nu, scale);
+ density_sigma2(iter, j) = ::exp(lndinvgamma(sigma2_st[j], nu, scale));
+ }
+
+ /////////////////////////////////////////////////
+ // 3.4. P| beta_st, D_st, y, bi, sigma2, s
+ /////////////////////////////////////////////////
+ double shape1 = 0;
+ double shape2 = 0;
+ for (int j =0; j<m; ++j){
+ shape1 = std::abs(P0(j,j) + nstate[j] - 1);
+ shape2 = P0(j,j+1) + 1; //
+ P(j,j) = stream.rbeta(shape1, shape2);
+ P(j,j+1) = 1 - P(j,j);
+ }
+ P(m, m) = 1;
+ }// end of reduced run
+ double pdf_sigma2 = log(prod(meanc(density_sigma2)));
+
+ //////////////////////
+ // 4. pdf.P:
+ //////////////////////
+ Matrix<> density_P(nstore, ns);
+ for (int iter = 0; iter < nstore; ++iter){
+ Matrix<int> nstate(ns, 1);
+ /////////////////////////////////////////////////
+ // 4.1. s| y, P, beta_st, sigma_st, D_st
+ /////////////////////////////////////////////////
+ Matrix<> F(ntime, ns);
+ Matrix<> pr1(ns, 1);
+ pr1[0] = 1;
+ Matrix<> py(ns, 1);
+ Matrix<> pstyt1(ns, 1);
+
+ for (int tt=0; tt<ntime ; ++tt) {
+ int nsubj_s = time_groupinfo(tt, 2);
+ Matrix<> Wbi(nsubj_s, 1);
+ for (int j=0; j<ns;++j){
+ Matrix<> Sigma = eye(nsubj_s)*sigma2_st[j];
+ Matrix<> WDW = Wt_arr[tt]*D_st[j]*t(Wt_arr[tt]);
+ Matrix<> Mu = Xt_arr[tt]*::t(beta_st(j,_));
+ py[j] = ::exp(lndmvn(Yt_arr[tt], Mu, WDW+Sigma));
+ }
+
+ if (tt==0) pstyt1 = pr1;
+ else {
+ pstyt1 = ::t(F(tt-1,_)*P);
+ }
+
+ Matrix<> unnorm_pstyt = pstyt1%py;
+ Matrix<> pstyt = unnorm_pstyt/sum(unnorm_pstyt);
+ for (int j=0; j<ns ; ++j) {
+ F(tt,j) = pstyt(j);
+ }
+ }
+ Matrix<int> state(ntime, 1);
+ Matrix<> ps = Matrix<>(ntime, ns);
+ ps(ntime-1,_) = F(ntime-1,_);
+ state(ntime-1) = ns;
+ Matrix<> pstyn = Matrix<>(ns, 1);
+ double pone = 0.0;
+ int tt = ntime-2;
+
+ while (tt >= 0){
+ int st = state(tt+1);
+ Matrix<> Pst_1 = ::t(P(_, st-1));
+ Matrix<> unnorm_pstyn = F(tt,_)%Pst_1;
+ pstyn = unnorm_pstyn/sum(unnorm_pstyn);
+
+ if (st==1){
+ state(tt) = 1;
+ }
+ else{
+ pone = pstyn(st-2);
+ if(stream.runif() < pone) state(tt) = st-1;
+ else state(tt) = st;
+ }
+ ps(tt,_) = pstyn;
+ --tt;
+ }// end of while loop
+
+ for (int j = 0; j<ns; ++j){
+ for (int i = 0; i<ntime; ++i){
+ if (state(i) == j + 1) {
+ nstate[j] = nstate[j] + 1;
+ }// end of if
+ }// end of int i<n
+ }
+ /////////////////////////////////////////////////
+ // 4.2. P| y, s, beta_st, sigma_st, D_st
+ /////////////////////////////////////////////////
+ double shape1 = 0;
+ double shape2 = 0;
+ for (int j =0; j<m; ++j){
+ shape1 = std::abs(P0(j,j) + nstate[j] - 1);
+ shape2 = P0(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);
+ }
+ P(m, m) = 1;
+ density_P(iter, ns-1) = 1; //without this, there will be an error
+ }// end of pdf.P MCMC run
+
+ double pdf_P = log(prod(meanc(density_P)));
+
+
+ //////////////////////////////////////////////////////////////////
+ // likelihood f(y|beta_st, D_st, sigma2_st, P_st)
+ //////////////////////////////////////////////////////////////////
+ Matrix<> F(ntime, ns);
+ Matrix<> pr1(ns, 1);
+ Matrix<> like(ntime, 1);
+ pr1[0] = 1;
+ Matrix<> py(ns, 1);
+ Matrix<> pstyt1(ns, 1);
+ for (int tt=0; tt<ntime ; ++tt) {
+ Matrix<> Wbi(nsubj, 1);
+ for (int j=0; j<ns;++j){
+ Matrix<> Sigma = eye(nsubj)*sigma2_st[j];
+ Matrix<> WDW = Wt_arr[tt]*D_st[j]*t(Wt_arr[tt]);
+ Matrix<> Mu = Xt_arr[tt]*::t(beta_st(j,_));
+ py[j] = ::exp(lndmvn(Yt_arr[tt], Mu, WDW+Sigma));
+ }
+
+ if (tt==0) pstyt1 = pr1;
+ else {
+ pstyt1 = ::t(F(tt-1,_)*P_st);
+ }
+ Matrix<> unnorm_pstyt = pstyt1%py;
+ Matrix<> pstyt = unnorm_pstyt/sum(unnorm_pstyt);
+ for (int j=0; j<ns ; ++j) {
+ F(tt,j) = pstyt(j);
+ }
+ like[tt] = sum(unnorm_pstyt);
+ }
+
+ loglike = sum(log(like));
+
+ //////////////////////
+ // log prior ordinate
+ //////////////////////
+ Matrix<double> density_beta_prior(ns, 1);
+ Matrix<double> density_sigma2_prior(ns, 1);
+ Matrix<double> density_D_prior(ns, 1);
+ Matrix<double> density_P_prior(ns, 1);
+ density_P_prior[ns-1] = 0; //
+
+ for (int j=0; j<ns ; ++j){
+ if (K == 1){
+ density_beta_prior[j] = log(dnorm(beta_st(j), b0[0], sqrt(B0inv[0])));
+ }
+ else{
+ density_beta_prior[j] = lndmvn(::t(beta_st(j,_)), b0, B0inv);
+ }
+ density_D_prior[j] = lndwish(Dinv_st[j], r0, R0);
+ density_sigma2_prior[j] = lndinvgamma(sigma2_st[j], c0/2, d0/2);
+ }
+ for (int j=0; j<m ; ++j){
+ density_P_prior[j] = log(dbeta(P_st(j,j), P0(j,j), P0(j,j+1)));
+ }
+
+ // compute marginal likelihood
+ double logprior = sum(density_beta_prior) + sum(density_sigma2_prior) + sum(density_D_prior) +
+ sum(density_P_prior);
+ logmarglike = (loglike + logprior) - (pdf_beta + pdf_P + pdf_sigma2 + pdf_D);
+
+ if (verbose >0){
+ Rprintf("\n ----------------------------------------------------------------------- \n");
+ Rprintf("\nlogmarglike = %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_sigma2);
+ Rprintf("pdf_D = %10.5f\n", pdf_D);
+ Rprintf("pdf_P = %10.5f\n", pdf_P);
+ }
+
+ }// end of marginal likelihood computation
+
+ delete[] nobst;
+ for(int k=0; k<ntime; k++) {
+ delete[] post_arr[k];
+ }
+ delete[] nobsk;
+ for(int k=0; k<nsubj; k++) {
+ delete[] posk_arr[k];
+ }
+ delete [] posk_arr;
+ delete [] Xk_arr;
+ delete [] Yk_arr;
+ delete [] Wk_arr;
+ delete [] Xt_arr;
+ delete [] tXt_arr;
+ delete [] Yt_arr;
+ delete [] cpWt_arr;
+ delete [] cpXt_arr;
+ delete [] Wt_arr;
+ delete [] tWt_arr;
+ delete [] tXYt_arr;
+ delete [] tXWt_arr;
+ delete [] tWXt_arr;
+ delete [] tWYt_arr;
+ delete [] bi;
+}// end of HMMpanelRE_impl
+
+extern "C" {
+ void HMMpanelRE(double* betadata, const int* betarow, const int* betacol,
+ double* sigmadata, double* Ddata, double *psout, double *sout,
+ const int* nsubj, const int* ntime,
+ const int* m, const int* nobs, const int* subjectid, const int* timeid,
+ const double* Ydata, const int* Yrow, const int* Ycol,
+ const double* Xdata, const int* Xrow, const int* Xcol,
+ const double* Wdata, const int* Wrow, const int* Wcol,
+ const double* YTdata, const double* XTdata, const double* WTdata,
+ const int* burnin, const int* mcmc, const int* thin, const int* verbose,
+ const int *uselecuyer, const int *seedarray, const int *lecuyerstream,
+ const double* betastartdata, const double* sigma2start,
+ const double *Pstart,
+ const double* b0data, const double* B0data,
+ const double* c0, const double* d0, const int* r0, const double* R0data,
+ const double* subject_groupinfodata, const double* time_groupinfodata,
+ double *logmarglikeholder, double *loglikeholder,
+ const int *chib){
+ // pull together Matrix objects
+ Matrix<> Y(*Yrow, *Ycol, Ydata);
+ Matrix<> X(*Xrow, *Xcol, Xdata);
+ Matrix<> W(*Wrow, *Wcol, Wdata);
+ Matrix<> YT(*Yrow, *Ycol, YTdata);
+ Matrix<> XT(*Xrow, *Xcol, XTdata);
+ Matrix<> WT(*Wrow, *Wcol, WTdata);
+ Matrix<> betastart(*Xcol, 1, betastartdata);
+
+ Matrix<> b0(*Xcol, 1, b0data);
+ Matrix<> B0(*Xcol, *Xcol, B0data);
+ Matrix<> R0(*Wcol, *Wcol, R0data);
+
+ Matrix<> Dstart = invpd(R0);
+ Matrix<int> subjectid_mat(*nobs, 1, subjectid);
+ Matrix<int> timeid_mat(*nobs, 1, timeid);
+ Matrix<> subject_groupinfo(*nsubj, 3, subject_groupinfodata);
+ Matrix<> time_groupinfo(*ntime, 3, time_groupinfodata);
+ const int mns = *m + 1;
+ Matrix<> beta_store(*betarow, *betacol);
+ Matrix<> sigma_store(*betarow, mns);
+ Matrix<> D_store(*betarow, *Wcol* *Wcol * mns);
+ double logmarglike;
+ double loglike;
+
+ if (*m == 0){
+ MCMCPACK_PASSRNG2MODEL(GaussianPanelRE_impl,
+ *nsubj, *ntime, *nobs,
+ subjectid_mat, timeid_mat,
+ Y, X, W,
+ *burnin, *mcmc, *thin, *verbose,
+ *sigma2start, Dstart,
+ b0, B0, *c0, *d0, *r0, R0,
+ time_groupinfo, subject_groupinfo,
+ beta_store, sigma_store, D_store,
+ logmarglike, loglike, *chib);
+
+
+ // store marginal likelihood
+ logmarglikeholder[0] = logmarglike;
+ loglikeholder[0] = loglike;
+
+ for (int i=0; i < (*betarow* *betacol); ++i){
+ betadata[i] = beta_store(i);
+ }
+ for (int i=0; i < (*betarow); ++i){
+ sigmadata[i] = sigma_store(i);
+ }
+ for (int i=0; i < (*betarow*mns* *Wcol* *Wcol); ++i){
+ Ddata[i] = D_store(i);
+ }
+ }
+ else {
+ Matrix <> P(mns, mns, Pstart);
+ Matrix<> P_store(*betarow, mns*mns);
+ Matrix<> s_store(*betarow, *ntime*mns);
+ Matrix<> ps_store(*ntime, mns);
+ MCMCPACK_PASSRNG2MODEL(HMMpanelRE_impl,
+ *nsubj, *ntime, *m, *nobs,
+ subjectid_mat, timeid_mat,
+ Y, X, W, YT, XT, WT,
+ *burnin, *mcmc, *thin, *verbose,
+ betastart, *sigma2start, Dstart,
+ b0, B0, *c0, *d0, *r0, R0,
+ P, time_groupinfo, subject_groupinfo,
+ beta_store, sigma_store, D_store, P_store, ps_store, s_store,
+ logmarglike, loglike, *chib);
+ // store marginal likelihood
+ logmarglikeholder[0] = logmarglike;
+ loglikeholder[0] = loglike;
+
+ for (int i=0; i < (*betarow* *betacol); ++i){
+ betadata[i] = beta_store(i);
+ }
+ for (int i=0; i < (*betarow*mns); ++i){
+ sigmadata[i] = sigma_store(i);
+ }
+ for (int i=0; i < (*betarow*mns* *Wcol* *Wcol); ++i){
+ Ddata[i] = D_store(i);
+ }
+ for (int i = 0; i<(*ntime *mns); ++i){
+ psout[i] = ps_store[i];
+ }
+ for (int i = 0; i<(*betarow* *ntime *mns); ++i){
+ sout[i] = s_store[i];
+ }
+
+ }
+ }// end of HMMpanelRE_CC
+
+}// end of extern C
+
+#endif /* HMMpanelRE_CC */
diff --git a/src/MCMCirtKdRob.cc b/src/MCMCirtKdRob.cc
index 6400c59..2bbfd61 100644
--- a/src/MCMCirtKdRob.cc
+++ b/src/MCMCirtKdRob.cc
@@ -42,6 +42,10 @@
#include <R.h> // needed to use Rprintf()
#include <R_ext/Utils.h> // needed to allow user interrupts
+#include <Rdefines.h>
+#include <Rinternals.h>
+
+
typedef Matrix<double,Row,View> rmview;
@@ -385,8 +389,9 @@ static void doubling(double (*logfun)(const double&,
x0 = delta1;
}
else {
- Rprintf("\nERROR: param not in {0,1,2,3} in doubling().\n");
- exit(1);
+ error("ERROR: param not in {0,1,2,3} in doubling().");
+ //Rprintf("\nERROR: param not in {0,1,2,3} in doubling().\n");
+ //exit(1);
}
L = x0 - w*U;
@@ -469,8 +474,9 @@ static void StepOut(double (*logfun)(const double&,
x0 = delta1;
}
else {
- Rprintf("\nERROR: param not in {0,1,2,3} in StepOut().\n");
- exit(1);
+ error("ERROR: param not in {0,1,2,3} in StepOut().");
+ //Rprintf("\nERROR: param not in {0,1,2,3} in StepOut().\n");
+ //exit(1);
}
@@ -628,8 +634,9 @@ static double shrinkageDoubling(double (*logfun)(const double&,
x0 = delta1;
}
else {
- Rprintf("\nERROR: param not in {0,1,2,3} in shrinkageDoubling().\n");
- exit(1);
+ error("ERROR: param not in {0,1,2,3} in shrinkageDoubling().");
+ //Rprintf("\nERROR: param not in {0,1,2,3} in shrinkageDoubling().\n");
+ //exit(1);
}
for (;;){
@@ -714,8 +721,9 @@ static double shrinkageStep(double (*logfun)(const double&,
x0 = delta1;
}
else {
- Rprintf("\nERROR: param not in {0,1,2,3} in shrinkageDoubling().\n");
- exit(1);
+ error("ERROR: param not in {0,1,2,3} in shrinkageDoubling().");
+ //Rprintf("\nERROR: param not in {0,1,2,3} in shrinkageDoubling().\n");
+ //exit(1);
}
for (;;){
diff --git a/src/MCMCregress.cc b/src/MCMCregress.cc
index bcad188..0f3438d 100644
--- a/src/MCMCregress.cc
+++ b/src/MCMCregress.cc
@@ -70,7 +70,7 @@ void MCMCregress_impl (rng<RNGTYPE>& stream, const Matrix<>& Y,
const Matrix<>& B0, double c0, double d0,
unsigned int burnin, unsigned int mcmc, unsigned int thin,
unsigned int verbose, bool chib,
- Matrix<>& result, double& logmarglike, double& loglike)
+ Matrix<>& result, double& logmarglike)
{
// define constants and form cross-product matrices
const unsigned int tot_iter = burnin + mcmc; //total iterations
@@ -151,7 +151,7 @@ void MCMCregress_impl (rng<RNGTYPE>& stream, const Matrix<>& Y,
for (unsigned int i = 0; i < X.rows(); ++i) {
loglike_sum += lndnorm(Y(i), eta(i), sigmastar);
}
- loglike = loglike_sum;
+ double loglike = loglike_sum;
// calculate log prior ordinate
double logprior = 0;
@@ -190,7 +190,7 @@ extern "C" {
const int *b0row, const int *b0col,
const double *B0data, const int *B0row,
const int *B0col, const double *c0, const double *d0,
- double* logmarglikeholder, double* loglikeholder, const int* chib)
+ double* logmarglikeholder, const int* chib)
{
// pull together Matrix objects
Matrix<> Y(*Yrow, *Ycol, Ydata);
@@ -200,13 +200,11 @@ extern "C" {
Matrix<> B0(*B0row, *B0col, B0data);
double logmarglike;
- double loglike;
Matrix<> storagematrix;
MCMCPACK_PASSRNG2MODEL(MCMCregress_impl, Y, X, betastart, b0, B0,
*c0, *d0, *burnin, *mcmc, *thin, *verbose,
- *chib, storagematrix, logmarglike, loglike);
+ *chib, storagematrix, logmarglike);
logmarglikeholder[0] = logmarglike;
- loglikeholder[0] = loglike;
const unsigned int size = *samplerow * *samplecol;
for (unsigned int i = 0; i < size; ++i)
sampledata[i] = storagematrix(i);
diff --git a/src/MCMCresidualBreakAnalysis.cc b/src/MCMCresidualBreakAnalysis.cc
new file mode 100644
index 0000000..e1003bc
--- /dev/null
+++ b/src/MCMCresidualBreakAnalysis.cc
@@ -0,0 +1,517 @@
+// MCMCresidualBreakAnalysis.cc
+//
+// Jong Hee Park
+// Department of Political Science
+// University of Chicago
+// jhp at uchicago.edu
+//
+// 03/03/2009
+
+#ifndef MCMCRESIDUALBREAKANALYSIS_CC
+#define MCMCRESIDUALBREAKANALYSIS_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;
+static double dinvgamma(double theta, double a, double b) {
+ double logf = a * log(b) - lngammafn(a) + -(a+1) * log(theta) +
+ -b/theta;
+ return exp(logf);
+}
+
+
+////////////////////////////////////////////
+// Start MCMCresidualBreakAnalysispoint function
+///////////////////////////////////////////
+template <typename RNGTYPE>
+void MCMCresidualBreakAnalysis_impl(rng<RNGTYPE>& stream,
+ const double m,
+ const Matrix<>& Y, Matrix<>& beta,
+ Matrix<>& Sigma, Matrix<>& P, Matrix<int>& s,
+ const double b0, const double B0,
+ const double c0, const double d0,
+ const Matrix<>& A0,
+ unsigned int burnin, unsigned int mcmc, unsigned int thin,
+ unsigned int verbose, bool chib,
+ Matrix<>& beta_store, Matrix<>& Sigma_store,
+ Matrix<>& P_store, Matrix<>& ps_store, Matrix<int>& s_store,
+ double& logmarglike, double&loglike)
+{
+ // define constants and form cross-product matrices
+ const int tot_iter = burnin + mcmc;
+ const int nstore = mcmc / thin;
+ const int n = Y.rows();
+ const int ns = m + 1;
+ const double B0inv = 1/B0;
+ Matrix<> sigma(ns, 1);
+
+ //MCMC loop
+ unsigned int count = 0;
+ for (int iter = 0; iter < tot_iter; ++iter){
+
+ //////////////////////
+ // 1. Sample beta and Sigma
+ //////////////////////
+ int beta_count = 0;
+ Matrix<int> nstate(ns, 1);
+
+ for (int j = 0; j<ns; ++j){
+ for (int i = 0; i<n; ++i){
+ if (s[i] == j + 1) {
+ nstate[j] = nstate[j] + 1;
+ }// 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);
+ double Bn = 1/(B0 + (double)nstate[j]/Sigma[j]);
+ double bn = Bn*(B0*b0 + sum(yj)/Sigma[j]);
+ beta(j) = stream.rnorm(bn, sqrt(Bn));
+
+ // SIGMA UPDATE
+ double shape = (c0 + (double)nstate[j])/2;
+ const Matrix<> ej(nstate[j], 1);
+ for (int i = 0; i<nstate[j]; ++i){
+ ej(i) = yj(i) - beta(j);
+ }
+ const Matrix<> SSE = crossprod (ej);
+ double scale =(d0 + SSE[0])/2;
+
+ Sigma[j] = 1/stream.rgamma(shape, scale);
+ sigma[j] = sqrt(Sigma[j]);
+ }
+
+ //////////////////////
+ // 2. Sample P
+ //////////////////////
+ double shape1 = 0;
+ double shape2 = 0;
+ P(ns-1, ns-1) = 1;
+
+ for (int j =0; j<(ns-1); ++j){
+ shape1 = A0(j,j) + (double)nstate[j] - 1;
+ shape2 = A0(j,j+1) + 1;
+ P(j,j) = stream.rbeta(shape1, shape2);
+ P(j,j+1) = 1 - P(j,j);
+ }
+
+ //////////////////////
+ // 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);
+ for (int tt=0; tt<n ; ++tt){
+ for (int j = 0; j< ns; ++j){
+ py[j] = dnorm(Y[tt], beta[j], sigma[j]);
+ }
+ if (tt==0) pstyt1 = pr1;
+ else {
+ pstyt1 = ::t(F(tt-1,_)*P);
+ }
+ Matrix<double> unnorm_pstyt = pstyt1%py;
+ const Matrix<double> pstyt = unnorm_pstyt/sum(unnorm_pstyt);
+ for (int j=0; j<ns ; ++j) F(tt,j) = pstyt(j);
+
+ }
+ 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;// 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<ns; ++i)
+ beta_store(count,i) = beta[i];
+ 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,_);
+ s_store(count,_) = s;
+
+ ++count;
+
+ }
+
+
+ // print output to stdout
+ if(verbose > 0 && iter % verbose == 0){
+ Rprintf("\n testpanelSubjectBreak iteration %i of %i \n", (iter+1), tot_iter);
+ for (int j = 0;j<ns; ++j){
+ Rprintf("\n The number of observations in state %i is %10.5f", j+1, static_cast<double>(nstate[j]));
+ }
+ Rprintf("\n beta \n");
+ for (int i = 0; i<ns; ++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){
+
+ Matrix<double> betast = meanc(beta_store);
+ Matrix<double, Row> beta_st(ns, 1);
+ for (int j = 0; j<ns; ++j){
+ beta_st[j] = betast[j];
+ }
+
+ Matrix<double> Sigma_st = meanc(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);
+ 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 double precision = 1.0/Sigma_store(iter, j);
+ const double Bn = 1/(B0 + (double)nstate[j]*precision);
+ double bn = Bn*(B0*b0 + sum(yj)*precision);
+ density_beta(iter, j) = dnorm(beta_st(j), bn, sqrt(Bn));
+ }
+ }// end of pdf.beta
+
+ double pdf_beta = log(prod(meanc(density_beta)));
+
+ //////////////////////
+ // 2. pdf.Sigma
+ //////////////////////
+ Matrix<double> density_Sigma(nstore, ns);
+ for (int iter = 0; iter<nstore; ++iter){
+ 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);
+ for (int tt=0; tt<n ; ++tt){
+ for (int j = 0; j< ns; ++j){
+ py[j] = dnorm(Y[tt], beta_st[j], sigma[j]);
+ }
+ if (tt==0) pstyt1 = pr1;
+ else {
+ pstyt1 = ::t(F(tt-1,_)*P);
+ }
+ Matrix<double> unnorm_pstyt = pstyt1%py;
+ const Matrix<double> pstyt = unnorm_pstyt/sum(unnorm_pstyt);
+ for (int j=0; j<ns ; ++j) F(tt,j) = pstyt(j);
+
+ }
+ ps(n-1,_) = F(n-1,_);
+
+
+ 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
+
+ int beta_count = 0;
+ Matrix<int> nstate(ns, 1);
+
+ for (int j = 0; j <ns ; ++j){
+ for (int i = 0; i<n; ++i){
+ if (s[i] == (j+1)) {
+ nstate[j] = nstate[j] + 1;
+ }// 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);
+ const Matrix<> ej(nstate[j], 1);
+ for (int i = 0; i<nstate[j]; ++i){
+ ej(i) = yj(i) - beta_st(j);
+ }
+ const Matrix<> SSE = crossprod (ej);
+ double scale =(d0 + SSE[0])/2;
+ double shape = (c0 + (double)nstate[j])/2;
+
+ Sigma[j] = stream.rigamma(shape, scale);
+ sigma[j] = sqrt(Sigma[j]);
+ density_Sigma(iter, j) = dinvgamma(Sigma_st[j], shape, scale);
+ }
+
+ double shape1 = 0;
+ double shape2 = 0;
+ P(ns-1, ns-1) = 1;
+
+ for (int j =0; j< (ns-1); ++j){
+ shape1 = A0(j,j) + (double)nstate[j] - 1;
+ shape2 = A0(j,j+1) + 1; //
+ P(j,j) = stream.rbeta(shape1, shape2);
+ P(j,j+1) = 1 - P(j,j);
+ }
+
+ }// 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){
+
+ 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);
+ for (int tt=0; tt<n ; ++tt){
+ for (int j = 0; j< ns; ++j){
+ py[j] = dnorm(Y[tt], beta_st[j], sqrt(Sigma_st[j]));
+ }
+ if (tt==0) pstyt1 = pr1;
+ else {
+ pstyt1 = ::t(F(tt-1,_)*P);
+ }
+ Matrix<double> unnorm_pstyt = pstyt1%py;
+ const Matrix<double> pstyt = unnorm_pstyt/sum(unnorm_pstyt);
+ for (int j=0; j<ns ; ++j) F(tt,j) = pstyt(j);
+
+ }
+ ps(n-1,_) = F(n-1,_);
+
+ 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;// stay
+ }
+ ps(tt,_) = pstyn;
+ --tt;
+ }// end of while loop
+
+ 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; //
+
+ }
+ 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){
+ for (int j = 0; j< ns; ++j){
+ py[j] = dnorm(Y[t], beta_st[j], sqrt(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);
+ }
+
+ 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; //
+ for (int j=0; j<ns ; ++j){
+ density_beta_prior[j] = log(dnorm(beta_st(j), b0, sqrt(B0inv)));
+ }
+ for (int j=0; j<ns ; ++j){
+ density_Sigma_prior[j] = log(dinvgamma(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("\nlogmarglike = %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
+
+}
+
+extern "C"{
+ void MCMCresidualBreakAnalysis(double *betaout, double *Sigmaout, double *psout,
+ const double *Ydata, const int *Yrow, const int *Ycol,
+ const int *m,
+ const int *burnin, const int *mcmc, const int *thin, const int *verbose,
+ const int *uselecuyer, const int *seedarray, const int *lecuyerstream,
+ const double *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 *chib){
+
+ // pull together Matrix objects
+ const Matrix <double> Y(*Yrow, *Ycol, Ydata);
+ const unsigned int tot_iter = *burnin + *mcmc;
+ const unsigned int nstore = *mcmc / *thin;
+ const int n = Y.rows();
+ const int ns = *m + 1;
+
+ // generate starting values
+ Matrix <> beta(ns, 1, betastart);
+ Matrix <> Sigma(ns, 1, Sigmastart);
+ Matrix <> P(ns, ns, Pstart);
+ Matrix <int> s(n, 1, statestart);
+ const Matrix <> A0(ns, ns, A0data);
+ double logmarglike;
+ double loglike;
+
+ // storage matrices
+ Matrix<> beta_store(nstore, ns);
+ Matrix<> Sigma_store(nstore, ns);
+ Matrix<> P_store(nstore, ns*ns);
+ Matrix<> ps_store(n, ns);
+ Matrix<int> s_store(nstore, n);
+ MCMCPACK_PASSRNG2MODEL(MCMCresidualBreakAnalysis_impl, *m,
+ Y, beta, Sigma, P, s,
+ *b0data, *B0data, *c0, *d0, A0,
+ *burnin, *mcmc, *thin, *verbose,
+ *chib,
+ beta_store, Sigma_store,
+ P_store, ps_store, s_store,
+ logmarglike, loglike);
+ logmarglikeholder[0] = logmarglike;
+ loglikeholder[0] = loglike;
+
+ // return output
+ for (int i = 0; i<(nstore*ns); ++i){
+ betaout[i] = beta_store[i];
+ }
+ for (int i = 0; i<(nstore*ns); ++i){
+ Sigmaout[i] = Sigma_store[i];
+ }
+ for (int i = 0; i<(n*ns); ++i){
+ psout[i] = ps_store[i];
+ }
+
+ }// end of MCMCresidualBreakAnalysis
+}//end extern "C"
+
+#endif
+
+
+
+
diff --git a/src/algorithm.h b/src/algorithm.h
index c5bf5ad..d82a4d1 100644
--- a/src/algorithm.h
+++ b/src/algorithm.h
@@ -45,6 +45,11 @@
#include "scythestat/matrix_random_access_iterator.h"
#endif
+#undef DO
+#undef DS
+#undef SO
+#undef SS
+
namespace scythe {
namespace {
typedef unsigned int uint;
diff --git a/src/datablock.h b/src/datablock.h
index 934e74a..0c16a24 100644
--- a/src/datablock.h
+++ b/src/datablock.h
@@ -42,10 +42,15 @@
#ifdef SCYTHE_COMPILE_DIRECT
#include "error.h"
+#include "pthread.h"
#else
#include "scythestat/error.h"
#endif
+#ifdef SCYTHE_PTHREAD
+#include "pthread.h"
+#endif
+
namespace scythe {
/* Convenience typedefs */
namespace { // local to this file
@@ -253,7 +258,13 @@ namespace scythe {
: data_ (0),
block_ (&nullBlock_)
{
+#ifdef SCYTHE_PTHREAD
+ pthread_mutex_lock (&ndbMutex_);
+#endif
block_->addReference();
+#ifdef SCYTHE_PTHREAD
+ pthread_mutex_unlock (&ndbMutex_);
+#endif
}
/* New block constructor: creates a new underlying block of a
@@ -278,7 +289,18 @@ namespace scythe {
: data_ (reference.data_ + offset),
block_ (reference.block_)
{
+#ifdef SCYTHE_PTHREAD
+ bool lock = false;
+ if (block_ == &nullBlock_) {
+ pthread_mutex_lock (&ndbMutex_);
+ lock = true;
+ }
+#endif
block_->addReference();
+#ifdef SCYTHE_PTHREAD
+ if (lock)
+ pthread_mutex_unlock (&ndbMutex_);
+#endif
}
/**** DESTRUCTOR ****/
@@ -287,7 +309,18 @@ namespace scythe {
*/
virtual ~DataBlockReference ()
{
+#ifdef SCYTHE_PTHREAD
+ bool lock = false;
+ if (block_ == &nullBlock_) {
+ pthread_mutex_lock (&ndbMutex_);
+ lock = true;
+ }
+#endif
withdrawReference();
+#ifdef SCYTHE_PTHREAD
+ if (lock)
+ pthread_mutex_unlock (&ndbMutex_);
+#endif
}
protected:
@@ -296,15 +329,33 @@ namespace scythe {
void referenceOther (const DataBlockReference<T_type>& ref,
uint offset = 0)
{
+#ifdef SCYTHE_PTHREAD
+ bool lock = false;
+ if (block_ == &nullBlock_ || ref.block_ == &nullBlock_) {
+ pthread_mutex_lock (&ndbMutex_);
+ lock = true;
+ }
+#endif
withdrawReference ();
block_ = ref.block_;
block_->addReference();
data_ = ref.data_ + offset;
+#ifdef SCYTHE_PTHREAD
+ if (lock)
+ pthread_mutex_lock (&ndbMutex_);
+#endif
}
void referenceNew (uint size)
{
- /* If we are the only referent to this data block, resize it.
+#ifdef SCYTHE_PTHREAD
+ bool lock = false;
+ if (block_ == &nullBlock_) {
+ pthread_mutex_lock (&ndbMutex_);
+ lock = true;
+ }
+#endif
+ /* If we are the only referent to this data block, resize it.
* Otherwise, shift the reference to point to a newly
* constructed block.
*/
@@ -322,12 +373,19 @@ namespace scythe {
data_ = block_->data();
block_->addReference();
}
+#ifdef SCYTHE_PTHREAD
+ if (lock)
+ pthread_mutex_unlock (&ndbMutex_);
+#endif
}
private:
/**** INTERNAL MEMBERS ****/
void withdrawReference ()
{
+ // All calls to withdrawReference are mutex protected and protecting
+ // this too can create a race condition.
+
if (block_->removeReference() == 0
&& block_ != &nullBlock_)
delete block_;
@@ -335,10 +393,17 @@ namespace scythe {
void referenceNull ()
{
+#ifdef SCYTHE_PTHREAD
+ pthread_mutex_lock (&ndbMutex_);
+#endif
withdrawReference();
block_ = &nullBlock_;
block_->addReference();
data_ = 0;
+
+#ifdef SCYTHE_PTHREAD
+ pthread_mutex_unlock (&ndbMutex_);
+#endif
}
@@ -349,6 +414,7 @@ namespace scythe {
private:
DataBlock<T_type>* block_;
static NullDataBlock<T_type> nullBlock_;
+ static pthread_mutex_t ndbMutex_;
}; // end class DataBlockReference
@@ -356,6 +422,13 @@ namespace scythe {
template <typename T>
NullDataBlock<T> DataBlockReference<T>::nullBlock_;
+#ifdef SCYTHE_PTHREAD
+ // mutex initialization
+ template <typename T>
+ pthread_mutex_t
+ DataBlockReference<T>::ndbMutex_ = PTHREAD_MUTEX_INITIALIZER;
+#endif
+
} // end namespace scythe
#endif /* SCYTHE_DATABLOCK_H */
diff --git a/src/error.h b/src/error.h
index b3d7df6..afebdac 100644
--- a/src/error.h
+++ b/src/error.h
@@ -57,6 +57,9 @@
#include <iostream>
#include <vector>
+#include <R.h> // needed to use Rprintf()
+#include <R_ext/Utils.h> // needed to allow user interrupts
+
/*! @cond */
#ifdef SCYTHE_DEBUG_LIB
#define SCYTHE_DEBUG_MSG(MSG) \
@@ -122,12 +125,12 @@
#define SCYTHE_THROW_30(EXCEP,MSG)
#endif
-#define SCYTHE_WARN(MSG) \
- { \
- std::cerr << "WARNING in " << __FILE__ << ", " \
- << __func__ << ", " << __LINE__ << ": " \
- << MSG << "\n"; \
- }
+// #define SCYTHE_WARN(MSG) \
+// { \
+// std::cerr << "WARNING in " << __FILE__ << ", " \
+// << __func__ << ", " << __LINE__ << ": " \
+// << MSG << "\n"; \
+// }
#define SCYTHE_CHECK_WARN(CHECK,MSG) \
{ \
@@ -593,9 +596,9 @@ namespace scythe
// The definition of our terminate handler described above
inline void scythe_terminate ()
{
- std::cerr << serr->what() << std::endl;
- std::cerr << std::endl;
- abort ();
+ // std::cerr << serr->what() << std::endl;
+ // std::cerr << std::endl;
+ Rf_error ("scythe error");
}
} // end namspace SCYTHE
diff --git a/src/matrix.h b/src/matrix.h
index 0432a17..02aa510 100644
--- a/src/matrix.h
+++ b/src/matrix.h
@@ -44,7 +44,7 @@
#ifndef SCYTHE_MATRIX_H
#define SCYTHE_MATRIX_H
-#include <climits>
+//#include <climits>
#include <iostream>
#include <iomanip>
#include <string>
@@ -1888,6 +1888,17 @@ namespace scythe {
(*this, i, 0, i, Base::cols() - 1));
}
+ /*! \brief Returns single element in matrix as scalar type
+ *
+ * This method converts a matrix object to a single scalar
+ * element of whatever type the matrix is composed of. The
+ * method simply returns the element at position zero; if error
+ * checking is turned on the method with throw an error if the
+ * matrix is not, in fact, 1x1.
+ *
+ * \throw scythe_conformation_error (Level 1)
+ */
+
/**** ASSIGNMENT OPERATORS ****/
/*
@@ -3166,7 +3177,7 @@ namespace scythe {
*/
inline void
save (const std::string& path, const char flag = 'n',
- const bool header = false)
+ const bool header = false) const
{
std::ofstream out;
if (flag == 'n') {
diff --git a/src/optimize.h b/src/optimize.h
index 6e9881d..20c76df 100644
--- a/src/optimize.h
+++ b/src/optimize.h
@@ -791,13 +791,13 @@ namespace scythe {
++count;
if (trace){
- std::cout << "BFGS iteration = " << count << std::endl;
- std::cout << "thetamin = " << (t(thetamin)) ;
- std::cout << "gradient = " << (t(fgrad)) ;
- std::cout << "t(gradient) * gradient = " <<
- (t(fgrad) * fgrad) ;
- std::cout << "function value = " << fun(thetamin) <<
- std::endl << std::endl;
+ // std::cout << "BFGS iteration = " << count << std::endl;
+ // std::cout << "thetamin = " << (t(thetamin)) ;
+ // std::cout << "gradient = " << (t(fgrad)) ;
+ // std::cout << "t(gradient) * gradient = " <<
+ // (t(fgrad) * fgrad) ;
+ //std::cout << "function value = " << fun(thetamin) <<
+ // std::endl << std::endl;
}
//std::cout << "Hessian = " << hesscdif(fun, theta) << "\n";
//std::cout << "H = " << H << "\n";
diff --git a/src/rng.h b/src/rng.h
index 11987ff..b9d2d5f 100644
--- a/src/rng.h
+++ b/src/rng.h
@@ -59,6 +59,9 @@
#include "scythestat/la.h"
#endif
+#include <R.h> // needed to use Rprintf()
+#include <R_ext/Utils.h> // needed to allow user interrupts
+
namespace scythe {
/* Shorthand for the matrix versions of the various distributions'
@@ -1040,8 +1043,7 @@ namespace scythe {
}
if (! finite(x)) {
- SCYTHE_WARN("Mean extremely far from truncation point. "
- << "Returning truncation point");
+ Rprintf("Mean extremely far from truncation point. Returning truncation point");
return below;
}
@@ -1096,8 +1098,7 @@ namespace scythe {
+ below;
}
if (! finite(x)) {
- SCYTHE_WARN("Mean extremely far from truncation point. "
- << "Returning truncation point");
+ Rprintf("Mean extremely far from truncation point. Returning truncation point");
return above;
}
@@ -1165,8 +1166,7 @@ namespace scythe {
- below) + below;
}
if (! finite(x)) {
- SCYTHE_WARN("Mean extremely far from truncation point. "
- << "Returning truncation point");
+ Rprintf("Mean extremely far from truncation point. Returning truncation point");
return below;
}
return x;
@@ -1235,8 +1235,7 @@ namespace scythe {
- below) + below;
}
if (! finite(x)) {
- SCYTHE_WARN("Mean extremely far from truncation point. "
- << "Returning truncation point");
+ Rprintf("Mean extremely far from truncation point. Returning truncation point");
return above;
}
return -1*x;
@@ -1277,7 +1276,7 @@ namespace scythe {
for (unsigned int i = 0; i < v; ++i) {
alpha = C * rnorm(Sigma.rows(), 1, 0, 1);
- A = A + (alpha * (t(alpha)));
+ A += (alpha * (t(alpha)));
}
return A;
@@ -1317,7 +1316,7 @@ namespace scythe {
for (ait = alpha.begin_f(); ait != alast; ++ait) {
*yit = rgamma(*ait, 1);
ysum += *yit;
- ++ait;
+ ++yit;
}
y /= ysum;
diff --git a/src/stat.h b/src/stat.h
index 5194539..8160073 100644
--- a/src/stat.h
+++ b/src/stat.h
@@ -224,7 +224,7 @@ namespace scythe {
T
median (const Matrix<T,PO,PS> &A)
{
- Matrix<T, PO, PS> temp(A);
+ Matrix<T, PO, Concrete> temp(A);
uint n = temp.size();
sort(temp.begin(), temp.end());
@@ -271,7 +271,7 @@ namespace scythe {
T
mode (const Matrix<T,PO,PS> &A)
{
- Matrix<T, PO, PS> temp(A);
+ Matrix<T, PO, Concrete> temp(A);
sort(temp.begin(), temp.end());
--
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