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

Question

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?

Answer 1

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.307557904886364251344772597476334481735556904662953985238`20.

and NProbability and NExpectation both give

0.307557904886364251344772597487160029314983328842092855138`20.

Answer 2

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