1
$\begingroup$

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.

$\endgroup$
  • 1
    $\begingroup$ Are you sure your example can't be simplified further? That's a lot of code to wade through! $\endgroup$ – MarcoB Jul 27 '17 at 19:50
  • $\begingroup$ @MarcoB Just a moment. I'll work on cutting it down. $\endgroup$ – FalafelPita Jul 27 '17 at 19:57
  • $\begingroup$ Adding the option Evaluated -> False to FindRoot might help. $\endgroup$ – bbgodfrey Jul 27 '17 at 22:53
  • $\begingroup$ @bbgodfrey Thanks. I tried, but it didn't fix things. $\endgroup$ – FalafelPita Jul 28 '17 at 0:13
  • $\begingroup$ This and this might be of interest. $\endgroup$ – J. M. is in limbo Jul 28 '17 at 13:38
2
$\begingroup$

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}

$\endgroup$
0
$\begingroup$

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]]}]]

enter image description here

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.

$\endgroup$
  • 1
    $\begingroup$ Use SetOptions[NIntegrate, EvaluationMonitor :> Sow[z1]] to set EvaluationMonitor for NIntegrate before calling NExpectation (which uses NIntegrate to do the actual integration) in order to obtain the evaluations used to by NExpectation. $\endgroup$ – bbgodfrey Jul 28 '17 at 19:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.