Why do Solve, NSolve, and FindRoot all fail on this simple, solvable actuarial equation?

Question

If we express conditional mortality as a vector of annual probabilities of death, like so

qx1={0.04772, 0.05854, 0.07519, 0.09659, 0.11762, 0.13904, 0.16124,
0.18363, 0.2041, 0.22319, 0.24276, 0.262, 0.27897, 0.29458, 0.31044,
0.32691, 0.34597, 0.36573, 0.38348, 0.39799, 0.40855, 0.41447,
0.41774, 0.42266, 0.43064, 0.43913, 0.44417, 0.44802, 0.45, 0.45};

we can compute the associated survival vector with

lx1=Drop[FoldList[#1 (1 + #2) &, 1, -qx1], -1]

and we can find the median life expectancy of this population with

Length[Select[lx1, # >= .5 &]]

We can combine these two equations into a single function that translates the conditional mortality curve directly into a median life expectancy. It would look like

medianLEfromQx[qx_]:=Length[Select[Drop[FoldList[#1 (1 + #2) &, 1, -qx], -1], # >= .5 &]]

And using the sample data above, we see that the median life expectancy for the population is 7 years. If we want to see the implications of reducing annual mortality by 10%, we can ask

medianLEfromQx[.9*qx1]

And we get an answer of 8 years. Fine.

Here lies the problem -- I can't get Mathematica to solve for that multiplier given a desired LE.

Solve[medianLEfromQx[x*qx1]==8,x]

doesn't work (returning an empty set)

NSolve[medianLEfromQx[x*qx1]==8,x]

doesn't work, also returning an empty set, and

FindRoot[medianLEfromQx[x*qx1]==8,{x,1}]

doesn't work either, with the error message that "The Function value {False} is not a list of numbers with dimensions {1} at {x} = {1.}"

What am I doing wrong here?

Answer 1

First of all, check here and here. This is a common pitfall, and questions related to this are posted literally weekly, so I am going to let you review those articles.

In short, if you try evaluating medianLEfromQx[x qx1] with x having no value, you'll see that it returns a number. This expression evaluates inside FindRoot even before FindRoot gets a chance to substitute a value for x. So you would have to make medianLEfromQx not evaluate except for truly numerical vector arguments.

You can do this by changing its definition to look like:

Clear[medianLEfromQx]
medianLEfromQx[qx_ /; (VectorQ[qx, NumericQ])] := ...

Now medianLEfromQx[x qx1] won't evaluate unless x has a numerical value.

---

Next, Solve and NSolve won't work on numerical blackboxes, only FindRoot will. Solve only works with symbolic equations with exact coefficients. NSolve is designed for solving polynomial equations (or equations that can be reduced to a polynomial equation) numerically, thus it also needs to see the structure of an equation and won't work with a numerical black box.

So the only candidate here is FindRoot.

---

However FindRoot isn't very appropriate here either. The methods it can use all assume that the function they're working with is a "nice and smooth one". Your function always returns integers, so it has a "step structure". The default FindRoot method would try to approximate the derivative of the function and would of course fail: the derivative is zero everywhere.

You can use Brent's method, but this isn't ideal either: FindRoot[medianLEfromQx[x*qx1] == 8, {x, 0, 2}, Method -> "Brent"]

Instead I would just plot the function and visually check the range of x which satisfies this equation.

Plot[medianLEfromQx[x qx1] - 8, {x, 0, 10}]

Answer 2

The way I understand your setup makes it seem like there's a simpler approach. The survival vector should be monotonic. So if you want the life expectancy to be, say, 8 years, then you want the 8th entry to be 1/2.

Clear[multiplier];
multiplier[le_Integer, qx_?(VectorQ[#, NumericQ] &)] := 
 Block[{x}, 
  First @ Sort @
    Select[
      x /. NSolve[Drop[FoldList[#1 (1 + #2) &, 1, - x * qx], -1][[le]] == 0.5, x, Reals],
      0 < # &]
  ]

Example

multiplier[8, qx1]
(*
  0.940737
*)

Since this is using approximate reals, there may be some boundary issues occasionally. You also might want to add some sanity checks. For instance,

Table[multiplier[x, qx1], {x, 2, 20}]
(*
  {10.4778, 5.50098, 3.3954, 2.27372, 1.62347, 1.21458, 0.940737, 
   0.748985, 0.610999, 0.508783, 0.430586, 0.36941, 0.320942, 0.28194,  
   0.249964, 0.223317, 0.200685, 0.181273, 0.164589}
*)

The first few multipliers would make the entries in qx1 greater than 1. In those cases, the function multiplier should print an error message I suppose.

---

Update 29 Dec 2014

In a comment, the OP was interested in adapting the above method to an interpolated mortality curve. Here's a way.

I'm unfamiliar with the standard way of interpolating mortality, but linear or exponential seem likely candidates. So one of these two, with the multiplier built into the InterpolatingFunction is the way to set it up:

Interpolation[
 Transpose[{Range@Length@qx1, Drop[FoldList[#1 (1 + #2) &, 1, -x*qx1], -1]}], 
 InterpolationOrder -> 1]

Exp @* Interpolation[
  Transpose[{Range@Length@qx1, Log @ Drop[FoldList[#1 (1 + #2) &, 1, -x*qx1], -1]}], 
  InterpolationOrder -> 1]
  (* Use Composition[Exp, Interpolation[<..>] in V9 or earlier *)

Now the set up I really want is that the interpolation should be a function of a vector of annual probabilities of death qx. So the definition I will use is

lxIF[qx_?(VectorQ[#1, NumericQ] &)] := 
  lxIF[qx] = 
   Evaluate[
     Exp @* Interpolation[
       Transpose[{Range@Length@qx, Log @ Drop[FoldList[#1 (1 + #2) &, 1, -#*qx], -1]}], 
       InterpolationOrder -> 1]
     ] &;
lx[x_, qx_?(VectorQ[#1, NumericQ] &)] := lxIF[qx][x];

A few things need mentioning. First,I replaced the multiplier x by the Function argument #. Second, we will be using this function many times, so used memoization to cache the interpolation in lxIF[qx] the first time it is computed so that it will be reused instead of recomputed. Finally, the function call lxIF[qx][x] replaces the argument # in lxIF[qx] by x and returns an InterpolatingFunction that is a function of life expectancy.

To calculate the probability of surviving 8 years for a multiplier x = 0.9, use

lxIF[qx1][0.9][8]
(*  0.516127  *)

To find the multiplier for the median life expectancy to be 8.5 years, use

FindRoot[lxIF[qx1][x][8.5] == 0.5, {x, 1.}]
(*  {x -> 0.833662}  *)

A general use function can be constructed thus:

multiplier2[le_?NumericQ, qx_?(VectorQ[#1, NumericQ] &)] := 
 Block[{x}, x /. FindRoot[lxIF[qx][x][le] == 0.5, {x, 1.}]]

Note that FindRoot works well here because interpolating functions have derivatives. Even though lxIF has a discontinuous derivative, it is strictly monotonic, which makes root-finding easy.

The mean life expectancy for a given multiplier x can be computed with

meanLE[x_?NumericQ, qx_?(VectorQ[#1, NumericQ] &)] := 
  NIntegrate[lxIF[qx][x][le], {le, 1, Length[qx]}];

FindRoot works on it, too.