Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question

migrated from stats.stackexchange.com Mar 15 '13 at 14:05

This question came from our site for people interested in statistics, machine learning, data analysis, data mining, and data visualization.

2 Answers 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.

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.

share|improve this answer
    
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 1: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 8 '13 at 15:49
2  
NMaxValue[] would allow for a more compact implementation of limit[]. –  J. M. Apr 11 '13 at 10:35
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

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

share|improve this answer
    
Thank you rm-rf for your answer. It was quite helpful! –  K-1 Mar 18 '13 at 1:30

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.