NIntegrate or NExpectation a vector function

Question

Question How do I use NIntegrate/NExpectation to integrate over a vector function for which each element is a function of a common solution to an underlying equation?

Background I know that NIntegrate can integrate vectors, but according to a comment to this question, the vector structure of the output must be particularly obvious to NIntegrate. The comment constrasts NIntegrate[{Sin[x],Cos[x]},{x,0,1}], which will work, with f[x_?NumericQ]:={Sin[x],Cos[x]}; NIntegrate[f[x],{x,0,1}], which does not work. In the example that works, the vector can be written in an obvious way because each element can be calculated completely separately from the others. I don't know how to make the vector structure apparent to NIntegrate when each element of my vector function needs the root of an equation that I want avoid solving separately for each element of the vector.

Details To be more specific, I need to solve for the smallest root greater than zero of a function of the incomplete Beta function, its first derivative, and a parameter. $$x^*(k) = \min\{x: x> 0 \text{ and } B(1,2,2) - B(x,2,2) - kx(1-x) = 0\}$$ for which $B()$ is the incomplete beta function, and $B^{\prime}() = x(1-x)$ is its first derivative.

Then, using this solution, I calculate a vector of objects of interest (for simplicity, I'm focusing on just two elements in my example) $\{\frac{B(x^*(k),2,2)}{(B(1,2,2)}, \frac{B(x^*(k),3,2)}{B(x^*(k),2,2)}\}$

Finally, I want to give the parameter a distribution and integrate the vector of interest over that distribution.

Clear[x,equation, startVals, modelOutcomes];
equation[x_?NumericQ, param_] := Beta[2,2] - Beta[x, 2, 2] - x*(1 - x)*param;

startVals = Range[1/10, 9/10, 1/10]; (*initial values for  FindRoot*)

modelOutcomes[param_?NumericQ] := 
  Module[{xSoln, result1, result2},

   (*find the smallest root greater than 0 *)
   xSoln = Min[
     Part[
      FindRoot[equation[xVal, param], {xVal, #, 0, 1}] & /@ startVals,
      All, 1, 2]
     ]; 

   (*calculate results*)
   result1 =  Beta[xSoln, 2, 2] / Beta[2, 2];
   result2 = Beta[xSoln, 3, 2] / Beta[xSoln, 2, 2];

   {result1, result2}
   ];

mu = 1/10; sigma = 1/1000;

NExpectation[modelOutcomes[z1] , z1 \[Distributed] LogNormalDistribution[mu, sigma]]
(* Integrand modelOutcomes[E^(1/10) z1] Piecewise[{{(500 E^(-500000 Log[z1]^2) Sqrt[2/\[Pi]])/z1,z1>0}},0]
is not numerical at {z1}={250.34090521920675`}*)

So as I have written the code, NIntegrate is unable to see and expect a vector output from the Integrand, telling me that the output is "not numerical."

If I alter modelOutcomes so that it just result1, then NExpectation works. If I alter it to return a {result1}, it gives me the same "not numerical" error. So the problem seems to be about getting NIntegrate to expect a vector integrand.

Answer 1

Maybe something like the following will suffice:

out1[x_?NumericQ] := First[{out1[x], out2[x]} = modelOutcomes[x]]

NExpectation[{out1[z1], out2[z1]}, z1 \[Distributed] LogNormalDistribution[mu,sigma]]

{0.0374898, 0.0550367}

Answer 2

We can write the integrand as a nested pure function, with the outer function taking, as an argument, the variable of integration, and the inner pure function taking the solution of the underlying equation and using it to produce an explicit list of the resulting objects of interest.

Clear[x, equation, startVals, modelOutcomes];
startVals = Range[1/10, 9/10, 1/10];
equation[x_?NumericQ, param_] := 
  Beta[2, 2] - Beta[x, 2, 2] - x*(1 - x)*param;

xSoln[param_?NumericQ] := 
  Min[Part[FindRoot[equation[xVal, param], {xVal, #, 0, 1}] & /@ 
     startVals, All, 1, 2]];

result1 = Beta[#, 2, 2]/Beta[2, 2] &;
result2 = Beta[#, 3, 2]/Beta[#, 2, 2] &;
modelOutcomes = {result1[#], result2[#]} &[xSoln[#]] &;

mu = 1/10; sigma = 1/1000;

NExpectation[modelOutcomes[z1], 
 z1 \[Distributed] LogNormalDistribution[mu, sigma]]
(*{0.0374898, 0.0550367}*)

This approach has a disadvantage, which is that NIntegrate is essentially separately integrating each element of the integrand vector. I can't demonstrate this with NExpectation, because it doesn't have the EvaluationMonitor option, but I can demonstrate it with NIntegrate. So for demonstration purposes I'll use a different function, but using the same nested pure functions approach to the vector integrand

Clear[g, f, samplePointsList];
g[y_] := y^2;
f := {Sin[#], Cos[#]} &[g[#]] &

Borrowing an approach from this question, I collect the evaluation points used by NIntegrate, and plot the point on the horizontal axis with its position in the order of evaluated points on the vertical axis,

samplePointsList = 
  Reap[NIntegrate[f[x], {x, 0, 1}, EvaluationMonitor :> Sow[x]
     ]][[2, 1]];
ListPlot[Transpose[{samplePointsList, 
   Range[Length[samplePointsList]]}]]

It looks like the evaluation happens twice at essentially the same set of points. What fraction of the evaluation points are unique?

N@Length[Union[samplePointsList]]/Length[samplePointsList] 
(*0.505792*)

So this approach facilitates the use of NIntegrate / NExpectation, but doesn't deal with the duplicated evaluation problem.