Gibbs Sampler for Bivariate Normal distribution

[GaryBAYLOR]

This is a realization of the Gibbs sampler for bivariate normal distribution from Gelman's book Bayesian Data Analysis (Third edition), Page 277.

Gibbs.bivariate.normal <- function(n = 1e3, burn = 1e2, rho = .8, y1 = 0, y2 = 0) {
    N <- n + burn
    theta1 <- numeric(N+1)
    theta2 <- numeric(N)
    theta1[1] <- 2.5

    for(i in 1:N) {
        theta2[i] <- rnorm(1, mean = y2 + rho *(theta1[i] - y1), sd = sqrt(1 - rho^2))
        theta1[i+1] <- rnorm(1, mean = y1 + rho *(theta2[i] - y2), sd = sqrt(1 - rho^2))
    }

    theta1 <- theta1[(burn+1):N]
    theta2 <- theta2[(burn+1):N]
    plot(theta1, theta2, pch = 19, cex = .3)
}

When rho = 0.7, the plot is like this When rho = -0.3, the plot is like this This is when rho = -0.9

Gibbs Sampler is really easy to implement and handy to use!