MCMCbinaryChange | R Documentation |
This function generates a sample from the posterior distribution of a binary model with multiple changepoints. The function uses the Markov chain Monte Carlo method of Chib (1998). The user supplies data and priors, and a sample from the posterior distribution is returned as an mcmc object, which can be subsequently analyzed with functions provided in the coda package.
MCMCbinaryChange(
data,
m = 1,
c0 = 1,
d0 = 1,
a = NULL,
b = NULL,
burnin = 10000,
mcmc = 10000,
thin = 1,
verbose = 0,
seed = NA,
phi.start = NA,
P.start = NA,
marginal.likelihood = c("none", "Chib95"),
...
)
data |
The data. |
m |
The number of changepoints. |
c0 |
|
d0 |
|
a |
|
b |
|
burnin |
The number of burn-in iterations for the sampler. |
mcmc |
The number of MCMC iterations after burn-in. |
thin |
The thinning interval used in the simulation. The number of MCMC iterations must be divisible by this value. |
verbose |
A switch which determines whether or not the
progress of the sampler is printed to the screen. If
|
seed |
The seed for the random number generator. If NA, current R system seed is used. |
phi.start |
The starting values for the mean. The default value of NA will use draws from the Uniform distribution. |
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 |
marginal.likelihood |
How should the marginal likelihood be
calculated? Options are: |
... |
further arguments to be passed |
MCMCbinaryChange
simulates from the posterior distribution
of a binary model with multiple changepoints.
The model takes the following form:
Y_t \sim
\mathcal{B}ernoulli(\phi_i),\;\; i = 1, \ldots, k
Where k
is the number of states.
We assume Beta priors for \phi_{i}
and for transition
probabilities:
\phi_i \sim \mathcal{B}eta(c_0, d_0)
And:
p_{mm} \sim \mathcal{B}eta{a}{b},\;\; m = 1, \ldots, k
Where
M
is the number of states.
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 prob.state
storage matrix that contains the probability of state_i
for
each period, and the log-marginal likelihood of the model
(logmarglike
).
Jong Hee Park. 2011. “Changepoint Analysis of Binary and Ordinal Probit Models: An Application to Bank Rate Policy Under the Interwar Gold Standard." Political Analysis. 19: 188-204. <doi:10.1093/pan/mpr007>
Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. “MCMCpack: Markov Chain Monte Carlo in R.”, Journal of Statistical Software. 42(9): 1-21. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.18637/jss.v042.i09")}.
Siddhartha Chib. 1995. “Marginal Likelihood from the Gibbs Output.” Journal of the American Statistical Association. 90: 1313-1321. <doi: 10.1080/01621459.1995.10476635>
MCMCpoissonChange
,plotState
,
plotChangepoint
## Not run:
set.seed(19173)
true.phi<- c(0.5, 0.8, 0.4)
## two breaks at c(80, 180)
y1 <- rbinom(80, 1, true.phi[1])
y2 <- rbinom(100, 1, true.phi[2])
y3 <- rbinom(120, 1, true.phi[3])
y <- as.ts(c(y1, y2, y3))
model0 <- MCMCbinaryChange(y, m=0, c0=2, d0=2, mcmc=100, burnin=100, verbose=50,
marginal.likelihood = "Chib95")
model1 <- MCMCbinaryChange(y, m=1, c0=2, d0=2, mcmc=100, burnin=100, verbose=50,
marginal.likelihood = "Chib95")
model2 <- MCMCbinaryChange(y, m=2, c0=2, d0=2, mcmc=100, burnin=100, verbose=50,
marginal.likelihood = "Chib95")
model3 <- MCMCbinaryChange(y, m=3, c0=2, d0=2, mcmc=100, burnin=100, verbose=50,
marginal.likelihood = "Chib95")
model4 <- MCMCbinaryChange(y, m=4, c0=2, d0=2, mcmc=100, burnin=100, verbose=50,
marginal.likelihood = "Chib95")
model5 <- MCMCbinaryChange(y, m=5, c0=2, d0=2, mcmc=100, burnin=100, verbose=50,
marginal.likelihood = "Chib95")
print(BayesFactor(model0, model1, model2, model3, model4, model5))
## plot two plots in one screen
par(mfrow=c(attr(model2, "m") + 1, 1), mai=c(0.4, 0.6, 0.3, 0.05))
plotState(model2, legend.control = c(1, 0.6))
plotChangepoint(model2, verbose = TRUE, ylab="Density", start=1, overlay=TRUE)
## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.