NIntegrate discrepancies

Question

I want to compute a numerical integral for a gaussian trivariate distribution. Say

NIntegrate[
  X* PDF[MultinormalDistribution[{0,0,1}, {{1, 1/10, 1/5}, {1/10, 1,1/10}, {1/5, 1/10, 1}}], {X, Y, Z}],
  {X, -Infinity, Infinity}, {Y, 0, Infinity}, {Z, -Infinity, -1}
]

(*Out: -0.00365287 *)

This is quite expensive to calculate, so I'd like to make use of parameterisations available for the Gaussian distribution. Specifically, the fact that:

$$ E[X\mid Y\geq y, Z\leq z]=E[E[X \mid Y,Z]\mid Y\geq y, Z\leq z]\\ =\beta_{XY;Z} E[Y\mid Y\geq y, Z\leq z]+\beta_{XZ;Y} E[Z \mid Y\geq y, Z\leq z] $$

where $\beta_{XY;Z}$ notes the regression coefficient of $X$ on $Y$ conditional on $Z$ (and $\beta_{XZ;Y}$ is analogously defined).

Applying this parameterisation to calculate the same integral from above gives

ρXY = 1/10;
ρXZ = 2/10;
ρYZ = 1/10;

βXYz = (ρXY - ρXZ*ρYZ)/(Sqrt[(1 - (ρXZ^2))] Sqrt[(1 - (ρYZ^2))]);
βXZy = (ρXZ - ρXY*ρYZ)/(Sqrt[(1 - ρXY^2)] Sqrt[(1 - ρYZ^2)]);

βXYz*NIntegrate[Y*PDF[BinormalDistribution[{0, 1}, {1, 1}, ρYZ], {Y, Z}], {Y, 0, Infinity}, {Z, -Infinity, -1}] 
+ 
βXZy*NIntegrate[Z*PDF[BinormalDistribution[{0, 1}, {1, 1}, ρYZ], {Y, Z}], {Y, 0, Infinity}, {Z, -Infinity, -1}]


(* Out: -0.00187282 *)

Is it possible that this discrepancy is coming from the precision of NIntegrate? Or am I just doing something wrong here?

NB. The two computational procedures coincide perfectly if ρXY= ρXZ= ρYZ

---

UPDATE: I am interested in speeding up the first integral, because what I would ultimately like to do is fix ρXY and ρXZ and calculate something like

   ContourPlot[ NIntegrate[
 X*PDF[MultinormalDistribution[{0,0,µ}, {{1, ρXY, ρXZ}, {ρXY, 1, ρYZ}, {ρXZ, ρYZ, 1}}], {X, Y, Z}], {X, -Infinity, Infinity}, {Y, a, Infinity}, {Z, -Infinity,b}],{µ,-2,2},{ρYZ,-.9,-9} ]

where the limits $a$ and $b$ are themselves numerical calculations (FindRoots). In other words, I would ideally like a symbolic expression of the first integral in terms of µ, ρYZ, $a$ and $b$

Answer 1

A one-time lengthy (one-minute) integral of the general integrand can reduce the dimension of NIntegrate by 2 to a fairly long expression (containing Erfi and Sqrt). Doing the Z integral symbolically seemed impossible. The dimension reduction speeds up NIntegrate by a factor of around 70 on my machine.

intXY = Integrate[
  X*PDF[MultinormalDistribution[{0, 0, μ},
     {{1, ρXY, ρXZ}, {ρXY, 1, ρYZ}, {ρXZ, ρYZ, 1}}], {X, Y, Z}],
  {Y, a, Infinity}, {X, -Infinity, Infinity}, 
  Assumptions -> {a, b, μ, ρXY, ρXZ, ρYZ} ∈ Reals, 
  GenerateConditions -> False]
(*  (E^(-(1/2) (Z - μ)^2) (<...>))/(4 Sqrt[2] π <...>)  *)

Check of OP's example: Something causes NIntegrate to use complex numbers, but the imaginary part comes out to be zero, which it should be. You can get twice the speed by applying Re to the integrand.

Block[{ρXY = 1/10, ρXZ = 2/10, ρYZ = 1/10,
   a = 0, b = -1, μ = 1},
 NIntegrate[Re@intXY, {Z, -Infinity, b}] // AbsoluteTiming]
(*  {0.014633, -0.00365287}  *)

No idea about appropriate parameters for the ContourPlot, but here goes:

Block[{ρXY, ρXZ, ρYZ, a, b, μ, NIntegrate},
  obj[
   ρXY_?NumericQ, ρXZ_?NumericQ, ρYZ_?NumericQ, a_?NumericQ, b_?NumericQ, μ_?NumericQ] = 
     NIntegrate[Re@intXY, {Z, -Infinity, b}];
  ];
ContourPlot[obj[1/10, 2/10, ρYZ, 0, -1, μ],
  {μ, -2, 2}, {ρYZ, -.9, .9}, 
  MaxRecursion -> 1] // AbsoluteTiming