I need to solve for the smallest root greater than zero of a function of the incomplete Beta function, its first derivative, and a few parameters. Then, using this solution, I calculate a vector of objects of interest.

Clear[f, x, alpha, beta, discount, gamma, probWDraw, impliedCapF0, startVals, modelOutcomes]; (*initialize some primitives*) discount = 97/100; gamma = 1/3; probWdraw = 15/100; (*define needed function*) f[x_, alpha_, beta_] = x^(alpha - 1)*(1 - x)^(beta - 1); impliedCapF0[x_, alpha_, beta_, cj_, gammaV_] := Beta[x, alpha, beta] + f[x, alpha, beta]*cj/gammaV; (*pick starting values for `FindRoot`*) startVals = Range[1/10, 9/10, 1/10]; modelOutcomes[alpha_, beta_, v_, cost1_?NumericQ, cost2_?NumericQ] := Module[ {capF0 = Beta[alpha, beta], tempSolnList, xSoln, capFxSoln, cj = cost1 + cost2*gamma, gammaV = gamma*v, result1, result2, result3, result4, result5}, (*find the smallest root greater than 0 *) tempSolnList = Part[ FindRoot[ capF0 - impliedCapF0[xVal, alpha, beta, cj, gammaV], {xVal, #, 0, 1}] & /@ startVals, All, 1, 2]; xSoln = Min[tempSolnList]; (*calculate results*) capFxSoln = Beta[xSoln, alpha, beta]/capF0; result1 = 1 - capFxSoln; result2 = capFxSoln*probWdraw; result3 = (1 - probWdraw)*capFxSoln; result4 = Beta[xSoln, alpha + 1, beta] / Beta[xSoln, alpha, beta]; result5 = discount*(v*xSoln + cost2); (*return result list*) {result1, result2, result3, result4, result5} ]; 

I want to give the parameters a joint density and integrate the vector of interest over that density.

mu = {8, 8}; sigma = {{1/800, 1/1000}, {1/1000, 1/800}}; NExpectation[ modelOutcomes[5, 3, 80000, z1, z2] , {z1, z2} \[Distributed] LogMultinormalDistribution[mu, sigma]] (* Integrand modelOutcomes[5,3,80000,E^8 z1,E^8 z2] Piecewise[{{(2000 E^(1/2 (-Log[z1] ((20000 Log[z1])/9-(16000 Log[z2])/9)-Log[z2] (-((16000 Log[z1])/9)+(20000 Log[z2])/9))))/(3 \[Pi] z1 z2),z1>0&&z2>0}},0]] is not numerical at {z1,z2}={3.`,3.`}*) 

Why am I being told the Integrand is not numerical? I'm using the ?NumericQ pattern test on the arguments that are supplied by NExpectation / NIntegrate. And if I plug in the values shown in the stack trace for the error, namely 3*Exp[8] for z1 and z2, the function has no problems producing a numerical result

modelOutcomes[5, 3, 80000, 3*Exp[8], 3*Exp[8]] (*{0.758906, 0.0361641, 0.20493, 0.405076, 48150.2}*) 

And NIntegrate doesn't seem to have a problem with the vector per se, because it integrates a constant vector without problem

NExpectation[ modelOutcomes[5, 3, 80000, 10000, 10000] , {z1, z2} \[Distributed] LogMultinormalDistribution[mu, sigma]] (*{0.788271, 0.0316643, 0.179431, 0.391449, 47730.3}*) 

What do I need to change so that NExpectation handles this without error?