Message Boards Message Boards

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

Posted 6 years ago
Attachments:
POSTED BY: daniele barducci
7 Replies
Anonymous User
Anonymous User
Posted 6 years ago

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

POSTED BY: Anonymous User

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

Attachments:
POSTED BY: daniele barducci

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

POSTED BY: Claude Mante
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}]
POSTED BY: Frank Kampas

try evaluating it with numerical arguments

POSTED BY: Frank Kampas

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

POSTED BY: daniele barducci

try testing your LKLH function.

POSTED BY: Frank Kampas
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