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!

Improve this question
$\endgroup$
3
  • $\begingroup$ The first contour plot gives you "densities" not "probabilities" at 50, 100, 150, 200, and 250. $\endgroup$ Commented Apr 8, 2019 at 1:22
  • $\begingroup$ Oops, my bad! I've edited the question. Thanks! $\endgroup$ Commented Apr 8, 2019 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$ Commented Apr 8, 2019 at 4:50

1 Answer 1

Reset to default
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 *)
Improve this answer
$\endgroup$
1
  • $\begingroup$ This method is so much easier (and much faster)!! Thank you so much @JimB. $\endgroup$ Commented Apr 8, 2019 at 12:58

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.