Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question
    
Not checking the actual input to SOlve. Always check that it's what you think it is. –  Daniel Lichtblau Mar 13 '14 at 21:51

2 Answers 2

up vote 7 down vote accepted

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

enter image description here

share|improve this answer
    
The plot will easily give you an x that satsfies the equation. If you need the full precise range of x, you can either write a bisection method, or extract it from Plot. I typed this up in a rush, let me know if you need more info. –  Szabolcs Mar 13 '14 at 20:50
    
Yes, thanks. I need to apply this in an automated way, so I implemented a brute-force-ish method that works. I wonder, if questions of this sort are posted weekly, if Wolfram might change the error messages produced by these functions to be more useful, and/or write about these issues more directly in their documentation. –  Michael Stern Mar 13 '14 at 21:52
    
@MichaelStern Well, they already have the "knowledge base" article I linked to. I'm not sure it's possible at all to fix this shortcoming, it's too deeply rooted in the language. I think everybody here fell into this trap at least once while using Mathematica. I certainly did. –  Szabolcs Mar 13 '14 at 22:00
2  
@MichaelStern Also sorry if the answer sounded a bit rough, I had to type it up in a hurry and in one go. I recommend you try the Brent method if all you need is finding one single value that satisifies the equation (i.e. if it doesn't matter whether the number returned is 0.86 or 0.89 for as long as both satisfy the equation). –  Szabolcs Mar 13 '14 at 22:11
    
Brent method proved very useful. –  Michael Stern Apr 15 '14 at 21:31

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.

share|improve this answer
    
@MichaelStern I based my answer on how I understood your description of the problem. If I've made a mistake, please let me know. –  Michael E2 Mar 13 '14 at 23:50
    
that's a perfectly good method for integer Median LEs, thank you. I'm trying to think of how to generalize it for non-integer LEs and Mean LEs. –  Michael Stern Apr 15 '14 at 19:58
    
@MichaelStern From the code in your post, the median LE is given by Length and therefore is an integer. How do you generalize it to non-integer LEs? –  Michael E2 Dec 29 '14 at 6:49
    
In practice, people often interpolate a monthly mortality curve and compute their results on that basis. It is also possible to interpolate a continuous curve. –  Michael Stern Dec 29 '14 at 12:17
    
@MichaelStern Thanks. I made an update, which I hope is along the lines of you described. One question: It seems to me that the domain of the continuous interpolation should be {0, Length[qx]} -- i.e. from 0 to 30 in the case qx1. But Interpolation will only go from the end points of the data, 1 to 30. Would interpolating Transpose[{Range[0, Length@qx1], FoldList[#1 (1 + #2) &, 1, -x*qx1]}], i.e. not dropping the last survival probability and adjusting the domain to 0 to 30, be correct? (If I'm way off beam, I'll just delete this answer.) –  Michael E2 Dec 29 '14 at 17:46

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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