5
$\begingroup$

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]]}]]];

tangent line to contour plot of PDF

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}];

enter image description here

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?

$\endgroup$
3
  • 1
    $\begingroup$ 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. $\endgroup$
    – JimB
    Sep 14, 2018 at 17:23
  • $\begingroup$ @JimB. Thank you, will give this a go. I may be able to eliminate the second one by reducing the WorkingPrecision specification. $\endgroup$
    – GerardF123
    Sep 14, 2018 at 19:45
  • 1
    $\begingroup$ I think you're missing the definition of contourval. It should be contourval = PDF[jointdist, coeffNorm]. $\endgroup$
    – JimB
    Sep 14, 2018 at 23:04

2 Answers 2

2
$\begingroup$
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.

$\endgroup$
3
$\begingroup$

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]

enter image description here

   0.307558
$\endgroup$
5
  • $\begingroup$ 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? $\endgroup$
    – GerardF123
    Sep 14, 2018 at 19:42
  • 1
    $\begingroup$ @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]. $\endgroup$
    – Stratus
    Sep 14, 2018 at 19:58
  • $\begingroup$ Actually I think I see what may be going wrong: g[x_,y_] is a function only of x. $\endgroup$
    – GerardF123
    Sep 14, 2018 at 21:13
  • $\begingroup$ And then when I use that formulation, the NIntegrate::slwcon error re-emerges. $\endgroup$
    – GerardF123
    Sep 14, 2018 at 21:17
  • 1
    $\begingroup$ 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. $\endgroup$
    – JimB
    Sep 14, 2018 at 23:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.