4

Question

The Emacs Calc package info page for Probability Distribution Functions has this to say about the normal distribution function:

The ‘utpn(x,m,s)’ function uses a normal (Gaussian) distribution with mean ‘m’ and standard deviation ‘s’. It is the probability that such a normal-distributed random variable would exceed ‘x’.

I am interested in solving for x in utpn(x,m,s) = p given p, m and s. The same info page has this to say about that:

While Calc does not provide inverses of the probability distribution functions, the a R command can be used to solve for the inverse. Since the distribution functions are monotonic, a R is guaranteed to be able to find a solution given any initial guess.

My question: how can I "solve for the inverse" as described above?

Update

I have provided my own answer to this question with many details. Credit to NickD's answer and db48x's answer for pointing me in the right direction.

This has also become Emacs bug#58929.

What I have tried

I have tried entering the formula utpn(x,10,1)=0.5 and invoking a R with an initial guess of 10. In this case I want the x that would cause utpn(x,10,1) to return 0.5. Note that the correct answer here is x=10 since I am passing m=10 for the "mean" arg to utpn, and by definition a random value from the normal distribution has an 0.5 chance of being larger than the mean.

I type:

` utpn(x,10,1)=0.5 RET
10 RET
a R RET
x RET

and I end up with this in the *Calculator* buffer:

2:  utpn(x, 10, 1) = 0.5
1:  [10., 0.]

This is correct. The 'a R' command has returned a vector whose first element is x and the second is the function's value at that x.

However, notice that I have given the answer in the initial guess. The problem is that the function doesn't do well when given other guesses.

` utpn(x,10,1)=0.5 RET
1.0 RET                  <---- a different guess
a R RET
x RET

results in this:

2:  utpn(x, 10, 1) = 0.5
1:  root(utpn(x, 10, 1) = 0.5, x, 1)

With this in the minibuffer: Newton's method failed to converge: 4.8639203627839e17

I'm finding that if the initial guess is close enough to the right answer, in this case 10, then a R works, otherwise not. E.g. guessing 8.6 fails but 8.7 works.

Since the distribution functions are monotonic, a R is guaranteed to be able to find a solution given any initial guess.

Have I found a bug here, or am I doing something wrong?

3 Answers 3

3

I think the manual is a bit too optimistic: these curves may be monotonic, which in theory means you can solve for any value of p between 0 and 1, but if you plot them, they approach 0 (or 1) asymptotically and while doing so they are VERY flat (almost horizontal). That's exactly the situation where Newton's method (which a R uses) becomes useless unless you choose a good starting point. So although the manual is so romantically optimistic, real life intrudes and dashes such hopes.

Here's what the curve looks like (obtained with Calc: [0..20] 'utpn(x, 10, 2) RET g f):

enter image description here

You need to choose a starting point as far away from the tails as possible: anywhere within a standard deviation from the solution would be OK, but if you are outside by a few standard deviations, Newton's method is going to at best converge very slowly.

For utpn(x, 10, 2) - 0.5 (which Calc implicitly equates to 0 in order to solve it), I get convergence for initial values between 7.23 and 12.79 (there is a slight asymmetry between left and right offsets from the mean, which I'm blaming on accidents of implementation of the algorithm and the finite arithmetic that any computer does, but I don't know for sure):

'utpn(x, 10, 2) - 0.5
8
a R x RET

shows [10.0, 0.0] as the result, but if I do 7.2 instead of 8, it fails to converge, i.e. Calc gives up before it has a chance to converge, although given round-off error, it might or might not converge (or even converge on a wrong answer), even if the number of iterations is increased.

BTW, although in theory, you could solve for p very close to 1 or very close to 0, the shape of the graph is going to cause almost insurmountable problems, so if possible, you should avoid those areas. If it is not possible, Newton is not going to help and you'll have to find another method and probably be satisfied with less accuracy.

I recommend that you consult a good Numerical Analysis book. My favorite is "Numerical Methods that (usually) work" by Forman S. Acton. It's an old (1970) book, but you can read a relatively recent (2015)review of it here.

4
  • Yes, that seems to be the "root" of the problem (pun intended). Incidentally, I find utpt to be even more likely to tickle Calc's root function into not working. Yet, I wrote a bisection method in Emacs Lisp that solves this easily. I now find a large class of algorithms described at en.wikipedia.org/wiki/Root-finding_algorithms -- and it looks like a relative newcomer may have advantages worth considering: en.wikipedia.org/wiki/ITP_method Commented Oct 31, 2022 at 23:42
  • 1
    A lot of things can go wrong with bisection or any other method as well (consult chapters 2, 7 and 14 of Acton's book); and where Newton works, it converges faster (but when it fails, it fails spectacularly). The best advice I think is to know as much about your problem before you try to solve it. Graphing the relevant function (or functions) should probably be your first step (but high dimensionality may cause problems). Acton says (in another book): "Sketching an equation before choosing an algorithm will prevent more computational disasters than any other effort."
    – NickD
    Commented Nov 1, 2022 at 1:59
  • 1
    Note that both answers (@db48x beat me to the punch - darnit !-) ) start by graphing the function. That's not a coincidence! Understanding the geometry is the first step towards a solution.
    – NickD
    Commented Nov 1, 2022 at 2:17
  • Thanks for your answer! In looking in to how to improve Calc to gracefully handle this case I found that it already had heuristics that help, which I have detailed in my own answer here: emacs.stackexchange.com/a/74383/21036. Separately, I think in the 50 years since Acton's book was published we can expect better of our computer assisted math packages. His approaches made a lot of sense when computers had the power of a 1990's era Casio watch, but today I think we can expect the computer to "look at the graph" and pick a good heuristic automatically, especially for simple cases. Commented Nov 1, 2022 at 18:34
3

Graph that function; it doesn’t have much slope at 1:

[0 .. 20]
' utpn(x,10,1) RET
g a C-u g p

enter image description here

Newton’s Method uses the slope to determine what guess to try next, so I guess it’s not that surprising that it would fail in this case. The number in the error message is the next guess that Newton’s method wants to make, based on the derivative at the previous guess. Since you know that the answer is 10, you can see how far off it has gotten.

3
  • 1
    I was a bit disappointed that you beat me to the punch :-) but c'est la vie - I upvoted your answer (which deserves to be accepted - hint to the OP). The only thing I'm not sure about (I haven't checked the Calc sources) is the error message: are you sure it's the value of the derivative? It's not close to 0 at all, it's huge (10^17). I think it's more likely the next value of x in the iteration - right? Basically, the almost horizontal tangent meets the x-axis wayyyy to the right, blowing up the computation.
    – NickD
    Commented Nov 1, 2022 at 2:13
  • 1
    Oh, you’re right. I was misreading the code; this error happens when the next value is not less than the limit value.
    – db48x
    Commented Nov 1, 2022 at 8:15
  • Thanks for this answer. It got me on the right track to a solution, which I've posted as an answer here: emacs.stackexchange.com/a/74383/21036 Commented Nov 1, 2022 at 18:29
1

TLDR

Give calc's a R (calc-find-root) [root] function a range, instead of a single numerical guess. Make sure the function evaluates to different signs on each end of the range. This gives calc enough information to fall back to a different root finding algorithm that will succeed.

For example, here I am attempting to solve utpn(x, 10, 1) = 0.5 for x, guessing that it might be x = 1.0.

` utpn(x,10,1)=0.5 RET
1.0 RET                  <---- the guess
a R RET
x RET

The above doesn't work, producing the error message "Newton's method failed to converge: 4.8639203627839e17".

Calc succeeds if the guess is the range [0..1000]:

` utpn(x, 10, 1) = 0.5 RET
[0..1000] RET                  <---- the better guess
a R RET
x RET

Success! Calc produces this:

2:  utpn(x, 10, 1) = 0.5
1:  [10., 0.]

Background

As pointed out in both NickD's answer and dbx48's answer, calc is having difficulty due to the extremely flat portions of the utpn function. Let's graph utpn(x, 10, 1) for x over the range [0..20]:

[0..200]
10
/
' utpn(x, 10, 1)
g f

utpn(x, 10, 1) across the range 0..201

You can see that outside the range of about [7..13] the function is very flat.

The problem with Newton's algorithm

Clues into Calc's behavior are found in its Root Finding documentation.

The initial guess can be a real number, in which case Calc searches for a real solution near that number, or a complex number, in which case Calc searches the whole complex plane near that number for a solution, or it can be an interval form which restricts the search to real numbers inside that interval.

Recall that in my original question I was giving the a R function a single guess:

` utpn(x, 10, 1) = 0.5 RET
1.0 RET                      <---- the guess
a R RET
x RET

On this the documentation says:

Calc tries to use a d to take the derivative of the equation. If this succeeds, it uses Newton’s method. If the equation is not differentiable Calc uses a bisection method.

So, let's see if utpn(x, 10, 1) - 0.5 is differentiable.

` utpn(x, 10, 1) - 0.5 RET
a d RET
x RET

Calc presents this derivative:

1:  0.398942280401 / exp((7.07106781188 - x / 1.41421356237)^2)

Let's plot that derivative over a range of [0 .. 20]:

[0..200]
10
/
' 0.398942280401 / exp((7.07106781188 - x / 1.41421356237)^2)
g f

Graph of derivative of utpn(x, 10, 1) over the range 0..20

The derivative is nearly zero outside of the range [7 .. 13]. Numerically, the derivative approaches zero asymptotically, because utpn itself asymptotically approaches zero.

Because calc thinks it can compute the derivative it will use (Newton's method)[https://en.wikipedia.org/wiki/Newton%27s_method) to find the root, and here things go astray.

Why? There is nothing mathematically wrong with using Newton's method with utpn. My guess: calc's limited precision arithmetic are exceeded during the computation, which "confuses" the algorithm. One intuition: Newton's method computes the next guess for x as x - f(x)/f'(x). Above calc produced an error message "Newton's method failed to converge: 4.8639203627839e17". Let's see the values involved at that point:

  1. utpn(4.8639203627839e17, 10, 1) - 0.5 evaluates to -0.5
  2. 0.398942280401 / exp((7.07106781188 - 4.8639203627839e17 / 1.41421356237)^2) yields an error: "Floating-point overflow occurred: exp(1.18288606478e35)"

I have no explanation for how Newton's algorithm has gone from the initial guess of 1.0 all the way out to 4.8639203627839e17 without realizing the root was at 10.0. Presumably calc's limited precision arithmetic is to blame for that, too.

Calc's fallback algorithm: Bisection

Calc's documentation claims it has code to detect when Newton's method is failing:

(If Newton’s method appears to be going astray, Calc switches over to bisection if it can, or otherwise gives up. In this case it may help to try again with a slightly different initial guess.)

For some reason this is not working for the particular case discussed here.

Calc does have an alternate algorithm, which uses a bisection approach. Read up on Bisection method on Wikipedia. Because the method does not use the derivative it is guaranteed to find a root if it has an interval where f(x) is of different sign at each end. The calc manual discusses this:

If the formula (or the difference between the sides of an equation) is negative at one end of the interval you specify and positive at the other end, the root finder is guaranteed to find a root. Otherwise, Calc subdivides the interval into small parts looking for positive and negative values to bracket the root. When your guess is an interval, Calc will not look outside that interval for a root.

Helping Calc use the Bisection method

Let's try the original problem giving an interval guess instead of a single number:

` utpn(x, 10, 1) = 0.5 RET
[0..1000] RET
a R RET
x RET

Success! Calc produces this:

2:  utpn(x, 10, 1) = 0.5
1:  [10., 0.]

Could wroot work too?

Calc has another root solving function, wroot, described in the manual as follows:

The H a R [wroot] command is similar to a R, except that if the initial guess is an interval for which the function has the same sign at both ends, then rather than subdividing the interval Calc attempts to widen it to enclose a root. Use this mode if you are not sure if the function has a root in your interval.

How might this fare with the original problematic case:

` utpn(x, 10, 1) = 0.5 RET
1.0 RET
H a R RET
x RET

No luck:

2:  utpn(x, 10, 1) = 0.5
1:  wroot(utpn(x, 10, 1) = 0.5, x, 1.)

With the message "Newton's method failed to converge: 4.8639203627839e17".

Presumably, what is going on here is that the extremely flat derivative causes the root finding algorithm to quickly reach very large x values, where the derivative becomes unstable.

Conclusion

The flat derivative of the utpn function is a worst case scenario for the Newton's method algorithm used by calc, given calc's limited precision arithmetic. Calc can sometimes detect this and switch to a (Bisection method)[https://en.wikipedia.org/wiki/Bisection_method], but in this case that heuristic is failing. To work around this, supply an explicit range to the "a R (calc-find-root) [root]" function, such that the sign of the result is different for each end of the interval, and the Bisection method should kick in and work.

4
  • In the Helping Calc use the bisection method, you are quoting the doc for H a R: ...except that if the initial guess is an interval for which the function has the same sign at both ends..., but you didn't give such an interval: you gave it a single value. If you give it an interval, it works.
    – NickD
    Commented Nov 1, 2022 at 19:42
  • BTW, Calc does arbitrary precision integer (and rational) arithmetic , but the floating point precision is limited (although far, far bigger than with typical programming languages). See C-h i g(calc)data types RET. If you can afford the time, you can bump up the precision and see how that works :-)
    – NickD
    Commented Nov 1, 2022 at 19:55
  • @NickD, yes, you're right about "H a R" being guaranteed to work with an interval with differing signs at each end. I had already demonstrated this for "a R". The point of trying "H a R" was to see if it could solve it given a single number. I've made this more clear by putting the "H a R" stuff in its own section. I have also tried this out to a precision of 1000...no real change in results. Rather hilariously, it tells me newton's method fails to convert and prints a 1000 digit number beginning with 486392036277920968.20377... Commented Nov 2, 2022 at 22:32
  • 1
    @NickD oh, I see what you mean. Even if I give utpn(x, 10, 1) - 0.5 a guess of [0..1] it correctly computes the solution at 10.0. Nice! Commented Nov 2, 2022 at 22:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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