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
- Why is the first method doing so poorly?
- Why is the second method failing altogether?
- Any guidance on how to tell NExpectation to be "smart" about the Boole condition?
Thanks
Regards
