# Graphical respresentation of MCMC (Figure 6.2)

# Ardo, UCL 2016
 

# Using <Mbeta> from script for MCMC 
# dim(Mbeta) is 2 X 1000 X 2
# chain 1 is Mbeta[1,,1:2] 
# chain 2 is Mbeta[2,,1:2]

# Load the chains:
load("Ch6_Mbeta.RData")
M <- dim(Mbeta)[2]
burn.in <- M/2

# Inits graphs:
opar <- par(mfrow=c(1,2), mex=0.8,mar=c(5,5,2,2))
lwd  <- 2
col  <- 2

##############
# Plot chains:
plot(c(-3,-1),c(-1,-7),type="n",
      ylab=expression(beta[13]),xlab=expression(beta[12]))
thin <- 10
indx <- 1
indy <- 2
for(m in seq(from=1,to=M-thin,by=thin)){
   # Chain 1:
   x1 <- Mbeta[1,m,indx]    
   y1 <- Mbeta[1,m,indy]
   points(x1,y1,cex=.8,col=1,pch=16)
   x0 <- Mbeta[1,m+thin,indx]
   y0 <- Mbeta[1,m+thin,indy]
   lines(c(x0,x1),c(y0,y1),lwd=lwd,col=1)
   # Chain 2:
   x1 <- Mbeta[2,m,indx]    
   y1 <- Mbeta[2,m,indy]
   points(x1,y1,cex=.8,col="grey",pch=16)
   x0 <- Mbeta[2,m+thin,indx]
   y0 <- Mbeta[2,m+thin,indy]
   lines(c(x0,x1),c(y0,y1),lwd=lwd,col="grey")
}

#################
# Plot densities:
lwd <- 4
chainx <- density(c(Mbeta[1,burn.in:M,indx], Mbeta[2,burn.in:M,indx]),bw="nrd0")#,bw=.15)
chainy <- density(c(Mbeta[1,burn.in:M,indy], Mbeta[2,burn.in:M,indy]),bw="nrd0")#,bw=.15)
# Setting up plot region:
plot(c(-3,-1),c(-1,-7),type="n",ylab=expression(beta[13]),
     xlab=expression(beta[12]), axes=FALSE)
Axis(side=1)
Axis(side=2)
# Plotting density first parameter (horizontal orientation)
# With some re-scaling...:
lines(chainx$x,-7+.5*chainx$y,lwd=lwd,col=1)
# Plotting density 2nd parameter (vertical orientation)
# With some re-scaling and only from -7 upwards...:
L  <- length(chainy$x)
xx <- -3+.5*chainy$y
yy <- chainy$x
x0 <- xx[1]
y0 <- yy[2]
for(i in 2:L){
 x00 <- xx[i]
 y00 <- yy[i]
 if(y0 > -7){lines(c(x0,x00),c(y0,y00),lwd=lwd,col=1)}
 x0 <- x00
 y0 <- y00
}
 
# Tidy up:
opar <- par(mfrow=c(1,2), mex=0.8,mar=c(5,5,2,2))