4
$\begingroup$

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, ∞}]
Improve this question
$\endgroup$

2 Answers 2

Reset to default
6
$\begingroup$

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}
Improve this answer
$\endgroup$
3
  • $\begingroup$ 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 $\endgroup$
    – OO_SE
    Commented Jul 9, 2013 at 11:03
  • $\begingroup$ 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). $\endgroup$
    – Andrew Moylan
    Commented Jul 9, 2013 at 18:45
  • $\begingroup$ @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). $\endgroup$
    – Stefan
    Commented Jul 9, 2013 at 18:49
5
$\begingroup$

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

Improve this answer
$\endgroup$
3
  • $\begingroup$ 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 $\endgroup$
    – Stefan
    Commented Jul 9, 2013 at 19:11
  • $\begingroup$ 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! $\endgroup$
    – OO_SE
    Commented Jul 10, 2013 at 9:41
  • $\begingroup$ Here is a more detailed exposition of the problem I'm trying to solve. Thanks $\endgroup$
    – OO_SE
    Commented Jul 10, 2013 at 13:40

Your Answer

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

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