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}}]

enter image description here

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}]


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.

  • $\begingroup$ you should change your implicit region definition from equality to < 1 $\endgroup$
    – ubpdqn
    Aug 2, 2014 at 11:39
  • $\begingroup$ 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. $\endgroup$
    – anthr
    Aug 2, 2014 at 12:01

1 Answer 1


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[2]];
reg3[a_, b_] := 
  x^2/a^2 + y^2/b^2 < 1 && 0 <= z < PDF[bn, {x, y}], {x, y, z}]

Visualizing to illustrate:

   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}]

enter image description here

  • $\begingroup$ Hugely useful - thank you so much for your help on these two questions! $\endgroup$
    – anthr
    Aug 3, 2014 at 12:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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