# Integrate over FindRoot solutions

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?

• Because Integrate[f2[X0, u], {u,-Infinity, Y0}] yields an analytic answer, replace it by that answer in FindRoot to save a lot of time. – bbgodfrey Jan 17 '16 at 15:33

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 InverseErfc[(E^(-(X0^2/2)) Erfc[(X0 - 2 Y0)/Sqrt])/Sqrt[2 π]])  *)


The warning perhaps calls for some numerical checks:

With[{range = Join[-Reverse@#, {0}, #] &[2^Range]},
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  *)

• Thanks for the reply. In fact, I have tried Solve function and it worked well as you indicated. However, I have several more complicated equations to be solved which Solve function doesn't handle. This is why I tried to experiment with FindRoot which is quite slow together with NIntegrate. – wlq Jan 17 '16 at 17:12