#install.packages("ellipse")
library(ellipse)

### Gibbs sampling when the target dist. is bivariate normal,
### correlation rho,standard marginals
### initial value is far away

rho <- 0.99  
wait <- 0.25 ### in seconds
nrep <- 250  ### number of iterations
lft <- -2; rht <- 2  ### location of "iteration bar"

set.seed(13)

th <- c(-10,-10)   ### far away to start!

plot(th[1],th[2],pch=19, col="blue", xlim=c(-10.5,2.5), ylim=c(-10.5,2.5))
points(ellipse(rho), type="l")  ### show where target is

points(c(lft,rht),c(-8.5,-8.5),type="l")  ### "iteration bar" 
text(lft,-9,"0")
text(rht,-9,as.character(nrep))

Sys.sleep(20)

opt <- matrix(NA, nrep, 2)
for (i in 1:nrep) {
  Sys.sleep(wait)
  opt[i,] <- th   ### store output at each iteration
  th[1] <- rnorm(1, rho*th[2], sqrt(1-rho^2))   ### 1st full conditiona1 
  th[2] <- rnorm(1, rho*th[1], sqrt(1-rho^2))   ### 2nd full conditional
  points(th[1],th[2], pch=19, col="blue")
  points(lft+(i/nrep)*(rht-lft),-8)
 
 }   

#### now show traceplots
par(ask=T)
par(mfrow=c(1,2))
plot(1:nrep, opt[,1], type="l", ylim=c(-10,3),xlab="theta1")
plot(1:nrep, opt[,2], type="l", ylim=c(-10,3),xlab="theta2")
par(ask=F); par(mfrow=c(1,1))