# Plotting 90% confidence ellipses for a gaussian pdf with small variances

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?

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

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]


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 *)

• This method is so much easier (and much faster)!! Thank you so much @JimB. – zack Apr 8 '19 at 12:58