9

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?

J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
K-1
  • 581
  • 1
  • 3
  • 12

3 Answers3

12

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.

Figure

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.

whuber
  • 20,544
  • 2
  • 59
  • 111
  • Thank you for suggesting a more flexible way of plotting this ellipsoid. The only part I cannot fully understand is the "Level t". What does it do exactly? – K-1 Mar 18 '13 at 01:31
  • @K-1 There are infinitely many confidence ellipses ("error ellipsoids"), each determined by its corresponding "confidence level," a number between $0$ and $1$. That's what t is based on: you can see its calculation at the beginning of the Block. – whuber Apr 08 '13 at 15:49
  • 2
    NMaxValue[] would allow for a more compact implementation of limit[]. – J. M.'s missing motivation Apr 11 '13 at 10:35
  • If you are interested in only two variables you should not slice the ellipsoid!! Wrong! You have to do marginalization, which is in this case even simpler. Like I explained here you have to modify the covariance matrix and perform this brilliant procedure again for reduced number of variables. – Vladimir Apr 19 '13 at 05:56
  • @Vladimir It all depends on what you want to know. I haven't advocated any particular interpretation of the slices, so there's no possible way in which slicing can be considered "wrong": you must have in mind some incorrect way in which you are thinking of misinterpreting such an analysis (and I won't even pretend to guess what that might be). Studying the likelihood function is useful and one way to learn about it is to look at slices. Slicing the covariance level surfaces is just a second-order approximation to slicing the likelihood itself. – whuber Apr 19 '13 at 11:55
  • Certainly, whuber. The comment was not for you, but for beginners who have that common mistake in interpretation. As the OP didn't really know all parts of procedure when he posted the question, he also couldn't specify exactly what he needed. The comment for you, whuber, would be to give interpretation on what you are doing as you are the one who has the knowledge :) And therefore, avoid wrong interpretations by beginners. – Vladimir Apr 19 '13 at 20:50
  • 1
    @Vladimir I appreciate your nuance, but something's missing: you still haven't even stated what your interpretation is! That is, for what purposes would one "do marginalization"? Until you're explicit and forthcoming about that, I don't see how your remark could help beginners. It seems to me that the only clear message it conveys is "Bad! Don't pay attention to this answer! Read my question instead!" – whuber Apr 19 '13 at 20:54
  • 1
    hahaha No, no, sorry, that's not what I wanted to say! I was more explicit to drag attention, only because comments are not often read. Honestly, I myself learned a lot from this answer and it is actually the answers for my question, too! (which is why I called it brilliant, as it is for any dimension). And now: if you are not interested in one variable and you want to see confidence regions for the rest of variables, you can either put some fixed value to that one variable and get new correlation matrix (by fitting again), or integrate probability function over that variable, eg. marginalize. – Vladimir Apr 19 '13 at 21:42
8

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

rm -rf
  • 88,781
  • 21
  • 293
  • 472
5

Since Ellipsoid[] is now built-in, and has a different syntax from the version that was once in MultivariateStatistics`​, let me present the way to render the confidence ellipsoid corresponding to a given covariance matrix:

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

With[{α = 0.05}, (* level of significance *)
     ell = Ellipsoid[{0, 0, 0}, InverseCDF[ChiSquareDistribution[3], 1 - α] c]];

Graphics3D[ell, Axes -> True]

confidence ellipsoid

If you want to render the ellipsoid's contours along its eigendirections, you can do this:

vecs = Eigenvectors[c];

RegionPlot3D[ell, Axes -> True, AxesStyle -> Opacity[1, Black], BaseStyle -> Transparent,
             Mesh -> 10, MeshFunctions -> Evaluate[Map[Function, vecs.{#1, #2, #3}]], 
             MeshStyle -> Opacity[1, Black], PlotStyle -> Opacity[1, White], 
             TicksStyle -> Opacity[1, Black]]

eigencontours

where some of the option settings (e.g. BaseStyle -> Transparent) was needed to work around a bug in RegionPlot3D[] involving meshline rendering.

J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574