Message Boards Message Boards

GROUPS:

2D likelihood. Calculate ellipses that contains a 95% of the area

Posted 1 year ago
1544 Views
|
7 Replies
|
0 Total Likes
|

Dear All,

I have a function of two variable (function[x,y] in the notebook attached) with which I build a 2D Likelihood from a Poissonian probability distribution function.

Given the structure of the function, the likelihood has isocontours which are ellipses in the x-y plane, as per figure. What I want to do find is the following. I would like to find the ellipses for which the integral of the LKLH[x,y] within it is 95%, e.g., of the total volume.

Following the example here (http://reference.wolfram.com/language/tutorial/IntegralsOverRegions.html) I tried to do something like

FindRoot[NIntegrate[
   If[LKLH[x, y] < a, 1, 0], {x, -500, 500}, {y, -500, 500}] == 
  0.95, {a, 1}]

but this doesn't work. Any idea on how one can work out this? It believe is pretty trivial, since it occurs anytime one has a 2D probability distribution function...

Thanks a lot to anyone

Attachments:
7 Replies

try testing your LKLH function.

I'm not sure what you mean with this suggestion. Could you please clarify? Thanks :-)

try evaluating it with numerical arguments

In[3]:= Poisson[nth_, nobs_] := nth^nobs Exp[-nth]/nobs!

In[4]:= function[x_, 
  y_] := (1 + 1.5 10^-3 x + 2.5 10^-3 x^2 + 1.7 10^-3 y + 
   1.8 10^-3 y^2 + 1.8 10^-3 x y)

In[8]:= LKLH[x_?NumericQ, y_?NumericQ] := 
 Poisson[function[x, y], function[0, 0]]/
 NIntegrate[
  Poisson[function[x, y], function[0, 0]], {x, -500, 500}, {y, -500, 
   500}]

In[14]:= LKLH[250, 250]

During evaluation of In[14]:= NIntegrate::dupv: Duplicate variable 250 found in NIntegrate[Poisson[function[250,250],function[0,0]],{250,-500,500},{250,-500,500}].

Out[14]= 1.6856*10^-164/NIntegrate[
 Poisson[function[250, 250], function[0, 0]], {250, -500, 
  500}, {250, -500, 500}]

Notice first that the problem can be simplified, since function[0, 0]=1, and one can compute K=NIntegrate[LKLH[x, y], {x, -100, 100}, {y, -100, 100}]. Now, we can define, with the same function:

Truc[nth_] := nth Exp[-nth];
LKLH[x_, y_] := Truc[function[x, y]]/K;

Since "function" is the equation of a conic, we can study it.

Pol = Expand@FullSimplify@(function[x, y] - L)
M = CoefficientList[%, {x, y}];
règles = CoefficientRules[Pol, {x, y}];

Put it in canonical form:

{a, b, c, d, e, 
   f} = {{2, 0}, {1, 1}, {0, 2}, {1, 0}, {0, 1}, {0, 0}} /. règles;
b^2 - 4 a c

Aq = {{a, b/2, d/2}, {b/2, c, e/2}, {d/2, d/2, f}};
MatrixForm@Aq
 Det[Aq] // Factor

Since the determinant is negative for L>= 1, its is the equation of an ellipse under this condition. Suppose now L>1; The solution of Pol==0 is the contour : function[x,y]==L. I think you can now write formally:

\[Integral]Truc[
    function[x, y]] \[DifferentialD]x \[DifferentialD]y = \!\(
\*SubsuperscriptBox[\(\[Integral]\), \(0\), \(+\[Infinity]\)]\(Truc[
    L] \[DifferentialD]L\)\)

Then, if it's true, you just have to compute:

Q = Quantile[GammaDistribution[2, 1], 0.95]

since Truc is indeed this Gamma distribution, Then, determine the ellipse, and plot it.

Ellips = Solve[Expand@FullSimplify@(function[x, y] - Q) == 0, {x, y}];
Print@Show[{ParametricPlot[{x, y /. Ellips[[1]]}, {x, -50, 50}, 
     PlotStyle -> Red], 
    ParametricPlot[{x, y /. Ellips[[2]]}, {x, -50, 50}, 
     PlotStyle -> Green]}, PlotRange -> All];
Glob = ContourPlot[function[x, y], {x, -70, 50}, {y, -70, 70}, 
  PlotPoints -> 50, Contours -> {1.003, 1.05, 1.5, 2, 5.5}, 
  ContourShading -> None, ContourLabels -> True, PlotRange -> All]
Print@Show[
   ContourPlot[function[x, y], {x, -70, 70}, {y, -70, 70}, 
    Contours -> {Q}, ContourShading -> None, ContourLabels -> True, 
    ContourStyle -> Red, PlotRange -> All]
   , Glob, PlotRange -> All];

I think the region of interest is the interior of the red ellipse (if there is no mistake in my reasoning!).

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract