## Markov chain Monte Carlo algorithm for a Bayesian (single) change point model
## Originally written by Murali Haran
KGUESS = 10 # our guess for k based on exploratory data analysis
## Note: this function is not written in the most efficient way since its
## purpose is primarily instructive
mhsampler <- function(dat, MCMCiterations=1000, kfixed=FALSE)
{
n <- length(dat)
## Set up 5 x MCMCiterations matrix to store Markov chain values.
## Each row corresponds to one of 5 parameters in order: theta,lambda,k,b1,b2
## Each column corresponds to a single state of the Markov chain
mchain <- matrix(NA, 5, MCMCiterations, dimnames=list(
c("theta", "lambda", "k", "b1", "b2"), NULL))
acc <- 0 # count number of accepted proposals (for k only; others are Gibbs
# sampled so they're always accepted)
## starting values for Markov chain
## This is somewhat arbitrary but any method that produces reasonable
## values for each parameter is usually adequate.
## For instance, we can use approximate prior means or approximate MLEs.
kinit <- floor(n/2) # approximately halfway between 1 and n
mchain[,1] <- c(1,1,kinit,1,1)
for (i in 2:MCMCiterations)
{
## most upto date state for each parameter
currtheta <- mchain[1,i-1]
currlambda <- mchain[2,i-1]
currk <- mchain[3,i-1]
currb1 <- mchain[4,i-1]
currb2 <- mchain[5,i-1]
## sample from full conditional distribution of theta (Gibbs update)
currtheta <- rgamma(1,shape=sum(Y[1:currk])+0.5,
scale=currb1/(currk*currb1+1))
## sample from full conditional distribution of lambda (Gibbs update)
currlambda <- rgamma(1,shape=sum(Y[(currk+1):n])+0.5,
scale=currb2/((n-currk)*currb2+1))
## sample from full conditional distribution of k (Metropolis-Hastings update)
propk <- sample(x=seq(2,n-1), size=1) # draw one sample at random from uniform{2,..(n-1)}
if (kfixed) {
currk <- KGUESS
} else {
## Metropolis accept-reject step (in log scale)
logMHratio <- sum(Y[1:propk])*log(currtheta)+sum(Y[(propk+1):n])*
log(currlambda)-propk*currtheta- (n-propk)*currlambda -
(sum(Y[1:currk])*log(currtheta)+sum(Y[(currk+1):n])*
log(currlambda)-currk*currtheta- (n-currk)*currlambda)
logalpha <- min(0,logMHratio) # alpha = min(1,MHratio)
if (log(runif(1))