Plotting the Bivariate Gaussian Density

...In which I decide to heed to a friend's advise and use MATLAB to generate some contour plots for the bivariate normal distribution and obtain pretty-looking results...

Cute(r) Plots for the Gaussian Density

The following plots of the bivariate normal (Gaussian) density function were generated in MATLAB. The script follows, together with as a small numerical example on how to compute the density.

The PDF

The probability density function (pdf) for the multivariate ($D$-variate) Gaussian is this: $$p(\mathbf{x}|\mu,\Sigma) = \frac{1}{(2\pi)^{D/2}|\Sigma|^{1/2}}\exp\left\{-\frac{1}{2}(\mathbf{x}-\mu)^T\Sigma^{-1}(\mathbf{x}-\mu)\right\}$$ where $\mu$ is the $D$x1 vector that describes the position of the distribution, and $\Sigma$ is a symmetric $D$x$D$ positive definite matrix (i.e. $\mathbf{u^T} \Sigma \mathbf{u}$ is a positive real number for any $D$-dimensional vector $\mathbf{u}$) that describes the shape of the distribution.

In our bivariate case $D$ is 2. The locus of points in the plane each of which gives the same fixed value of the density is an ellipse.

If this sounds like a mouthful, don't despair, a 2-$D$ numeric example will clarify things. It's really much less complicated than it looks.

An Example

Suppose for example that $\mu = \left[ \begin{array} {c} 1 \\ 0 \end{array} \right]$ and $\Sigma = \left[ \begin{array}{cc} 0.8 & 0.7 \\ 0.7 & 1.3 \end{array} \right] $. Evaluate the density at some point $[x, y]$, one step at a time.

Consider the expression in the exponent first and ask MATLAB to compute the inverse of $\Sigma$ (by asking it to evaluate inv(sigma), see code that follows in a minute). MATLAB returns $\Sigma^{-1} = \left[ \begin{array}{cc} 2.36 & -1.27 \\ -1.27 & 1.45 \end{array} \right] $. Thus the expression in the exponent is $$ \left[ \begin{array} {r} x - 1, y - 0 \end{array} \right] \left[ \begin{array}{cc} 2.36 & -1.27 \\ -1.27 & 1.45 \end{array} \right] \left[ \begin{array} {c} x - 1 \\ y - 0 \end{array} \right] $$ So if you were to evaluate this at $[x, y] = [3, 1]$, for instance, then $$ \left[ \begin{array} {r} 2, 1 \end{array} \right] \left[ \begin{array}{cc} 2.36 & -1.27 \\ -1.27 & 1.45 \end{array} \right] \left[ \begin{array} {c} 2 \\ 1 \end{array} \right] = \left[ \begin{array} {r} 2, 1 \end{array} \right] \left[ \begin{array} {c} 3.45 \\ -1.09 \end{array} \right] = 5.82 $$

The expression that multiplies $e^{-5.82}$ is just a constant ($| \Sigma |$ is the determinant of $\Sigma$). To plot the contours, we don't really need this exact constant, but let's compute it anyway - it's really quick in MATLAB. When asked to compute det(sigma), MATLAB returns 0.55. Thus, we finally have: $$ \frac{1} {2 \pi \sqrt{0.55}} \exp(-5.82) = \text{close enough to zero} $$ So this bivariate density evaluates to about 0 at $[3, 1]^T$, which is in the black region in the figure above.

Different Types of the Covariance Matrix

The previous example looked at a full covariance matrix all of whose off-diagonal entries were non-zero. The random variables were dependent and the iso-probability contour lines were ellipsoids whose axes were not aligned with the coordinate axes in any special way.

If the random variables $X$ and $Y$ are independent, then their covariance is zer0 and the off-diagonal entries in $\Sigma$ are equal to zero.

In general, the covariance matrices in the bivariate case take three forms: $$ \Sigma_{spher} = \left[ \begin{array}{cc} \sigma^2 & 0 \\ 0 & \sigma^2 \end{array} \right], \Sigma_{diag} = \left[ \begin{array}{cc} \sigma_1^2 & 0 \\ 0 & \sigma_2^2 \end{array} \right], \Sigma_{full} = \left[ \begin{array}{cc} \sigma_{11}^2 & \sigma_{12}^2 \\ \sigma_{21}^2 & \sigma_{22}^2 \end{array} \right] $$

And this is what their pictures look like:

Notice that the spherical covariance matrices are multiples of the identity matrix. The variables are independent and the equi-probability contour lines are perfect circles. On the other hand, diagonal covariance matrices permit distinct non-zero entries on the main diagonal. The variables are still independent but scaled differently and the equi-probability contour lines are ellipses (hyper-ellipsoids in higher dimensional spaces) and the principal axes are aligned with the coordinate axes.

The full covariance matrices I considered in the beginning of this blog are symmetric positive definite. Variables there are dependent and the equi-probability contour lines are ellipsoids whose principal axes are not aligned with the coordinate axes.

The Code

Finally, here's the code that plotted the figure:

  % Plots the Bivariate Normal Density Function for Various Parameters Values
  % author: David Dobor 

  clear all;close all;

  %% Define a Gaussian with a Full Covariance Matrix 
  mu = [1;2];
  sigma = [0.8 0.7;0.7 1.3];
  % Define the grid for visualization
  [X,Y] = meshgrid(-3:0.07:4,-2:0.07:5);
  % Define the constant
  const = (1/sqrt(2*pi))^2;
  const = const/sqrt(det(sigma));
  % Compute Density at every point on the grid
  temp = [X(:)-mu(1) Y(:)-mu(2)];
  pdf = const*exp(-0.5*diag(temp*inv(sigma)*temp'));
  pdf = reshape(pdf,size(X));

  % plot the result 
  figure(1)
  surfc(X, Y, pdf, 'LineStyle', 'none');

  %add some info to the plot
  mu_str = '$$\mu = \left[ \begin{array} {c} 1 \\  2 \end{array} \right]$$';
  text('Interpreter','latex',...
  'String',mu_str,...
  'Position',[-2.6 4.2],...
  'FontSize',11, 'Color', 'white')

  sigma_str = '$$\Sigma = \left[ \begin{array}{cc} 0.8 & 0.7 \\ 0.7 & 1.3\end{array} \right]$$';
  text('Interpreter','latex',...
  'String',sigma_str,...
  'Position',[-2.6 3],...
  'FontSize', 11, 'Color', 'white')

  % Add title and axis labels
  set(gca,'Title',text('String','Bivariate Normal Density, Full Covariance',...
  'FontAngle', 'italic', 'FontWeight', 'bold'),...
  'xlabel',text('String', '$\mathbf{X}$', 'Interpreter', 'latex'),...
  'ylabel',text('String', '$\mathbf{Y}$', 'Interpreter', 'latex'))

  % Adjust the view angle
  view(0, 90);
  colormap('hot')

And here is the code that produced the 3-$D$ mesh-grid view of the density:

  %plot the surface in 3-d (no need to recompute meshgrid but done anyway)
  figure(2)
  [X,Y] = meshgrid(-3:0.1:4,-2:0.1:5);
  const = (1/sqrt(2*pi))^2;
  const = const/sqrt(det(sigma));
  temp = [X(:)-mu(1) Y(:)-mu(2)];
  pdf = const*exp(-0.5*diag(temp*inv(sigma)*temp'));
  pdf = reshape(pdf,size(X));

  mesh(X, Y, pdf);

  % Add title and axis labels
  set(gca,'Title',text('String','A 3-D View of the Bivariate Normal Density',...
  'FontAngle', 'Italic', 'FontWeight', 'bold'),...
  'xlabel',text('String', '$\mathbf{X}$', 'Interpreter', 'latex'),...
  'ylabel',text('String', '$\mathbf{Y}$', 'Interpreter', 'latex'),...
  'zlabel',text('String', 'density', 'FontAngle', 'Italic', 'FontWeight', 'bold'))
  view(-10, 50)
  colormap('winter')

Copy and paste the code into your MATLAB environment and experiment by changing the parameters of the distribution and behold the different shapes.

BONUS: Exercise: Alter the parameters in the code so that the output is an ellipsoid that is slanted from the top left to the bottom right of the figure. What is the covariance matrix in this case?