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.

Context

I am interested in computing numerically the number of extrema at a given threshold for random fields. These numbers are expectations of MultinormalDistributions.

Problem

This integral should give the differential number of maxima of a scale invariant Gaussian random field ($\gamma=0$ i.e. white noise) at threshold zero.

NExpectation[(-(x*y) + z^2)*Boole[x*y - z^2 <= 0], 
{x, z, y} \[Distributed] 
MultinormalDistribution[{0, 0, 0}, {{3/8, 0, 1/8}, {0, 1/8, 0}, {1/8, 0, 3/8}}]]

Its value should be 1/4/Sqrt[3]= 0.144338 whereas Mathematica returns 0.138245 which is wrong by a factor of 5 %.

Strangely enough if I split the integral in two as

NExpectation[(-(x*y) + z^2)*Boole[x + y >= 0 && x*y - z^2 <= 0], 
 {x, z, y} \[Distributed] 
 MultinormalDistribution[{0, 0, 0},{{3/8, 0, 1/8}, {0, 1/8, 0}, {1/8, 0, 3/8}}]] +
NExpectation[(-(x*y) + z^2)*Boole[x + y <= 0 && x*y - z^2 <= 0], 
  {x, z, y} \[Distributed] 
MultinormalDistribution[{0, 0, 0}, {{3/8, 0, 1/8}, {0, 1/8, 0}, {1/8, 0, 3/8}}]]

Mathematica returns the correct answer to 5 digits.

I understand integrating numerically Boole is likely to cause problems but I am nonetheless surprised by the poor result.

On a related topic

 NExpectation[-(x*y) + z^2 \[Conditioned] (x*y - z^2 < 0), 
 {x, z, y} \[Distributed] 
 MultinormalDistribution[{0, 0, 0}, {{3/8, 0, 1/8}, {0, 1/8, 0}, {1/8, 0, 3/8}}]]

seems to take a while before returning unevaluated (even though it should probably be the correct way to proceed?).

Questions

  1. Why is the first method doing so poorly?
  2. Why is the second method failing altogether?
  3. Any guidance on how to tell NExpectation to be "smart" about the Boole condition?

Thanks

Regards

share|improve this question
    
Did you read reference.wolfram.com/mathematica/tutorial/… already? –  belisarius Jul 31 '12 at 12:23
add comment

1 Answer

up vote 4 down vote accepted

Numerical integration in Mma is a big topic. I suggest reading at least this.

The following gives the correct answer:

NExpectation[(-(x*y) + z^2) Boole[x*y - z^2 <= 0], {x, z, y} \[Distributed] 
  MultinormalDistribution[{0, 0, 0}, {{3/8, 0, 1/8}, {0, 1/8, 0}, {1/8, 0, 3/8}}], 
     Method -> {"NIntegrate", {Exclusions -> True}}]
share|improve this answer
    
Thanks a lot! Will look at your link. Do you happen to know why the Conditioned formulation fails altogether? –  chris Jul 31 '12 at 12:57
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.