# changepoint-sim.R
#
# Simulate and estimate a simple changepoint example using MCMCpack
#
# Patrick T. Brandt
#
# 20100209: Initial version
#

# Load the MCMCpack library
library(MCMCpack)

# Simulate a time series of counts with two changes in the mean
set.seed(10805)
n <- 50
y <- c(rpois(n, 3), rpois(n, 1), rpois(n, 10))

# Plot the data
par(mfrow=c(1,2))
plot(ts(y), ylab="Counts")
abline(v=50)
abline(v=100)
plot(cumsum(ts(y)), ylab="Cumulative Counts", xlab="Time", type="l")
abline(v=50)
abline(v=100)

# Estimate the models via MCMC sampling
m1 <- MCMCpoissonChange(y ~1, m=1, c0=1, d0=1, marginal.likelihood="Chib95")
m2 <- MCMCpoissonChange(y ~1, m=2, c0=1, d0=1, marginal.likelihood="Chib95")
m3 <- MCMCpoissonChange(y ~1, m=3, c0=1, d0=1, marginal.likelihood="Chib95")

#Summarize the means for each changepoint
summary(m1)
summary(m2)
summary(m3)

# Compute the Bayes factors for the models
print(BayesFactor(m1, m2, m3))

# Plot the locations of the changepoints and their densities
plotState(m2)