How to use a 3×3 covariance matrix to plot an error ellipsoid?

Question

I have a 3×3 error covariance in Mathematica, but I don't know how to use it for plotting the error ellipsoid. It would be great if you can show me how I can do that for the below covariance matrix:

CovMat= {{88.5333, -33.6, -5.33333}, 
         {-33.6, 15.4424, 2.66667}, 
         {-5.33333, 2.66667, 0.484848}}

eigenvalues= {0.0098, 0.4046, 104.7}

eigenvectors= {{0.93, 0.36, -0.03}, {-0.36, 0.9, -0.23}, {-0.06, 0.23, 0.97}}

And as the last question, how can I project this ellipsoid onto 2D planes?

Answer 1

Ellipsoids can be easily plotted using the MultivariateStatistics` package. The eigenvalues of your covariance matrix denote the lengths of the axes and the eigenvectors, their orientation. Here's how it can be done for your CovMat:

CovMat = {{88.5333, -33.6, -5.33333}, 
          {-33.6, 15.4424, 2.66667}, 
          {-5.33333, 2.66667, 0.484848}};

e = Ellipsoid[{0, 0, 0}, Sequence @@ Eigensystem[CovMat]];
Graphics3D[e, BoxRatios -> 1, Axes -> True, Boxed -> False, AxesLabel -> {"x", "y", "z"}]

You can get the different projections by altering the ViewPoint:

Graphics3D[e, BoxRatios -> 1, Axes -> True, Boxed -> False, 
    AxesLabel -> {"x", "y", "z"}, ViewPoint -> #, 
    ImageSize -> 200] & /@ Permutations[{0, 0, Infinity}, {3}] // Row[#, Spacer[10]] &

If you want better control over the appearance of your ellipsoid, you should do the plotting yourself with ParametricPlot3D (which is what Ellipsoid does behind the scenes). For example:

ParametricPlot3D[
   Transpose[#2].(# {Cos[t] Cos[u], Sin[t] Cos[u], Sin[u]}), {t, 0, 2 π}, {u, -(π/2), π/2}, 
    Boxed -> False, BoxRatios -> 1, AxesLabel -> {"x", "y", "z"}, 
   ColorFunction -> "AvocadoColors", Mesh -> False] & @@ Eigensystem@CovMat

Answer 2

For full control over the plot and the analysis, it is useful to know how to do the calculations yourself. They include:

• Finding the contour level corresponding to the desired confidence (without this the result scarcely can be called a "confidence" ellipse!);

• Determining the limits (and aspect ratios) of the plot;

• Establishing a useful mesh on the ellipsoid (showing contours along the eigendirections).

These steps should be apparent in the lines of the following code, which takes for input a number a (the confidence will be $1-a$) and the covariance matrix, here called c:

limit[ci_, n_, t_] := Abs[n.{x, y, z}] /. 
  Last[NMaximize[{(n.{x, y, z})^2, {x, y, z}.ci.{x, y, z} <= t}, {x, y, z}]];
Block[{t = InverseCDF[ChiSquareDistribution[3], 1 - 0.05], ci = Inverse[c], x0, y0, z0, mf},
 {x0, y0, z0} = 1.05 limit[ci, #, t] & /@ IdentityMatrix[3];
 mf = Function[{x, y, z}, #.{x, y, z}] & /@ Eigenvectors[c];
 ContourPlot3D[{x, y, z}.ci.{x, y, z} == t, {x, -x0, x0}, {y, -y0, y0}, {z, -z0, z0}, 
  AxesLabel -> {x, y, z}, 
  ContourStyle -> Opacity[0.8],
  MeshFunctions -> mf , MeshStyle -> Opacity[0.2], 
  BoxRatios -> {x0, y0, z0}]]

limit finds the extreme absolute values of any linear form (given by a three-vector n) along the ellipse specified by matrix ci at level t. You can observe its use within the subsequent Block where, by applying it to the three vectors $(1,0,0)$, $(0,1,0)$, and $(0,0,1)$ (the rows of IdentityMatrix), we obtain the extreme values of $x$, $y$, and $z$ in the plot (and expand them by five percent to give a small margin).

mf computes the distances along the eigenvectors.

InverseCDF computes the proper contour level for the desired confidence.

Note that the confidence ellipsoid is a contour of the inverse of the covariance matrix.

---

In another answer, @rm-rf has given some expedient ways to plot projections. You can also use the technique given above to draw slices through the ellipsoid: simply fix one of $x$, $y$, or $z$ (perhaps by making it controllable by Manipulate) and invoke ContourPlot instead of ContourPlot3D, changing the arguments in obvious ways. This would be a nice way to obtain conditional confidence ellipses.