gibbs_v1<-function(x = c(0, 0, 5), y = c(0, 0, 2.5), rho = 0.5,
        n = 10000, burn = 10,cont=T)
{
# function to produce the sample Gibbs MCMC with contour 
# lines and lines showing values
# for bivariate Normal
  FAC <- sqrt(1 - rho^2)
# draw the contour lines
  if(cont)
      contour(a, b, normgrid(a, b, mu = c(x[2], y[2]),
                   sig = c(x[3], y[3]), rho = rho),
                   levels = c(0.1, 0.05, 0.01, 0.005, 0.001))
# replicate the x and y data to give 1 stream for each output
  x2 <- rep(x, n)
  dim(x2) <- c(3, n)
  x2 <- t(x2)
  y2 <- rep(y, n)
  dim(y2) <- c(3, n)
  y2 <- t(y2)
# burn in
  for(i in 1:burn)
  {
    ynew <- x2[,2] + rnorm(n, (rho * y2[, 3])/x2[, 3]
                 * (x2[, 1] - x2[,2]), FAC * y2[,3])
    xnew <- y2[,2] + rnorm(n, (rho * x2[, 3])/y2[, 3]
                 * (ynew - y2[,2]), FAC * x2[, 3])
    x2[, 1] <- xnew
    y2[, 1] <- ynew
  }
  points(x2, y2)
  list(x2, y2)
}

gibbs_v<-function(x = c(0, 0, 5), y = c(0, 0, 2.5), rho = 0.5,
        n = 10000, burn = 10)
{
# vectorised Gibbs MCMC for bivariate normal
  FAC <- sqrt(1 - rho^2)
# replicate the x and y data to give 1 stream for each output
  x2 <- rep(x, n)
  dim(x2) <- c(3, n)
  x2 <- t(x2)
  y2 <- rep(y, n)
  dim(y2) <- c(3, n)
  y2 <- t(y2)
# burn in
  for(i in 1:burn)
  {
    ynew <- x2[,2] + rnorm(n, (rho * y2[, 3])/x2[, 3]
              * (x2[, 1] - x2[,2]), FAC * y2[,3])
    xnew <- y2[,2] + rnorm(n, (rho * x2[, 3])/y2[, 3]
              * (ynew - y2[,2]), FAC * x2[, 3])
    x2[, 1] <- xnew
    y2[, 1] <- ynew
  }
  list("x"=x2[,1], "y"=y2[,1])
}

gibbs <- function(x = c(10, 0, 5), y = c(-10, 0, 2.5), rho = 0.5, n = 4, burn = 0)
{
# simple Gibbs for bivariate Normal
  FAC <- sqrt(1 - rho^2)
  out <- c(x[1], y[1])
  for(i in 1:n + burn)
  {
    ynew <- rnorm(1, (rho * y[3])/x[3] * x[1], FAC * y[3])
    xnew <- rnorm(1, (rho * x[3])/y[3] * ynew, FAC * x[3])
    if(i > burn)
    {
      x[1] <- xnew
      y[1] <- ynew
      out <- c(out, c(x[1], y[1]))
    }
  }
  dim(out) <- c(2,n+1)
  out
}

"binorm"<-
function(xy, mu = c(0, 0), sig = c(1, 1), rho = 0)
{
# bivariate normal probabilities for contour plotting
  const <- 1/(2 * pi * sqrt(1 - rho^2))
  Q <- (xy[1] - mu[1])^2/sig[1]^2 - (2 * rho * (xy[1] - mu[1]) * (xy[2] -
    mu[2]))/(sig[1] * sig[2]) + (xy[2] - mu[2])^2/sig[2]^2
  prob <- const * exp( - Q/(2 * (1 - rho^2)))
  prob
}

"normgrid"<-
function(x, y, mu = c(0, 0), sigma = c(1, 1), rho = 0)
{
# normal grid generator for contour plotting
  a <- length(x)
  b <- length(y)
  d <- a * b
  grid <- rep(0, d)
#
  dim(grid) <- c(b, a)
  for(i in 1:a)
  {
    for(j in 1:b)
    {
      grid[j, i] <- binorm(c(y[j], x[i]), mu, sigma, rho)
    }
  }
  grid
}

f <- function(theta, mu = 3.508, sigma = 0.01, n = 10)
{
#normal likelihood for sim_full
  z <- exp(-n/2/sigma^2 *(mu - theta)^2)
  z
}
 
sim_full <- function (theta=3.5, n = 10, w = 1000, l=f, p=dnorm)
{
# Metropolis Hastings MCMC with symmetric transition probability
# w is burn  in n is number of samples theta is start value
# l is likelihood function p is prior
  z <- rep(0, n)
  q <- rep(0, w + n)
  for(i in 1 : (w + n))
  {
    new <- theta + rnorm(1, 0, 0.01)
    prob <- (p(new) * l(new, n=n)) 
    prob <- prob / (p(theta) * l(theta, n=n))
    alpha <- min(prob, 1)
    q[i] <- alpha
    if(runif(1) < alpha) theta <- new
    if(i > w) z[i - w] <- theta
  }
  z
}