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 at 21:51
add comment

2 Answers 2

up vote 4 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 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 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 at 22:00
1  
@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 at 22:11
    
Brent method proved very useful. –  Michael Stern Apr 15 at 21:31
add comment

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.

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 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 at 19:58
add comment

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.