Integrate over FindRoot solutions

Question

I have a function of bivariate normal PDF and its marginals defined as

f[x_] := PDF[NormalDistribution[0, 1], {x}]
f2[x_, y_] := 
  PDF[MultinormalDistribution[{0, 0}, {{1, ρ}, {ρ, 1}}], {x, y}] /. {ρ -> 1/2}

The function fw is quite complex so I solve it numerically using FindRoot for given X0 and Y0.

fw[X0_?NumericQ, Y0_?NumericQ] := 
   Y /. 
     FindRoot[
       Integrate[f2[X0, u], {u, -∞, Y0}]f[X0] == Integrate[f2[X0, u], {u, -∞, Y}], 
       {Y, Y0}]

I would like to calculate the expected value of the function fw over the Gaussian curve but for a reason this doesn't work.

NIntegrate[fw[XD, YD] f2[XD, YD], {XD, -∞, ∞}, {YD, -∞, ∞}]

I suppose NIntegrate is either too slow or it cannot supply numerical values to fw. Do you have any ideas how to fix this?

Answer 1

As @bbgodfrey commented, if the integrals in the equation in the OP's FindRoot command can be evaluated before passing the equation to FindRoot, one can save a lot of time. It seems there is still more to be done. I found FindRoot struggles to find an accurate root in some areas of the domain of the equation. It turns out one can use Solve to solve the equation and avoid FindRoot altogether.

fwEQ = Integrate[f2[X0, u], {u, -Infinity, Y0}] f[X0] ==
  Integrate[f2[X0, u], {u, -Infinity, Y}];
fwII[X0_, Y0_] = Y /. First@Solve[fwEQ, Y]

Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information. >>

(*  1/2 (X0 -Sqrt[6] InverseErfc[(E^(-(X0^2/2)) Erfc[(X0 - 2 Y0)/Sqrt[6]])/Sqrt[2 π]])  *)

The warning perhaps calls for some numerical checks:

With[{range = Join[-Reverse@#, {0}, #] &[2^Range[10]]},
   Table[
    fwEQ /. Equal -> Subtract /. {X0 -> x, Y0 -> y} /. Y -> fwII[x, y],
    {x, range}, {y, range}]
   ] // Abs // Max
(*  0  *)

The solution seems valid, so time to integrate.

Mathematica does not return an exact solution using Integrate, so we'll resort to a numerical approximation:

NIntegrate[fwII[XD, YD] f2[XD, YD], {XD, -Infinity, Infinity}, {YD, -Infinity, Infinity}]
(*  -1.08854  *)