## Bivariate Gaussian distribution for varying covariance matrix
library(animation)
library(lattice)
pmgauss <- function(x, mu, Sigma) {
    dim <- dim(mu)
    apply(x, 1, function(xj) {
        if (length(xj) != dim) 
            stop("x and mu have incompatible dimensions\n")
        1/((2 * pi)^dim * sqrt(det(Sigma)^2)) * exp(-0.5 * 
            t(xj - mu) %*% solve(Sigma) %*% (xj - mu))
    })
}
plot.bigauss <- function(mu, Sigma) {
    x <- seq(-4, 4, 0.1)
    y <- x
    z <- outer(x, y, function(a, b) pmgauss(as.matrix(cbind(a, 
        b)), mu, Sigma))
    par(mfrow = c(1, 2))
    persp(x, y, z, phi = 30, theta = 30, xlim = c(-3, 3), 
        ylim = c(-3, 3), r = 8, box = FALSE, border = NA, 
        shade = 0.5, col = "yellow")
    mtext(side = 3, line = -2, bquote(Sigma == bgroup("[", 
        atop(list(.(Sigma[1, 1]), .(Sigma[2, 1])), list(.(Sigma[1, 
            2]), .(Sigma[2, 2]))), "]")))
    contour(x, y, z, cex = 0.9)
}
mu <- array(c(0, 0))
Sigmas <- list(matrix(c(0.5, 0, 0, 0.5), ncol = 2), matrix(c(1, 
    0, 0, 1), ncol = 2), matrix(c(1, 0.5, 0.5, 1), ncol = 2), 
    matrix(c(1, 0.8, 0.8, 1), ncol = 2))
i = 1
while (i <= ani.options("nmax") & i <= length(Sigmas)) {
    plot.bigauss(mu, Sigmas[[i]])
    ani.pause()
    i = i + 1
}
## R version 2.12.0 (2010-10-15)
## Platform: i686-pc-linux-gnu (32-bit)
## Other packages: animation 2.0-3, lattice 0.19-13