How to plot a covariance error ellipse

A typical way to visualize two-dimensional gaussian distributed data is plotting a confidence ellipse. Lets assume we have data \(D\sim\mathcal{N}(\mu, \Sigma)\) and want to plot an ellipse representing the confidence \(p\) by calculating the radii of the ellipse, its center and rotation. The following plot shows randomly drawn data and the ellipses for \(p\in\{0.9, 0.95, 0.99\}\):

Derivation

If the data is uncorrelated and therefore has zero covariance, the ellipse is not rotated and axis aligned. The radii of the ellipse in both directions are then the variances.

Geometrically, a not rotated ellipse at point \((0, 0)\) and radii \(r_x\) and \(r_y\) for the x- and y-direction is described by

\[\left(\frac{x}{r_x}\right)^2 + \left(\frac{y}{r_y}\right)^2 = 1\]

In our case, the radius in each direction is the standard deviation \(\sigma_x\) and \(\sigma_y\) parametrized by a scale factor \(s\), known as the Mahalanobis radius of the ellipsoid:

\[\left(\frac{x}{\sigma_x}\right)^2 + \left(\frac{y}{\sigma_y}\right)^2 = s\]

The goal must be to determine the scale \(s\) such that confidence \(p\) is met. Since the data is multivariate gaussian distributed, the left hand side of the equation is the sum of squares of gaussian distributed samples, which follows a \(\chi^2\) distribution. A \(\chi^2\) distribution is defined by the degrees of freedom and since we have two dimensions, the number of degrees of freedom is also two. We now want to know the probability that the sum and therefore \(s\) has a certain value under a \(\chi^2\) distribution. Since we have a confidence interval \(p\) we formulate the problem to determine the probability having \(s\) smaller than a certain value \(z\):

\[P(s < z) = p\]

We can solve this equation using a \(\chi^2\) table of a statistical text book, using Matlab function s=chi2inv(p, 2) or simply:

\[s = -2\log(1-p)\]

For the example of \(p=0.9\) we get \(s=4.6057\), for \(p=0.95\) we get \(s=5.99146\) and for \(p=0.99\) we get \(s=9.21034\). The ellipse can then be drawn with radii \(\sigma_x\sqrt{s}\) and \(\sigma_y\sqrt{s}\).

Generaliztion for given covariance

In the general case, covariances \(\sigma_{xy}\) and \(\sigma_{yx}\) are not zero and therefore the ellipse-coordinate system is not axis-aligned. However, the previously derived approach is still usable, but instead of using the variance as a spread indicator, we use the eigenvalues of the covariance matrix \(\Sigma=\left(\begin{array}{cc}\sigma_x^2&\sigma_{xy}\\\sigma_{yx}&\sigma_y^2\end{array}\right)\). The eigenvalues represent the spread in the direction of the eigenvectors, which are the variances under a rotated coordinate system. By definition a covariance matrix is positive definite therefore all eigenvalues are positive and can be seen as a linear transformation to the data.

The actual radii of the ellipse are \(\sqrt{\lambda_1}\) and \(\sqrt{\lambda_2}\) for the two eigenvalues \({\lambda_1}\) and \({\lambda_2}\) of the scaled covariance matrix \(s\cdot\Sigma\). It can easily be shown that \(\text{eigenvalue}(s\cdot\Sigma)=s\cdot\text{eigenvalue}(\Sigma)\), which simplifies the scaling of the eigenvalues afterwards.

Matlab Implementation

Having all these information together, an implementation of the procedure in Matlab can be made with a linear combination like so:

function plotErrorEllipse(mu, Sigma, p)

    if nargin == 2
      p = 0.95; 
    end

    s = -2 * log(1 - p);

    [V, D] = eig(Sigma * s);

    t = linspace(0, 2 * pi);
    a = (V * sqrt(D)) * [cos(t(:))'; sin(t(:))'];

    plot(a(1, :) + mu(1), a(2, :) + mu(2));

JavaScript Implementation

For JavaScript the case is a little more complicated since we don't have access to linear algebra functions natively and must calculate the eigenvalues ourselves. Furthermore we don't want to draw the ellipse with many small lines, but use the canvas ellipse function or the SVG ellipse tag, where we need to calculate the parameters explicitly.

For a covariance matrix of the form \(\Sigma=\left(\begin{array}{cc}a&b\\b&d\end{array}\right)\), we can compute the eigenvectors (\(V\)) and eigenvalues (\(D\)) like this:

\[V =\left(\begin{array}{cc} -\frac{\sqrt{(a - d)^2 + 4 b^2} - a + d}{2 c} & \frac{\sqrt{(a - d)^2 + 4 b^2} + a - d}{2 c}\\ 1& 1 \end{array}\right)\]

\[D =\left(\begin{array}{cc} \frac{1}{2} (-\sqrt{(a - d)^2 + 4 b^2} + a + d)& 0\\ 0& \frac{1}{2} (\sqrt{(a - d)^2 + 4 b^2} + a + d) \end{array}\right)\]

After normalizing the column vectors in V, we choose the eigenvector with the larger eigenvalue and calculate its angle to the global x-axis. This is the rotation of the ellipse. All that is needed then is calculating the radii of the ellipse.

function plotErrorEllipse(ctx, mu, Sigma, p) {

  p = p || 0.95;

  var s = -2 * Math.log(1 - p);

  var a = Sigma[0][0];
  var b = Sigma[0][1];
  var c = Sigma[1][0];
  var d = Sigma[1][1];

  var tmp = Math.sqrt((a - d) * (a - d) + 4 * b * c);
  var V = [
    [-(tmp - a + d) / (2 * c), (tmp + a - d) / (2 * c)],
    [1, 1]
  ];
  var sqrtD = [
    Math.sqrt(s * (a + d - tmp) / 2),
    Math.sqrt(s * (a + d + tmp) / 2)
  ];

  var norm1 = Math.hypot(V[0][0], 1);
  var norm2 = Math.hypot(V[0][1], 1);
  V[0][0] /= norm1;
  V[1][0] /= norm1;
  V[0][1] /= norm2;
  V[1][1] /= norm2;

  var ndx = sqrtD[0] < sqrtD[1] ? 1 : 0;

  var x1 = mu[0] + V[0][ndx] * sqrtD[ndx];
  var y1 = mu[1] + V[1][ndx] * sqrtD[ndx];

  var x2 = mu[0] + V[0][1 - ndx] * sqrtD[1 - ndx];
  var y2 = mu[1] + V[1][1 - ndx] * sqrtD[1 - ndx];

  ctx.ellipse(
          mu[0], mu[1],
          Math.hypot(x1 - mu[0], y1 - mu[1]),
          Math.hypot(x2 - mu[0], y2 - mu[1]),
          Math.atan2(y1 - mu[1], x1 - mu[0]),
          0, Math.PI * 2,
          false);
}