## Markov chain Monte Carlo algorithm for a Bayesian (single) change point model ## read in the data ## chptdat = read.table("chpt.dat",header=T) ## Y = chptdat\$Ener ## chptdat = read.table("coal.dat",header=T) ## Y = chptdat\$Deaths 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(NUMIT=1000,dat=Y) { n = length(dat) cat("n=",n,"\n") ## set up ## NUMIT x 5 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, NUMIT) acc = 0 # count number of accepted proposals (for k only) ## 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:NUMIT) { ## 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)} ## 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))