Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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, ∞}]
share|improve this question
add comment

2 Answers

up vote 5 down vote accepted

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}
share|improve this answer
    
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 –  EOO Jul 9 '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 9 '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 9 '13 at 18:49
add comment

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

share|improve this answer
    
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 –  Stefan Jul 9 '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! –  EOO Jul 10 '13 at 9:41
    
Here is a more detailed exposition of the problem I'm trying to solve. Thanks –  EOO Jul 10 '13 at 13:40
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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