# Integrating a bivariate distribution over a region bounded by a straight line

Summary of problem: I'm using Mathematica version 8 to try to integrate the bivariate distribution over a region bounded by a straight line. The two random variables are uncorrelated. When I use Nintegrate, I get the following error message:

"NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small."

Do I simply have to live with this problem or is there something I can do to get rid of the error message and gain assurance that my result is correct?

Model: I give an example using the Bivariate Normal, though I find that this problem also exists with other distributions linked using the product copula.

Assume that the bivariate distribution is given by:

mu1 = 0; sig1 = 0.1; mu2 = 0; sig2 = 0.05; corr = 0; jointdist =

MultinormalDistribution[{mu1, mu2}, {{sig1^2, 0}, {0, sig2^2}}];


Definition of the region: We have a point on the Cartesian plane which we will use to determine the straight line that defines the region. This is calculated as

coeffNorm = {0.0485058025031513, -0.006616445600455179};


The straight line (defining the region to be integrated) is the projection of the tangent to the joint PDF going through that point. This is calculated and visualized as follows (on the Cartesian plane):

PDF[jointdist, {x, y}];

myFxy[x_, y_] :=
31.830988618379067 E^(
1/2 (-x (0. + 99.99999999999999 x) -
y (0. + 399.99999999999994 y)));

xVec = D[myFxy[x, y], x] /. {x -> coeffNorm[[1]], y -> coeffNorm[[2]]};

yVec = D[myFxy[x, y], y] /. {x -> coeffNorm[[1]], y -> coeffNorm[[2]]};

Show[RegionPlot[
PDF[jointdist, {x, y}] >= contourval, {x, -0.15, 0.15}, {y, -0.05,
0.05}], Plot[-xVec/yVec*(x - coeffNorm[[1]]) +
coeffNorm[[2]], {x, -7, 3}],
Graphics[Point[{coeffNorm[[1]], coeffNorm[[2]]}]]];


The region over which the PDF to be integrated is visualized as follows (though it extends to infinity in both directions):

RegionPlot[
y <= -xVec/yVec*(x - coeffNorm[[1]]) + coeffNorm[[2]], {x, -0.15,
0.15}, {y, -0.05, 0.05}];


Integral: So the probability that I wish to calculate is given by:

NIntegrate[
PDF[jointdist, {x, y}]*
Boole[y <= -xVec/yVec*(x - coeffNorm[[1]]) +
coeffNorm[[2]]], {x, -Infinity, Infinity}, {y, -Infinity,
Infinity}, WorkingPrecision -> 20];


Error message: I get a numerical output (a value of 0.30755790488620519025), but an error message also. I've extended the WorkingPrecision to 20 as a result of getting the error, but this doesn't seem to help. I get the following errors:

"NIntegrate::precw: The precision of the argument function (31.831 E^(1/2 (-x (0. +100. x)-y (0. +400. y))) Boole[y<=-0.00661645+1.83277 (-0.0485058+x)] ) is less than WorkingPrecision (20.). >>
NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small. >>"

Can you suggest how I can improve my code to get rid of these errors? I would not have expected the second one, given that the Normal distribution is quite well behaved. Regarding the first one, does it mean that I need to improve the working precision of all of the calculations that I use to derive the numbers going into the specification of the distribution, or is there something else I should do? How accurate is the integral that is being calculated?

• If you Rationalize all of the constants, you'll get rid of the first warning but not the second. In other words, use sig1 = 1/10; and Rationalize[coeffNorm,0], etc. – JimB Sep 14 '18 at 17:23
• @JimB. Thank you, will give this a go. I may be able to eliminate the second one by reducing the WorkingPrecision specification. – GerardF123 Sep 14 '18 at 19:45
• I think you're missing the definition of contourval. It should be contourval = PDF[jointdist, coeffNorm]. – JimB Sep 14 '18 at 23:04

expr = y <= -xVec/yVec*(x - coeffNorm[[1]]) + coeffNorm[[2]];


With the default value of 0 for the option MinRecursion

NIntegrate may miss sharp peaks of integrands

Setting a larger value for this option forces a finer subdivision of the integration region:

NIntegrate[PDF[jointdist, {x, y}] Boole[expr], {x, -∞, ∞}, {y, -∞, ∞}, MinRecursion -> 5]


0.3075579043682307

We get the same result using NProbability and NExpectation with the option Method -> {"NIntegrate", MinRecursion -> 5}:

NProbability[expr, Distributed[{x, y}  jointdist],
Method -> {"NIntegrate", MinRecursion -> 5}]


0.3075579043682307

NExpectation[Boole@expr, Distributed[{x, y}  jointdist],
Method -> {"NIntegrate", MinRecursion -> 5}]


0.3075579043682307

If you Rationalize the numeric inputs as suggested by JimB

{mu1, mu2, sig1, sig2, coeffNorm} = Rationalize[{mu1, mu2, sig1, sig2, coeffNorm}, 0];


and add the option WorkingPrecision ->20, NIntegrate gives

0.30755790488636425134477259747633448173555690466295398523820.

and NProbability and NExpectation both give

0.30755790488636425134477259748716002931498332884209285513820.

EDITED

As suggested by JimB I added the -y to gx[x_,Y_].

 mu1 = 0; sig1 = 0.1; mu2 = 0; sig2 = 0.05; corr = 0;jointdist =
MultinormalDistribution[{mu1, mu2}, {{sig1^2, 0}, {0, sig2^2}}];
cplot = ContourPlot[{PDF[jointdist, {x, y}]}, {x, -0.15,
0.15}, {y, -0.05, 0.05}, Contours -> 10];
myFxy[x_, y_] :=
31.830988618379067 E^(1/
2 (-x (0. + 99.99999999999999 x) -
y (0. + 399.99999999999994 y)));
coeffNorm = {0.0485058025031513, -0.006616445600455179};

xVec = D[myFxy[x, y], x] /. {x -> coeffNorm[[1]], y -> coeffNorm[[2]]};
yVec = D[myFxy[x, y], y] /. {x -> coeffNorm[[1]], y -> coeffNorm[[2]]};
gx[x_, y_] = -xVec/yVec*(x - coeffNorm[[1]]) + coeffNorm[[2]] - y
rgplot = RegionPlot[gx[x, y] >= 0, {x, -0.15, 0.15}, {y, -0.05, 0.05}];
Show[cplot, rgplot]

NProbability[
gx[x, y] >= 0 && -Infinity < x < Infinity && -Infinity < y <
Infinity, {x, y} \[Distributed] jointdist]


   0.307558

• Thank you for this @Diogo. ImplicitRegion is not available in Mathematica v8. Do you know of a solution that uses v8 vocabulary? Another comment is that the straight line in your plot looks as though it is vertical, but it shouldn't be. Do you understand why this may be? – GerardF123 Sep 14 '18 at 19:42
• @GerardF123 you are welcome. I don't know why in my plot it's a straigtht line, I checked the equations and the are the same, maybe it's a bug. In v8 you can use NProbability[ gx[x, y] >= 0 && -Infinity < x < Infinity && -Infinity < y < Infinity, {x, y} \[Distributed] jointdist]. – Diogo Sep 14 '18 at 19:58
• Actually I think I see what may be going wrong: g[x_,y_] is a function only of x. – GerardF123 Sep 14 '18 at 21:13
• And then when I use that formulation, the NIntegrate::slwcon error re-emerges. – GerardF123 Sep 14 '18 at 21:17
• If you add ListPlot[{coeffNorm}] to Show, you'll see that the line you've drawn is not the desired line as it doesn't go through coeffNorm. Part (or all) of the issue is that y isn't part of the right-hand side of gx[x_,y_]. You want gx[x_, y_] := -xVec/yVec*(x - coeffNorm[[1]]) + coeffNorm[[2]] - y. – JimB Sep 14 '18 at 23:00