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