Speed of convergence for NIntegrate

Question

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, ∞}]

Stefan · Accepted Answer · 2013-07-09T15:02:59.990

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 — OO_SE, Jul 09 '13 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). — Andrew Moylan, Jul 09 '13 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). — Stefan, Jul 09 '13 at 18:49

score 5 · Answer 2 · edited Jul 09 '13 at 19:07

5

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[1][#]} & /@ Range[3]), 
 PlotRange -> All, PlotPoints -> 50]

enter image description here

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

edited Jul 09 '13 at 19:07

Stefan

5,347
25
32

answered Jul 09 '13 at 18:59

Andrew Moylan

4,260
22
25

nice graphical analysis! (+1) from me...just installed the @halirutan java-script thing...you may have a look...works like a breeze. http://meta.mathematica.stackexchange.com/a/1044/5478 – Stefan Jul 09 '13 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! – OO_SE Jul 10 '13 at 09:41
Here is a more detailed exposition of the problem I'm trying to solve. Thanks – OO_SE Jul 10 '13 at 13:40

Speed of convergence for NIntegrate

2 Answers2