1
$\begingroup$

Motivation / Context

I'm working on an economic model which, for present purposes, we can think of as a mapping from parameters to results. Some of the parameters are fixed and given, while others are assumed to be random variables with a known distribution. I would like to calculate the expected value of the results of the model.

Structure of Problem

I've written a simpler outline of my project, which preserves the structure of my approach. var1 is constant that I choose, while var2 has a Log Normal distribution. Given var1 and var2, I have an equation in one variable that I solve using NSolve. I want to calculate the expectation of soln (the solution to the equation) over the distribution of var2 using NExpectation. I can get this to work when the equation solved by NSolve can be solved symbolically, but not when it cannot. I think this happens because NSolve is evaluated before being mapped over whatever table of values is used for numerical integration.

Working Example of Problem

(*Example of the problem*)
Clear[dBeta, solnEquation, solnFun, modelMoment];

dBeta[k_Integer, x_, var1_] = 
  Derivative[k, 0][Beta[#1, #2, 2] &][x, var1];
solnEquation[var1_, var2_, x_] := 
  Sum[Binomial[5, k]*dBeta[k, x, var1]*var2^k, {k, 0, 5}];

solnFun[var1_, var2_] := 
  x /. NSolve[
      solnEquation[var1, var2, x] == Beta[var1, 2] && 0 < x < 1, x, 
      Reals][[1]][[1]];

modelMoment[var1_, mu_, sigma_] := 
  NExpectation[solnFun[var1, var2], 
   Distributed[var2, LogNormalDistribution[mu, sigma]]];

solnFun[2, 1/10]
testMoment[2, 1/10, 1/20]

(*Variation with no problem *)
Clear[solnEquation, solnFun]
solnEquation[var1_, var2_, x_] := var1 - var2 + x;

solnFun[var1_, var2_] := 
  x /. NSolve[solnEquation[var1, var2, x] == 0, x, Reals][[1]][[1]];

solnFun[2, 1/10]
testMoment[2, 1/10, 1/20]

output below

0.328152 
This system cannot be solved with the methods available to NSolve. 
{-20 var2^3+10 var2^2 (1-2 x)+5 var2 (1-x) x+Beta[x,2,2]==1/6} is
neither a list of replacement rules nor a valid dispatch table, and 
so cannot be used for replacing. 
... 
"Further output of \!\(\* StyleBox[RowBox[{\"ReplaceAll\", \"::\", \"reps\"}], \"MessageName\"]\) will 
be suppressed during this calculation." 
NExpectation[  x /. -20 var2^3 + 10 var2^2 (1 - 2 x) + 5 var2 (1 - x) x + 
    Beta[x, 2, 2] == 1/6,   var2 \[Distributed] LogNormalDistribution[1/10, 1/20]]

-1.9
-0.893447

In both cases, solnFun can be evaluated successfully. In the first case, the NExpectation attempt fails, in the second case, it works. I think the difference between the two examples is that the NSolve inside of SolnFun is evaluated before any specific values for var2 are given to it by NExpectation. If a closed form solution exists for symbolic equation, then no error messages appear. But in my project NSolve cannot provide a solution until after it has a specific value for var2.

Question

Is there some way I can get NSolve to be evaluated after the equation has been mapped over whatever table of values NExpectation uses to perform the numeric integration?

I'm new to MMA, so in addition to help with this main problem, I would appreciate any other tips or references on better programming style or approaches.

$\endgroup$
  • 3
    $\begingroup$ Could you provide some code that reproduces the problem (i.e. a simple example without a "symbolic closed form expression")? It sounds like a NumericQ issue. $\endgroup$ – Pragabhava Mar 4 '17 at 2:44
  • $\begingroup$ @Pragabhava I changed solnEquation so that it now reproduces the problem. I was hoping to find a simpler example completely different from my project, but failed to do so after about an hour of trying various things. So instead I just tried to simplify the function I am working with in my project. $\endgroup$ – FalafelPita Mar 4 '17 at 4:06
2
$\begingroup$

Pragabhava's guess is correct: it is indeed the well-known _?NumericQ issue.

Basically you need to define solnFun as solnFun[var1_?NumberQ, var2_?NumberQ] := .... Incorporating also a workaround for NSolve bug, your code can be rewritten as follows:

Clear[dBeta, solnEquation, solnFun];

dBeta[k_Integer, x_, var1_] = Derivative[k, 0][Beta[#1, #2, 2] &][x, var1];
solnEquation[var1_, var2_, x_] = Sum[Binomial[5, k]*dBeta[k, x, var1]*var2^k, {k, 0, 5}];
equation[var1_, var2_, x_] := solnEquation[var1, var2, x] == Beta[2, 2];

solnFun[var1_?NumberQ, var2_?NumberQ] := 
 x /. NSolve[FunctionExpand[equation[2, 1/10, x]] && 0 < x < 1, x, Reals][[1]][[1]];

modelMoment[var1_, mu_, sigma_] := 
  NExpectation[solnFun[var1, var2], Distributed[var2, LogNormalDistribution[mu, sigma]]];

Now

modelMoment[2, 1/10, 1/20]
0.328152
$\endgroup$

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.