# Graphical respresentation of MCMC (Figure 6.2) # Ardo, UCL 2016 # Using 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))