Joint pdf for bivariate normal (e.g., Gaussian) Takes means (muX, muY), variances (sigmaX, sigmaY), and correlation (rho) Default parameters gives a “standard” isotropic Gaussian

f = function(x, y, muX=0, muY=0, sigmaX=1, sigmaY=1, rho=0)
{
  K =  1 / (sigmaX * sigmaY * sqrt(1 - rho^2))
  return(K * exp(-0.5 * ((x - muX)^2/sigmaX^2 + (y - muY)^2/sigmaY^2 -
                           2*rho*(x - muX)*(y - muY)/(sigmaX*sigmaY))))
}

List of x, y coordinates

x = seq(-4, 4, 0.1)
y = seq(-4, 4, 0.1)

Create an isotropic Gaussian (uncorrelated and independent)

f.iso = outer(x, y, FUN=f)

Plot wireframe surface

Plot contours

contour(x, y, f.iso, asp=1, )

Create an anisotropic Gaussian (x and y are correlated) Notice that “outer” function lets me pass the parameter “rho” to f

f.aniso = outer(x, y, FUN=f, rho=-0.75)

Plot wireframe surface

persp(x, y, f.aniso, phi=50, theta=20)

Plot contours

contour(x, y, f.aniso, asp=1)

The “rgl” library will let you display interactive plots (you have to install this package first, use install.packages(“rgl”) )

library(rgl)
persp3d(x, y, f.iso, col='darkred')
persp3d(x, y, f.aniso, col='darkred')

Uncorrelated does not imply independence

Here is an example of two random variables (x and y) that are both Gaussian and uncorrelated, but they are NOT independent See http://en.wikipedia.org/wiki/Normally_distributed_and_uncorrelated_does_not_imply_independent

n = 1000
x = rnorm(n)
# w is either -1 or 1, with 0.5 probability for each
w = sample(c(-1,1), n, replace=TRUE)
y = x * w

Histograms of x and y show the marginal densities are Gaussian

hist(x, freq=FALSE)

hist(y, freq=FALSE)

But plotting them jointly in 2D, we see they are not a bivariate Gaussian

plot(x, y)