############################################################################################ # Example: COAL MINING DISASTERS, POISSON CHANGEPOINT PROBLEM # ############################################################################################ # # # This archive is part of the free distribution of data and statistical software code for # # "Bayesian Methods: A Social and Behavioral Sciences Approach" by Jeff Gill (c) 2002. # # You are free to use, modify, distribute, publish, etc. provided attribution. Please # # forward bugs, complaints, comments, and useful changes to: jgill@ucdavis.edu. # # # ############################################################################################ # # # This is an estimation of the coal mining disasters Poisson changepoint. A very well- # # studied dataset, slightly simplified here. # # # ############################################################################################ # WRITE OUT THE DATA coal.mining.disasters <- c(4,5,4,1,0,4,3,4,0,6,3,3,4,0,2,6,3,3,5,4,5,3,1,4,4,1,5,5,3,4,2,5,2,2,3,4,2,1, 3,2,2,1,1,1,1,3,0,0,1,0,1,1,0,0,3,1,0,3,2,2, 0,1,1,1,0,1,0,1,0,0,0,2,1,0,0,0, 1,1,0,2,3,3,1,1,2,1,1,1,1,2,4,2,0,0,0,1,4,0,0,0,1,0,0,0,0,0,1,0,0,1,0,1) # A SIMPLE, CUSTOMIZED GIBBS SAMPLER FOR THIS PROBLEM gibbs3 <- function(theta.matrix,y,reps) { # alpha <- 4; beta <- 1; gamma <- 1; delta <- 2 alpha <- 2; beta <- 2; gamma <- 2; delta <- 2 for (i in 2:(reps+1)) { lambda <- rgamma(1,alpha+sum(y[1:theta.matrix[(i-1),3]]),(beta+theta.matrix[(i-1),3])) phi <- rgamma(1,gamma+sum(y[theta.matrix[(i-1),3]:length(y)]),(delta+length(y)-theta.matrix[(i-1),3])) m <- round(rexp(1,lambda+phi) + beta*lambda + phi*delta + length(y)*phi/delta) theta.matrix <- rbind(theta.matrix,c(lambda,phi,m)) } theta.matrix } # RUN THE GIBBS SAMPLER start <- matrix(c(1,1,100),1,3) coal.out <- gibbs3(start,coal.mining.disasters,1000) summary(coal.out) # THIS IS A FUNCTION TO PLOT THE GIBBS SAMPLER 1 DIMENSION AT A TIME UPDATING plot.walk.multi <- function(walk.mat,sim.rm) { walk.mat <- walk.mat[,-sim.rm] points(walk.mat[1,1],walk.mat[1,2]) for(i in 1:(nrow(walk.mat)-1)) { segments(walk.mat[i,1],walk.mat[i,2],walk.mat[(i+1),1],walk.mat[i,2],col="blue") segments(walk.mat[(i+1),1],walk.mat[i,2],walk.mat[(i+1),1],walk.mat[(i+1),2],col="blue") } } # DO THE PLOT FOR A SMALL NUMBER OF UPDATES. THIS IS REALLY FUN TWO WATCH. start <- matrix(c(1,1,100),1,3) coal.out <- gibbs3(start,coal.mining.disasters,100) par(mfrow=c(1,2),mar=c(3,2,2,2),oma=c(3,3,1,1)) plot(coal.out,xlim=range(coal.out[,1]),ylim=range(coal.out[,3]),type="n",xlab="",ylab="",font=2) # Main title mtext("Gibbs Sampling Path\nk versus Lambda",side=3,line=0.50,cex=1.2,font=2) mtext(outer=F,side=1,cex=1.1,"Lambda",line=2,font=2) mtext(outer=F,side=2,"k",line=2.25,font=2) plot.walk.multi(coal.out,2) plot(coal.out,xlim=range(coal.out[,2]),ylim=range(coal.out[,3]),type="n",xlab="",ylab="",font=2) mtext("Gibbs Sampling Path\nk versus Phi",side=3,line=0.50,cex=1.2,font=2) mtext(outer=F,side=1,cex=1.1,"Phi",line=2.25,font=2) mtext(outer=F,side=2,cex=1.1,"k",line=2,font=2) plot.walk.multi(coal.out,1) summary(coal.out)