rm(list=ls())
quant025=function(x){quantile(x,0.025)}
quant975=function(x){quantile(x,0.975)}

# full conditional of m
full = function(m,lambda,phi,y,n,alpha,beta,gamma,delta){
  (alpha-1+ifelse(m>1,sum(y[1:m]),0))*log(lambda)-(beta+m)*lambda+
  (gamma-1+ifelse(m<n,sum(y[(m+1):n]),0))*log(phi)-(delta+n-m)*phi
}

set.seed(2125235)
ns     = c(50,50,50)
lambda = c(10,5,2)
y = rpois(ns[1],lambda[1])
y = c(y,rpois(ns[2],lambda[2]))
y = c(y,rpois(ns[3],lambda[3]))
n = length(y)

# hyperparameters
# ---------------
a1 = 0.001
b1 = 0.001
a2 = 0.001
b2 = 0.001
a3 = 0.001
b3 = 0.001

#########################################
# Gibbs Sampler
#########################################
set.seed(123456)
m1     = 50       
m2     = 100
M0     = 1000      
M      = 10000   
niter  = M+M0
draws  = matrix(0,niter,5)
for (iter in 1:niter){
  s1 = sum(y[1:m1])
  s2 = sum(y[(m1+1):m2])
  s3 = sum(y[(m2+1):n])
  lambda1 = rgamma(1,a1+s1,b1+m1)	
  lambda2 = rgamma(1,a2+s2,b1+m2-m1)	
  lambda3 = rgamma(1,a3+s3,b1+n-m2)	
  fulls  = rep(0,m2)
  for (j in 1:m2)
    fulls[j] = full(j,lambda1,lambda2,y[1:m2],m2,a1,b1,a2,b2)
  fulls = exp(fulls-max(fulls))
  m1     = sample(1:m2,size=1,prob=fulls)
  fulls  = rep(0,n-m1)
  for (j in 1:(n-m1))
    fulls[j] = full(j,lambda2,lambda3,y[(m1+1):n],n-m1,a2,b2,a3,b3)
  fulls = exp(fulls-max(fulls))
  m2 = sample((m1+1):n,size=1,prob=fulls)
  draws[iter,] = c(lambda1,lambda2,lambda3,m1,m2)
}
draws = draws[(M0+1):niter,]

# Posterior summaries
# -------------------
pdf(file="workedexample-poissonwith2changepoints.pdf",width=12,height=8)
plot(y,type="h",xlab="",ylab="")
points(ns[1],0,col=2,pch=16,cex=2)
points(ns[1]+ns[2],0,col=2,pch=16,cex=2)

ind   = seq(1,M,by=M/1000)
par(mfrow=c(3,5))
plot(ind,draws[ind,1],xlab="iteration",ylab="",main=expression(lambda[1]),type="l")
abline(h=lambda[1],col=2)
plot(ind,draws[ind,2],xlab="iteration",ylab="",main=expression(lambda[2]),type="l")
abline(h=lambda[2],col=2)
plot(ind,draws[ind,3],xlab="iteration",ylab="",main=expression(lambda[3]),type="l")
abline(h=lambda[3],col=2)
plot(ind,draws[ind,4],xlab="iteration",ylab="",main=expression(m[1]),type="l")
abline(h=ns[1],col=2)
plot(ind,draws[ind,5],xlab="iteration",ylab="",main=expression(m[2]),type="l")
abline(h=ns[1]+ns[2],col=2)
for (i in 1:5)
  acf(draws[,i],main="")
for (i in 1:3){
  hist(draws[,i],prob=TRUE,xlab="",main="")
  points(lambda[i],0,pch=16,col=2)
}
plot(table(draws[,4])/M,type="h",xlab="",main="",ylab="Probability")
points(ns[1],0,pch=16,col=2)
plot(table(draws[,5])/M,type="h",xlab="",main="",ylab="Probability")
points(ns[1]+ns[2],0,pch=16,col=2)
dev.off()