# 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:
Answer
7 Replies
Sort By:
Posted 1 year ago
 try testing your LKLH function.
Answer
Posted 1 year ago
 I'm not sure what you mean with this suggestion. Could you please clarify? Thanks :-)
Answer
Posted 1 year ago
 try evaluating it with numerical arguments
Answer
Posted 1 year ago
 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}]
Answer
Posted 1 year ago
 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!).
Answer
Posted 1 year ago
 Attachments:
Answer
Posted 1 year ago
Answer
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments