# Speed of convergence for NIntegrate

I'm trying to optimise numerically a function that entails computing the expected value of a truncated trivariate normal distribution and this is taking extremely long -I also get warned about NIntegrate::slwcon

I've seen posts (such as this) about rewriting the problem as a differential equation, however I couldn't rewrite the problem successfully.

Here's a simplified version of my problem (focusing on the expectation, which is what slows down the optimisation and produces NIntegrate::slwcon ):

JointD[X_,Y_,Z_]:= PDF[MultinormalDistribution[
{0,0,0}, {{1, 0.5, 0.5}, {0.5, 1, 0.5}, {0.5, 0.5, 1}}], {X,Y,Z}]

NIntegrate[X*JointD[X,Y,Z],{X, -∞, ∞}, {Y, 1, ∞}, {Z, 2, ∞}]


WorkingPrecision of NIntegrate

If we read the warning, it states:

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

Thanks for the advice. So let's increase the WorkingPrecision:

JointD[X_, Y_, Z_] :=
PDF[MultinormalDistribution[{0, 0, 0},
{{1, 0.5, 0.5}, {0.5, 1, 0.5}, {0.5, 0.5, 1}}], {X, Y, Z}]

NIntegrate[
X*JointD[X, Y, Z], {X, -∞, ∞}, {Y, 1, ∞}, {Z, 2,∞},
WorkingPrecision -> 10] // Timing

==> {0.328000, 0.01853405413}


Strategies for singularities:

There are several strategies in NIntegrate when you're dealing with singularities:

1) DuffyCoordinates: simplify or eliminate certain types of singularities

2) AdaptiveMonteCarlo:in cases when only a rough integral estimate is needed

etc.

NIntegrate[
X*JointD[X, Y, Z], {X, -∞, ∞}, {Y, 1, ∞}, {Z,
2, ∞},Method -> "DuffyCoordinates"
] // Timing

==> {1.476000, 0.0185346}

NIntegrate[
X*JointD[X, Y, Z], {X, -∞, ∞}, {Y, 1, ∞}, {Z, 2, ∞},
Method -> "AdaptiveMonteCarlo", PrecisionGoal -> 3
] // Timing

==> {0.236000,0.0185348}

• Thanks. Increasing working precision alone does get rid of the message but does not seem to improve the time significantly (at least once you embed it in a more complicated numerical optimisation routine). I'll try the two strategies you mention, and will def read the tutorial on integration strategies Jul 9, 2013 at 11:03
• WorkingPrecision -> 10 amounts to a decrease of working precision from ~16 (the value of \$MachinePrecision). It also decreases the value of PrecisionGoal -> Automatic, to ~4 I think (from a default of ~6). Jul 9, 2013 at 18:45
• @AndrewMoylan the number 10 was not set in stone. it was heuristics to satisfy the warning issue...to be mor e specific, especially in terms of what precision is really needed, there was to less information (for me). Jul 9, 2013 at 18:49

NIntegrate::slwcon doesn't always indicate an error, and I think you can safely Quiet it for this integral.

Roughly speaking, slwcon is emitted when NIntegrate repeatedly subdivides a particular integration subregion and finds the relative error gets worse (or stays the same) instead of getting better.

In your example, I suspect this is happening as NIntegrate subdivides along the X-axis, since there is interesting structure away from X==0:

Plot3D[Evaluate[Table[X*JointD[X, Y, Z], {Z, 2, 3, 0.5}]], {Y, 1,
3}, {X, -2, 2}, Mesh -> None,
PlotStyle -> ({Opacity[0.5], ColorData[#]} & /@ Range),
PlotRange -> All, PlotPoints -> 50] Thus my suggestion would be to transform the integrand such that the peak of the function is near zero:

NIntegrate[(X - 3/2)*
JointD[X - 3/2, Y, Z], {X, -∞, ∞}, {Y,
1, ∞}, {Z, 2, ∞}]

=> 0.0185346


This helps a bit with speed too. Unfortunately multidimensional integrals are pretty much always expensive. But here is one suggestion to explore if performance is crucial: Perhaps you can transform to spherical coordinates (and maybe even do the angular integrals in closed form leaving only a radial integral).

• nice graphical analysis! (+1) from me...just installed the @halirutan java-script thing...you may have a look...works like a breeze. meta.mathematica.stackexchange.com/a/1044/5478 Jul 9, 2013 at 19:11
• Interesting approach. However it turns out that performance is indeed crucial because the problem I actually want to solve is a non-linear constrained optimisation problem over the limits of the integral. I tried transforming the integrand as you suggest (i.e. shifting the peak as a function of the truncation points) but this does not improve performance sufficiently. Your last suggestion is hence relevant. I don't know how to transform the problem to shperical coordinates, so I'll try -and probably open a new post asking about it. Thanks! Jul 10, 2013 at 9:41
• Here is a more detailed exposition of the problem I'm trying to solve. Thanks Jul 10, 2013 at 13:40