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

Posted 4 months ago
615 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
Sort By:
Posted 4 months ago
Posted 4 months ago
 I'm not sure what you mean with this suggestion. Could you please clarify? Thanks :-)
Posted 4 months ago
 try evaluating it with numerical arguments
Posted 4 months 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}] 
Posted 4 months 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!).
 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 rulescan 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