Repository navigation
Numerical precision of normal_lcdf #1284
Description
Activity
Thanks for reporting.
This is a known problem with pretty much every one of our
lcdfandlccdffunctions, which are implemented aslog(foo_cdf(...))andlog1m(foo_cdf(...)). What's needed are iterative algorithms for the cdfs that work on the log scale with stable derivatives. Often what happens when you get out into the tails,log(foo_cdf(...))remains finite, but one or more of its derivatives, which have the term1 / foo_cdf(...), explode.I don't think derivatives are the problem here, but computing log(1+erf(z)), using erf from boost.
R stats package pnorm.c https://github.com/SurajGupta/r-source/blob/a28e609e72ed7c47f6ddfbb86c85279a0750f0b7/src/nmath/pnorm.c has numerically more accurate case for computing lcdf in tails with trivial derivatives. both pnorm.c and Boost erf are using the same rational Chebyshev approximation, but pnorm.c handles the cdf and lcdf far in tails with more care (see lines 199-204 and 219-260). Boost has cdf for normal but not lcdf, so it seems that the solution would be the follow the approach used in pnorm.c
I tested that using
log(pnorm(x,log.p=F))in R results also to loss of precision (with same discrete levels in lcdf values, but slightly different location of discontinuities. So it really is the case that lcdf specific computation is needed (and the pnorm.c has the reference implementation)How does the phi_approx compare to this? It may not be as precise, but it should be better as it stays continuous...of course i don’t expect great precision in the far tails, but maybe that’s better than the discrete steps you see. To be clear, I mean the common logisitic approximation to the cdf which is described in the stan manuals.
log(Phi_approx(x))underflows to 0 in this range, ie, it's much worse.No...you should use the log inv logit function and define it accordingly.
You wrote phi_approx and I tested phi_approx. Can you be more specific and how your suggestion compares to the algorithm used in pnorm.c?
You should always avoid calculating on the natural scale and log at the very end if you target log.
So this is phi approx:
inv_logit(0.07056 * pow(x, 3.0) + 1.5976 * x)
And in stan we have a log inv logit...so you should try
Log_inv_logit(0.07056 * pow(x, 3.0) + 1.5976 * x)
I hope this helps.
I know about "avoid calculating on the natural scale and log at the very end if you target log." and that's why I recommend the approach from pnorm.c which does things in log scale and has been demonstrated to produce good results.
So this is phi approx:
The confusion came as there is also a function called Phi_approx().
I know this approximation, but I doubt this is not good approximation for normal_lcdf compared to what pnorm.c does. I think that in this issue, we should focus on accurate computation for normal_lcdf and continue in the original discourse thread discussion about whether @alecksphillips could use this approximation instead of normal_lcdf.Phi_approx often works because often there is no difference even between logit or probit (Phi_approx has x^3 term so it's not just logit).
Sure...the issue should be on the pnorm.c implementation and my note is better made in discourse. For practical purposes the log inv logit solution may give good enough results. If the approximation gives already good fits then this issue can be prioritized adequatley.
Log_inv_logit(0.07056 * pow(x, 3.0) + 1.5976 * x)
Relative error compared to solution in pnorm.c in the example range is about 4550000%
Sure...but this is far in the tails! Does it really matter?
I think it matters for this issue, and it's separate thing whether it matters in a certain model, which should be discussed somewhere else.
lcdf specific computation is needed (and the pnorm.c has the reference implementation)
This is the way forward. We need an algorithm for the derivatives, too, if differentiating the algorithm doesn't work. We can't directly import GPL-ed code---is there a reference to the algorithm for this and others like it?
The other thing to worry about is the complementary CDF. That probably needs a different algorithm to retain precision.
The references for the algorithm are mentioned in pnorm.c
DESCRIPTION
The main computation evaluates near-minimax approximations derived
from those in "Rational Chebyshev approximations for the error
function" by W. J. Cody, Math. Comp., 1969, 631-637. This
transportable program uses rational functions that theoretically
approximate the normal distribution function to at least 18
significant decimal digits. The accuracy achieved depends on the
arithmetic system, the compiler, the intrinsic functions, and
proper selection of the machine-dependent constants.REFERENCE
Cody, W. D. (1993).
ALGORITHM 715: SPECFUN - A Portable FORTRAN Package of
Special Function Routines and Test Drivers".
ACM Transactions on Mathematical Software. 19, 22-32.EXTENSIONS
The "_both" , lower, upper, and log_p variants were added by
Martin Maechler, Jan.2000;
as well as log1p() and similar improvements later on.That 1969 paper has the rational Chebyshev approximation. SPECFUN code https://www.netlib.org/specfun/ doesn't mention any license
I couldn't find that extras by Martin Maechler would be reported in any article or report.
upper is same as the complementary CDF
upper and log variants don't seem to have any article reference
boost seems to have the same Chebyshev approximations, but doesn't have upper and log variants.
I checked that for the range z>-20 (z=(x-mu)/sigma) we can use the simple approach used in log_ndtr in scipy https://github.com/scipy/scipy/blob/9109a70aad860f19bb2eeb7147d61efebedfb7f6/scipy/special/cephes/ndtr.c
This would be easy to add and scipy license is BSD3. It would be possible to further extend to z<=-20, but that requires a bit more work.For the range z<6 we would use what we have now, and for z>6 we would have code equivalent to -(exp(normal_lcdf(-x | 0, 1)) which is based on the approximation log(1+x) \approx x. Naturally the exact code is not that, but I hope that gives the idea of the computation. I tested that even direct -(exp(normal_lcdf(-x | 0, 1)) produces the similar accuracy as pnorm.
I've been working with @alecksphillips on a new implementation of normal_lcdf(x | mu, sigma) to solve this issue. I've been using the Abramowitz and Stegun numerical approximations that are referenced in several sources but the new algorithm still suffers from the same issue. However, I did manage to derive some simple equations for the derivative. The only issue is it goes to infinity as x->0+ so I've had to bridge the gap with Taylor polynomials.
It would be good to understand the R solution as it seems stable well beyond z<-20. I'll have a look at the approach used in scipy too.
Based on the error rates given for Abramowitz and Stegun approximations in https://en.wikipedia.org/wiki/Error_function#Numerical_approximations, they are not as good as the Rational Chebyshev approximations by Cody (see the references above). Boost has error rates for erf implementation here https://www.boost.org/doc/libs/1_64_0/libs/math/doc/html/math_toolkit/sf_erf/error_function.html. I think boost erf is good choice and there would not be need to change erf approximation. We would need to only adjust for log(erf) in z>6 as in scipy and if you want to add something for z<-20, something simpler than what scipy uses would be good (scipy uses iterative algorithm which has error rate based stopping rule so that would be most accurate, but I've seen several non-iterative proposals which might be sufficient).
Hi All,
We eventually came up with a solution and were able to implement it in user-defined C++ files (attached), which has served its purpose in our application.
Most of the inaccuracy was resolved by switching from using erf(x) to erfc(x), since erfc(x)->0 as x -> +infinity (allowing the 0's to go into the exponent, thereby preserving numerical precision). For the derivative function I couldn't find a numerical approximation that gave decent accuracy for the region 0.1<x<2.2. I've therefore used 3 Taylor expansions of the analytic expression to bridge this gap.
I was now wondering if anyone could advise on the best procedure to implement this function in the stan math library?
Thanks! The place to start is here:
https://github.com/stan-dev/stan/wiki/Contributing-New-Functions-to-Stan
It's a lot to digest, so please don't hesitate to ask question on our forums---there's a Developer category:
- added a commit that references this issue
on Nov 20, 2019
Description
The normal_lcdf(x | mu, sigma) function suffers from precision problems at x >> mu/sigma when compared to the equivalent pnorm(x, mu, sigma, log.p=T) in base R.
Example
R code:
Expected Output:
See plot comparing Stan and R's output here: https://discourse.mc-stan.org/t/numerical-precision-of-normal-lcdf/9685
Current Version:
v2.19.1