1
$\begingroup$

Mathematica noob here, and I'm having some issues plotting 90% confidence ellipses of a small-variance 2D Gaussian PDF. I have it working perfectly fine for larger covariance matrices, but absolutely no luck with these tiny ones! Here's my code:

Sigma={{0.0093041, -0.00126552}, {-0.00126552, 0.00020695}};
xv={X,Y};
x0={66.2782, 0.791879};
pdf = 1/Sqrt[(2 \[Pi])^Length[xv] Det[Sigma]]Exp[-(xv - x0).Inverse[Sigma].(xv - x0)/2]
ContourPlot[pdf, {X, 66, 66.6}, {Y, .76, .83}, PlotRange -> All, PlotPoints -> 100]

Which gives me a standard multivariate Gaussian, albiet with contours at densities of 50-250. Which I guess is fine assuming such a small range. Also, just as a quick check, the integral does indeed equal 1:

NIntegrate[pdf, {X, 60, 70}, {Y, .5, .9}](*=1*)

Now, trying to compute the 90% confidence regions.. Following the solution here: https://mathematica.stackexchange.com/a/20481/58658, I do the following:

f[x_, y_, t_: 0] := With[{z = pdf /. X -> x /. Y -> y}, z Boole[z >= t]];
r = FindRoot[NIntegrate[f[x, y, t], {x, 60, 70}, {y, .7, .9}, AccuracyGoal -> 3,
  PrecisionGoal -> 6] - 0.90, {t, 0, 2}, Method -> "Brent", AccuracyGoal -> 3, PrecisionGoal -> 6];

Which results in numerous errors "The integrand has evaluated to non-numerical values for all sampling points in the region..", or non-sensical answers as I play around with the integration ranges. I don't really see how to make this work, seeing as all the pdf densities evaluate to be much greater than 1 in such a small range!! Does anyone have any experience with issues like this? Or any clever roundabout ways of doing this?

Thanks so much in advance!

$\endgroup$
  • $\begingroup$ The first contour plot gives you "densities" not "probabilities" at 50, 100, 150, 200, and 250. $\endgroup$ – JimB Apr 8 at 1:22
  • $\begingroup$ Oops, my bad! I've edited the question. Thanks! $\endgroup$ – zack Apr 8 at 1:29
  • $\begingroup$ From the answer below you'll see that it really has nothing to do with the size of the variances. $\endgroup$ – JimB Apr 8 at 4:50
3
$\begingroup$

One can do it "exactly". A $100(1-\alpha)%$ confidence ellipsoid (for a multivariate normal) is given by

$$(X-\mu)' \Sigma^{-1} (X-\mu)=\chi^2_k(1-\alpha)$$

where $k$ is the dimension of the multivariate normal and $\chi^2_k(1-\alpha)$ is the $\chi^2$ value with $k$ degrees of freedom that has $1-\alpha$ of the probability to the left of it. So for a bivariate normal the 90% contour you requested is

Sigma={{0.0093041, -0.00126552}, {-0.00126552, 0.00020695}};
xv={x,y};
x0={66.2782, 0.791879};
pdf = 1/Sqrt[(2 π)^Length[xv] Det[Sigma]]Exp[-(xv - x0).Inverse[Sigma].(xv - x0)/2]
c90 = (2 π)^(-1) Det[Sigma]^(-1/2) Exp[-InverseCDF[ChiSquareDistribution[2], 0.9]/2]
ContourPlot[pdf, {x, 66, 66.6}, {y, .75, .83}, 
  PlotPoints -> 100, Contours -> {c90}, ContourShading -> None]

Bivariate normal contour

As a check...

NIntegrate[
 Boole[(xv - x0).Inverse[Sigma].(xv - x0) <= InverseCDF[ChiSquareDistribution[2], 0.90]] pdf,
 {x, 66, 66.6}, {y, 0.75, 0.83}]
(* 0.9 *)
$\endgroup$
  • $\begingroup$ This method is so much easier (and much faster)!! Thank you so much @JimB. $\endgroup$ – zack Apr 8 at 12:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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