0
$\begingroup$

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

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.

Improve this question
$\endgroup$
  • $\begingroup$ you should change your implicit region definition from equality to < 1 $\endgroup$ – ubpdqn Aug 2 '14 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 '14 at 12:01
1
$\begingroup$

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_] := 
 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}]

enter image description here

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

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.