# Joint pdf for bivariate normal # 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 persp(x, y, f.iso, phi=50, theta=20) # 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) library(rgl) persp3d(x, y, f.iso, col='darkred') persp3d(x, y, f.aniso, col='darkred') # 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)