# Using ImplicitRegion to define an ellipse around Multinormal distribution for integration

I have both a 2D 'MultinormalDistribution', and also a single xy point, and I would like to be able to calculate the probability of this point (given the multinormal distribution) and also plot an ellipse with the same probability as this point around the mean of the distribution.

So far I have defined the Multinormal (xy) distribution like:

bvn = MultinormalDistribution[{-3.89764, 1.29137}, {{0.08369444426^2, 0}, {0, 0.0175124^2}}]


I also have the fixed point in this xy plane, at -3.98529, 1.30482. Previously (thanks to help from users ubpdqn and Karsten 7 in a previous question) I approached this by essentially doing two-tailed Gaussian integrals in both of the normal distributions that make up the multinormal distribution, which gave me the probability inside a square, represented visually with:

Plot3D[PDF[bvn, {x, y}], {x, -4.5, -3.5}, {y, 1.2, 1.4}, PlotRange -> All,    MeshFunctions -> {#1 &, #2 &}, Mesh -> {{xl = -3.98529, xu = -3.80999}, {yl = 1.27792, yu = 1.30482}}] I realize now that this square isn't exactly what I want, and instead I want to define an ellipse of constant probability around the mean, consistent with this xy point.

ubpdqn suggested several approaches to performing the integration in the square, including:

NIntegrate[PDF[bvn, {x, y}], {x, xl, xu}, {y, yl, yu}]


and:

reg = ImplicitRegion[xl < x < xu && yl < y < yu, {x, y}]
NIntegrate[PDF[bvn, {x, y}], {x, y} \[Element] reg]


I have attempted to modifiy this implicit region such that it becomes an ellipse centered on the mean of the multinormal distribution (-3.89764, 1.29137), and bounded by the xy point like this:

xl = -3.98529;
yl = 1.27792;
xc = -3.89764;
yc = 1.29137;
reg2 = ImplicitRegion[((x - xc)^2/xl^2) + ((y - yc)^2/yl^2) == 1, {x,y}];
NIntegrate[PDF[bvn2, {x, y}], {x, y} \[Element] reg2]


Although this isn't yielding a sensible answer so it is clear to me that I'm doing something wrong. If anyone could suggest a modification to my approach I would be very grateful.

• you should change your implicit region definition from equality to < 1 Aug 2 '14 at 11:39
• Thanks again for your help. I guess I still have an error somewhere because changing to "ImplicitRegion[(x - xc)^2/(xd^2) + (y - yc)^2/(yd^2) < 1, {x, y}]" results in the integration having an imaginary component, although when I plot this region with ContourPlot it looks correct. Aug 2 '14 at 12:01

I post this for illustrative purposes. For thie particular distribution, there are two issues (i) the complex number can be handled with Chop (ii) the very small variance relative to the ellipical region.

In the following I use a standard binormal distribution and use NIntegrate on region:

bn = MultinormalDistribution[{0, 0}, IdentityMatrix];
reg3[a_, b_] :=
ImplicitRegion[
x^2/a^2 + y^2/b^2 < 1 && 0 <= z < PDF[bn, {x, y}], {x, y, z}]


Visualizing to illustrate:

Manipulate[Column[{
Show[Plot3D[PDF[bn, {x, y}], {x, -3, 3}, {y, -3, 3}, Mesh -> False,
PlotStyle -> Opacity[0.4], PerformanceGoal -> "Quality",
PlotRange -> All, ImageSize -> 400],
RegionPlot3D[reg3[a, b], PlotStyle -> {Blue, Opacity[0.6]}]],
Chop[NIntegrate[PDF[bn, {x, y}], {x, y} \[Element] reg[a, b]]]
}, Alignment -> Center], {a, 1, 3}, {b, 1, 3}] • Hugely useful - thank you so much for your help on these two questions! Aug 3 '14 at 12:38