## 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))