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.
Evaluated -> False
toFindRoot
might help. $\endgroup$