Message Boards Message Boards


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

Posted 10 months ago
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 ( I tried to do something like

   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

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]]/
  Poisson[function[x, y], function[0, 0]], {x, -500, 500}, {y, -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}};
 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:

    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]
   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!).

Hi, Thanks all. I'm not sure the last soluction works. Beside the fact that the LKLH that you define is not normalized to one, the contour you are showing are for constant function[x,y], and not for constant integral of LKLH[x,y] within an ellipses. In fact if you plot the 1, 2, 3, 4, 5... sigma contour they are not "closer to each other. LKLH is not even used there, but maybe I misunderstood the solution. In any case I managed to do something. I attach it here since it might be useful to someone else. The "trick" was to define a table for the values of the integral of LKLH within for different elippses "radius" instad than using FindRoot. Clearly one needs to know which value of "radius" can give results close to the desired ones, but in a real problem this is always the case. Then I interpolate the table and find the function that gives me the "radius" for different confidence level values. Then I plot them


i haven't followed the notebook "closely". but if mathematica can't Integrate or NIntegrade f[x,y] to show the (volume), then how will you (use mathematica) to determine if any ellipse is within %95 of the correct answer which was not yet found?

In[1113]:= Clear[x, y, t]

In[1114]:= x[t_] := t;

In[1115]:= y[t_] := t;

In[1116]:= Integrate[
 1 + 1.5 10^-3 x[t] + 2.5 10^-3 x[t]^2 + 1.7 10^-3 y[t] + 
  1.8 10^-3 y[t]^2 + 1.8 10^-3 x[t] y[t], t]

Out[1116]= t + 0.0016 t^2 + 0.00203333 t^3

or on the other hand,

In[1117]:= Integrate[
 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, {x, 0, 1}, {y, 0, 1}]

Out[1117]= 1.00348

what i'm not seeing here is why you don't use a double integral (for volume) which is solvable, or set up an implicit equation and use implicit rules

can either or both of x and y be independant? can they both be written that way? are functions allows as parameters to this function or only numbers? it's possible for some functions of f[x,y] to be written in the form of a family of ellipses (with multple as argument). are you intending that f[x,y] be any function, or is f[x,y] the only function you need evaluated?

i think a little more needs to be know about your specifications before one can conclude regions

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

Group Abstract Group Abstract