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.

Please consider the following :

covMatrice = {{34.925, -10.21}, {-10.21, 22.462}};

COG = Mean@(Nfixations["d"][[1]])[[All, 1]];

fixations = {{19.4688, 17.4281}, {18.0563, 21.7156}, {13.0219, 24.7219}, 
         {22.9594,25.5219}, {28.5406, 24.6719}, {27.0688, 17.1656}, 
         {27.6781,16.325}, {28.9281, 10.7719}, {16.025, 13.6625}, 
         {19.1313, 17.1094}};

      eigVec = Eigenvectors[covMatrice],
      eigVal = Eigenvalues[covMatrice]},
                White, Rectangle @@ frmXY, 
                Black, Disk[#, .5] & /@ fixations, 
                Red, Line[(# + COG) & /@
                           {0, 0},

enter image description here

How could I draw an ellipse representing on the EigenSystem given that Neither Disk or Circle enable to implement an "orientation" ?

A rough example of my desired output as drawn using PPT :

enter image description here

The documentation suggests using rotate to get all possible ellipses. (See the answers below for code!) –  Simon Jan 27 '12 at 3:25
@Simon Rotate works fine if you know the angle. As the OP only has the new coordinates, GeometricTransformation is the better choice. –  rcollyer Jan 27 '12 at 3:27
If I may be indulged some numerical advice: it might not always be a good idea to form a covariance matrix from your data for the express purpose of seeing ellipses (or ellipsoids) for that matter. You should look into the singular value decomposition, for instance. –  J. M. Jan 27 '12 at 17:36
Just careful, it's the roots of the eigenvalues of the covariance matrix what you want to be using as semi-axis lengths –  Rojo Nov 19 '12 at 7:07

4 Answers 4

up vote 7 down vote accepted

You can use Rotate to draw the ellipse too. Note that knowing the eigenvectors is the same as knowing the orientation of the ellipse, so there is no necessity to favour GeometricTransformation over Rotate. So, to orient your ellipse along the first eigenvector (corresponding to the largest eigenvalue), it is as simple as:

Graphics[Rotate[Disk[meanVec, eigVals], ArcTan @@ eigVecs[[1]]]]

enter image description here

where meanVec is the mean (here, I've taken it to be {0,0})

Going by your comment under rcollyer's answer, here's a example reproducing your desired figure with simple shifting (change the center of the disk) and rotation of the disk. This approach will be simpler to follow (as the transformations are spelled out), if you do not understand what GeometricTransformation does. Modifying rcollyer's module,

Module[{mat = #, avg = Mean@fixations, eigVals, eigVecs},
   {eigVals, eigVecs} = Eigensystem@mat;
   Graphics[{{Black, Disk[#, .5] & /@ fixations}, 
      {Directive[Opacity[0.1], Red, EdgeForm[Gray]], 
          Rotate[Disk[avg, eigVals/5], ArcTan @@ eigVecs[[1]]]}

enter image description here

I figured I needed to give you an upvote since you clearly answer the question. :P –  rcollyer Jan 27 '12 at 20:38

The simplest approach is to rotate an ellipse aligned with the coordinate axes to axes specified by the eigenvectors. Using your matrix,

{evals, evects} = Eigensystem@covMatrice
=> {{40.6549, 16.7321}, {{-0.872057, 0.489405}, {-0.489405, -0.872057}}}

Graphics[ GeometricTransformation[ Disk[{0,0}, evals], evecs\[Transpose] ] ]


rotated ellipse

Edit: From the comments, the previous code was incomplete in that the ellipse needed to be centered over the other graphics objects. GeometricTransformation functions by taking as its second argument a TransformationFunction and applying it to the Graphics objects passed in via its first argument. An alternate form of GeometricTransformation accepts a matrix representing the transformation, instead, and that is what the above code uses. In essence, every point, r, in the object is turned into m.r, where m is the matrix. The third form accepts a matrix and a translation vector, v, transforming each point into m.r + v. The following code uses that form to center the ellipse at the mean position of the fixations:

Module[{mat = #, avg = Mean@fixations, g, evals, evecs},
  {evals, evecs} = Eigensystem@mat;
  g = GeometricTransformation[#, {evecs\[Transpose], avg}] &;
    {Black, Disk[#, .5] & /@ fixations},
    {Directive[Opacity[0.1], Red,  EdgeForm[Gray]], g@Disk[{0, 0}, evals/5]},
    {Red, Line[(# + avg) & /@ Riffle[evecs evals/5, {{0, 0}}] ]}
] & @ covMatrice

which produces

rotated ellipse translated over data points

I consider this method superior to Rotate in that no other calculations are needed.

thank you it takes the good shape but I can`t place it correctly with my other Graphics objects there is a systematic shift toward the {0,0}. I don`t understand GeometricTransformation and can`t figure out how to adjust the coordinates of the resulting graphics object. –  500 Jan 27 '12 at 5:03
@500, I've added the correction. Sorry about that. You wanted to use the third form of GeometricTransformation not the second, which is what I was using. –  rcollyer Jan 27 '12 at 14:48

If you know the angle, try using Rotate, simple example:

Graphics[{Rotate[Disk[{0, 0}, {1, 2}], \[Pi]/4]}]

(note: the two argument form of Rotate rotates around the center of the primitive)

Here is a more complete example with a Manipulate:

   EdgeForm[{AbsoluteThickness[3], Gray}],
   Rotate[Disk[{cx, cy}, {rx, ry}], a, {cx, cy}],
   {Red, Rotate[
     Line[{{cx, cy}, {cx + rx, cy}}], a, {cx, cy}]},
   {Red, Rotate[
     Line[{{cx, cy}, {cx, cy + ry}}], a, {cx, cy}]},
   {Blue, Line[{{0, 0}, {cx, cy}}]}
   }, PlotRange -> 7, Axes -> True, GridLines -> Automatic],
 {{a, 2}, 0, 2 \[Pi]},
 {{rx, 3}, .5, 3},
 {{ry, 1}, .5, 3},
 {{cx, 3}, -3, 3},
 {{cy, 3}, -3, 3}

Mathematica graphics

(If you don't know the angle, @rcollyer's answer seems the more useful one).

Thank You very useful demo. –  500 Jan 27 '12 at 4:30

Another approach to solving your problem would be to use the MultivariateStatistics package, which includes the EllipsoidQuantile function. Here's an example:


fixations = {{19.4688, 17.4281}, {18.0563, 21.7156}, {13.0219, 
24.7219}, {22.9594, 25.5219}, {28.5406, 24.6719}, {27.0688, 
17.1656}, {27.6781, 16.325}, {28.9281, 10.7719}, {16.025, 
13.6625}, {19.1313, 17.1094}};

e50 = EllipsoidQuantile[fixations, 0.5]; (* use the median ellipse for this example *)

Now extract the ellipse axes and plot it:

major1 = e50[[1]] + e50[[3]]\[Transpose].{e50[[2, 1]], 0}; (* construct the major and minor axes for plotting *)
major2 = e50[[1]] - e50[[3]]\[Transpose].{e50[[2, 1]], 0};
minor1 = e50[[1]] + e50[[3]]\[Transpose].{0, e50[[2, 2]]};
minor2 = e50[[1]] - e50[[3]]\[Transpose].{0, e50[[2, 2]]};

 ListPlot[fixations, AspectRatio -> Automatic, 
  PlotRange -> {All, {10, 30}}],
 Graphics[{Red, Thick, e50}],
 Graphics[{Red, Thick, Line[{major1, major2}]}],
 Graphics[{Red, Thick, Line[{minor1, minor2}]}]

enter image description here

Your example data doesn't have many points which makes it difficult (for me at least) to visualize the answer. Assuming the data are Normally distributed, you can generate some extra points this way:

mean = e50[[1]];
covMatrice = {{34.925, -10.21}, {-10.21, 22.462}};
\[ScriptCapitalD] = MultinormalDistribution[mean, covMatrice];
data = RandomVariate[\[ScriptCapitalD], 1000];

enter image description here

e50 = EllipsoidQuantile[data, 0.5];
major1 = e50[[1]] + e50[[3]]\[Transpose].{e50[[2, 1]], 0};
major2 = e50[[1]] - e50[[3]]\[Transpose].{e50[[2, 1]], 0};
minor1 = e50[[1]] + e50[[3]]\[Transpose].{0, e50[[2, 2]]};
minor2 = e50[[1]] - e50[[3]]\[Transpose].{0, e50[[2, 2]]};

 ListPlot[data, AspectRatio -> Automatic,PlotRange -> {All, {10, 30}}],
 Graphics[{Red, Thick, e50}],
 Graphics[{Red, Thick, Line[{major1, major2}]}],
 Graphics[{Red, Thick, Line[{minor1, minor2}]}]

enter image description here


This site is currently not accepting new answers.

Not the answer you're looking for? Browse other questions tagged .