############################################################################################
#
# Example 7.15 (pages 284) – Figure 7.8 (page 285)  -  Slice sampling 
#
############################################################################################
#
# Author : Hedibert Freitas Lopes
#          Graduate School of Business
#          University of Chicago
#          5807 South Woodlawn Avenue
#          Chicago, Illinois, 60637
#          Email : hlopes@ChicagoGSB.edu
#
############################################################################################
f  =  function(theta){(2+theta)^14*(1-theta)^3*theta^5}
Sy  =  function(y){
  ind  =  f(theta)>y
  return(c(min(theta[ind]),max(theta[ind])))
}
theta = seq(0,1,length=1000)

set.seed(1244)
M   = 5000
th  = 0.4
y   = runif(1,0,f(th))
ths = c(th,y) 
for (i in 1:M){
  int = Sy(y)
  th  = runif(1,int[1],int[2])
  y   = runif(1,0,f(th))
  ths = rbind(ths,c(th,y))
}
xs  =  ths
xss = matrix(0,2*M,2)
xss[1,] = xs[1,1:2]
xss[2,2] = xs[1,2]
for (i in 2:M){
  xss[2*(i-1),1] = xs[i,1]
  xss[2*i-1  ,1] = xs[i,1]
  xss[2*i-1,2] = xs[i,2]
  xss[2*i  ,2] = xs[i,2]
}

a = 0.0005716342
xss[,2]=a*xss[,2]
par(mfrow=c(2,2))
plot(theta,a*f(theta),xlab="",ylab="",type="l",main="(a)",axes=F)
axis(1);axis(2)
lines(xss[1:20,])
plot(theta,a*f(theta),xlab="",ylab="",type="l",main="(b)",axes=F)
axis(1);axis(2)
lines(xss[1:500,])
plot(ths[1:500,1],type="l",xlab="iterations",ylab="",main="(c)",ylim=c(0,1),axes=F)
axis(1);axis(2)
hist(ths[,1],xlab="",ylab="",prob=T,main="(d)",breaks=seq(0,1,length=30))
lines(theta,a*f(theta))


c(mean(ths[,1]),sqrt(var(ths[,1])),
quantile(ths[,1],0.025),quantile(ths[,1],0.975))