Bugfix/issue #1284 numerical precision of normal lcdf by PhilClemson · Pull Request #1411 · stan-dev/math · GitHub

Summary

Wrote a new piecewise numerical approximation of the normal_lcdf function and its derivatives. This is based on algorithms derived from W. J. Cody, Math. Comp. 23(107):631-638 (1969) and Abramowitx and Stegun (1964), in addition to original Taylor expansions that have been derived to bridge the gap where the numerical approximations are unstable.

Tests

Two new valid_values tests have been added to test/prob/normal_cdf_log_test.hpp to test the function at values far from mu/sigma, which fail for the old implementation and pass in the new implementation.

I also had to alter other tests which were based on the limits of the old function (expected -infinity in regions where actual approximations of the the function are now defined). The new limits are based on the maximum values for a double. In addition, some cases tested against a seemingly arbitrary absolute error (e.g. <10^-4), rather than a relative error. This makes sense when the function converges asymptotically or when the absolute value is small. However, in the case of the derivative of the log function this does not converge but increases fairly linearly as y -> -infinity. As such it is expected that the absolute error will be much larger than 10^-4 for these values, so the tests now check for the relative error before considering the absolute error test.

Side Effects

None that I am aware of.

Checklist

Math issue Numerical precision of normal_lcdf #1284
Copyright holder: University of Liverpool

The copyright holder is typically you or your assignee, such as a university or company. By submitting this pull request, the copyright holder is agreeing to the license the submitted work under the following licenses:
- Code: BSD 3-clause (https://opensource.org/licenses/BSD-3-Clause)
- Documentation: CC-BY 4.0 (https://creativecommons.org/licenses/by/4.0/)
the basic tests are passing
- unit tests pass (to run, use: ./runTests.py test/unit)
- header checks pass, (make test-headers)
- docs build, (make doxygen)
- code passes the built in C++ standards checks (make cpplint)
the code is written in idiomatic C++ and changes are documented in the doxygen
the new changes are tested

…sue (causes test failures).

…stable/2017-11-14)

Thanks for submitting. The unit tests for distributions are failing and Jenkins is pointing to a .zip file. @serban-nicusor-toptal do you know how to find that file?

Would it be better to factor this into a separate log_Phi() function? There are several other functions that could benefit from this.

normal_lccdf
skew_normal_lpdf
exp_mod_normal_lpdf
exp_mod_normal_lcdf
exp_mod_normal_lccdf
ordered_probit_lpdf
I guess only normal_lccdf is in scope of this PR, others should have their own issues.

@bob-carpenter Hey, the job ran on-demand so the machine is already gone. I will try to restart this job now

Would it be better to factor this into a separate log_Phi() function? There are several other functions that could benefit from this.
* `normal_lccdf`

* `skew_normal_lpdf`

* `exp_mod_normal_lpdf`

* `exp_mod_normal_lcdf`

* `exp_mod_normal_lccdf`

* `ordered_probit_lpdf`
  I guess only `normal_lccdf` is in scope of this PR, others should have their own issues.

Yes normal_lccdf could benefit from the same approximations (assumed it was better to do one change at a time). I believe the other functions would require independent numerical approximations though.

Any step in the right direction is appreciated. A log_Phi function would definitely be helpful. As far as I know, all of our log cdfs are just implemented as logs of the cdf. The derivative of log F(x) is F'(x) / F(x), so if the CDF F(x) and its derivative F'(x) are not stable, it's hard to make the log cdf stable for derivatives.

Thanks for working through all this Phil! I had a couple comments about how things are arranged in the code.

I think my big question is how do these piecewise definitions affect derivatives? Like two functions join together at zero. Do the gradients of the functions on the right and the left match? And do they match the value of the analytical gradient given to ops partials?

Also what about higher order derivatives? Do we have 2nd derivatives through all our piecewise things? 3rd order I think is the highest we'd ever likely go -- any idea on those? Feel free to ask questions back. I'd have to think a bit how the third order stuff would work.

I understand the last implementation was fragile and would blow up a lot, so it's fine that this will have limitations too (we just want this to blow up less). I just want to be clear on the limitations before I do the deep dive on the tests.

In the same way as above would it be possible to flatten these conditionals on the piece-wise defined derivatives? It looks like they're all defined in terms of scaled_diff.

Also would it be possible to write them in order of what interval they represent using similar notation? Like it's difficult to parse sometimes scaled_diff > 1.5 sometimes scaled_diff < 2.2 or whatever. If it's necessary to write it that way it's fine just checking :D.

Yep that makes a lot more sense. I think the way to do this would be

if (scaled_diff > 2.2) {
...
} else if (scaled_diff > 1.5) {
...
} else if (scaled_diff > 0.8) {
...
} else if (scaled_diff > 0.1) {
...
} else if (scaled_diff > -4.0) {
...
} else {
...
}

although there will be some code duplication in the 4th and 5th statements (but I'm guessing this is ok if it improves code readability?).

Did a teeny review but I'm not to familiar with how to review this sort of thing so I'll leave the main review to Ben

For these we prefer if you use the is_* methods available in Stan ala

https://github.com/stan-dev/math/blob/develop/stan/math/prim/scal/err/is_not_nan.hpp

Ah, sorry wasn't aware!

We prefer if declarations happen on different lines

T_partials_return t = 0.0;
T_partials_return t2 = 0.0;
T_partials_return t4 = 0.0;

Are the t values here only used within the if scope? The compiler may have an easier time optimizing things if it knows these temps lifetimes exist only within the if scope

Nice spot, yes they should be moved to within the if braces.

(if below is a non-issue feel free to close it)

Should the constants here all have the same number of decimals? It feels like if we are going after things that are over/underflowing and near zero that could matter. Though this is out of my scope so feel free to ignore

Yep this is a valid comment. I derived these constants in Matlab which automatically rounded them to 5 sf.

Can we remove the dbl from the name? Unless there is a dncdf_log_something

Yes we can rename it. This is a residual from the original implementation where I had the derivative defined in a separated function called "dncdf_log".

You can use return_type_t here instead of typename stan::return_type<...>::type

Okay, didn't know that! :)

No worries! I'm working on a PR right now that should make these sorts of things easier to find on the math docs site

Am I ok to do a new commit with these changes?

Ah, sorry wasn't aware!

Yep that makes a lot more sense. I think the way to do this would be

if (scaled_diff > 2.2) {
...
} else if (scaled_diff > 1.5) {
...
} else if (scaled_diff > 0.8) {
...
} else if (scaled_diff > 0.1) {
...
} else if (scaled_diff > -4.0) {
...
} else {
...
}

although there will be some code duplication in the 4th and 5th statements (but I'm guessing this is ok if it improves code readability?).

Nice spot, yes they should be moved to within the if braces.

Yep this is a valid comment. I derived these constants in Matlab which automatically rounded them to 5 sf.

Yes we can rename it. This is a residual from the original implementation where I had the derivative defined in a separated function called "dncdf_log".

Okay, didn't know that! :)

Thanks for working through all this Phil! I had a couple comments about how things are arranged in the code.

I think my big question is how do these piecewise definitions affect derivatives? Like two functions join together at zero. Do the gradients of the functions on the right and the left match? And do they match the value of the analytical gradient given to ops partials?

Also what about higher order derivatives? Do we have 2nd derivatives through all our piecewise things? 3rd order I think is the highest we'd ever likely go -- any idea on those? Feel free to ask questions back. I'd have to think a bit how the third order stuff would work.

I understand the last implementation was fragile and would blow up a lot, so it's fine that this will have limitations too (we just want this to blow up less). I just want to be clear on the limitations before I do the deep dive on the tests.

The log cdf function is pretty tricky as the only useful symmetry you can rely on is erf(x)=-erf(-x). In general I made sure the error at the joins was no more than ~10^-6.

I checked the respective functions in R and scipy and both relied on similar piecewise approximations. However, neither of these had implemented a numerical approximation of the derivatives. The higher-order derivatives could also be defined but I imagine the approximations would look even more horrible. Do we have an idea of where they might be applied?

Am I ok to do a new commit with these changes?

Yeah, definitely. How this works is after we post reviews, you patch things until you're happy with it and then poke the reviewer again (me in this case) and I'll come back and have a look.

I checked the respective functions in R and scipy and both relied on similar piecewise approximations.

Yeah I figure there's probably a lot of nastiness under the hood in a lot of these complex functions. I suppose there's probably a lot of these things that are just hidden from us in standard libraries and stuff.

The higher-order derivatives could also be defined but I imagine the approximations would look even more horrible.

Since this is a one argument function, let's pull in Bob's autodiff tester and see what it tells us for the derivatives. Instructions here: https://github.com/stan-dev/math/wiki/Automatic-Differentiation-Testing-Framework, examples here: https://github.com/stan-dev/math/blob/develop/test/unit/math/mix/scal/fun/bessel_first_kind_test.cpp

It should be pretty easy to add tests with this. Do you mind just playing around and seeing what you find? Like run the tests in the tails, at the joining points in values, at the joining points in derivatives, etc. I think if we find problems out in the super-far tails it's fine. We just want to convince ourselves things worked at least where they worked before (and hopefully farther out too).

Ask if it's not clear how to use that tester. It's some fancy C++.

Bob's autodiff test thing basically assumes we're set up good prim tests to us (evaluating precision and checking edge cases). Once we have those it tries to automate all the gradients. Ofc. in this situation the finite differences and stuff might just break down before the new numerics but we can at least try :D.

... The higher-order derivatives could also be defined but I imagine the approximations would look even more horrible. Do we have an idea of where they might be applied?

The gradients in reverse-mode get used in all of our current algorithms (HMC, L-BFGS, and ADVI). We want to use higher-order autodiff to calculate Hessians for higher-order algorithms and for Laplace approximations and as nested components of our ODE and other solvers that need Jacobians at various points which will be more efficient in forward mode.

It should be pretty easy to add tests with this. Do you mind just playing around and seeing what you find? Like run the tests in the tails, at the joining points in values, at the joining points in derivatives, etc. I think if we find problems out in the super-far tails it's fine. We just want to convince ourselves things worked at least where they worked before (and hopefully farther out too).

Sure I'll definitely have a go at this. Completely understand the need to check these things.

The gradients in reverse-mode get used in all of our current algorithms (HMC, L-BFGS, and ADVI). We want to use higher-order autodiff to calculate Hessians for higher-order algorithms and for Laplace approximations and as nested components of our ODE and other solvers that need Jacobians at various points which will be more efficient in forward mode.

Ah okay, just wanted to be sure there was a useful application before we delve down the rabbit hole!

…mprovements to the numerical approximations.

https://github.com/PhilClemson/math into bugfix/issue-#1284-Numerical-precision-of-normal_lcdf

I've added in all of the changes in the latest commit. After playing around with the autodiff tester I found some regions where it was failing some of the tests, so I added some new approximations (Taylor expansions and fits of the residuals) to improve the accuracy.

The only remaining issue is happening right on the boundaries in the piecewise function. It passes all but one of the autodiff tests, here's an example of the output:

[----------] 1 test from mathMixScalFun
[ RUN ] mathMixScalFun.lcdf_derivatives
./test/unit/math/expect_near_rel.hpp:50: Failure
The difference between 0 and relative_diff is 2, which exceeds tol, where
0 evaluates to 0,
relative_diff evaluates to -2, and
tol evaluates to 0.01.
expect_near_rel_finite(-2.7823299091694418, 0.0031959719798173507, tolerance = 0.01); relative diff = -2
in: expect_near_rel; require items x1(i) = x2(i): grad hessian grad_H_fd[i] == grad_H_ad[i]
[ FAILED ] mathMixScalFun.lcdf_derivatives (0 ms)
[----------] 1 test from mathMixScalFun (0 ms total)

I get relative diff = -2 in each case, but I'm not sure what it means?

Thanks for doing this! Did you add the mix tests to the branch? I only see test code for rev and prim.

Let me know when you have all the the tests up and I'll go about checking things and see what I can make of it. @bob-carpenter without looking at the code, does that output above (#1411 (comment)) indicate any specific sort of failure to you?

Thanks for doing this! Did you add the mix tests to the branch? I only see test code for rev and prim.

No sorry, I'd added the tests to the prim test file. I've now included the autodiff tests in a new mix test file.

./runTests.py test/unit/math/mix/scal/prob/normal_cdf_log_test

Runs for me (and passes). What's an example of one of the tests that was failing?

I left out the tests it was failing on to avoid confusion. It will fail when using the unscaled value of one of the borders in the function for the derivative. For example:

stan::test::expect_ad(f(0.0, 1.0), 0.1 * stan::math::SQRT_2);

I doubt it was passing those tests before. Usually that means the finite diff is failing. When the difference between numbers is big, the relative difference approaches 2. rel_diff(a, b) = abs(a - b) / (0.5 * (abs(a) + abs(b))) Suppose without loss of generality that a >> b >> 0. Then it's easy to see the numerator goes to a and the denominator to 0.5 * a and the rel diff to 2. Given that this is happening at third-order derivatives, I'd suggest just ignoring it and setting tolerances this way: stan::test::ad_tolerances tols; tols.grad_hessian_grad_hessian_ = 1e1; ... expect_ad(tols, ...);

…

On Nov 14, 2019, at 7:19 AM, PhilClemson ***@***.***> wrote: I've added in all of the changes in the latest commit. After playing around with the autodiff tester I found some regions where it was failing some of the tests, so I added some new approximations (Taylor expansions and fits of the residuals) to improve the accuracy. The only remaining issue is happening right on the boundaries in the piecewise function. It passes all but one of the autodiff tests, here's an example of the output: [----------] 1 test from mathMixScalFun [ RUN ] mathMixScalFun.lcdf_derivatives ./test/unit/math/expect_near_rel.hpp:50: Failure The difference between 0 and relative_diff is 2, which exceeds tol, where 0 evaluates to 0, relative_diff evaluates to -2, and tol evaluates to 0.01. expect_near_rel_finite(-1.8418968633584807, 0.16197045242952307, tolerance = 0.01); relative diff = -2 in: expect_near_rel; require items x1(i) = x2(i): grad hessian grad_H_fd[i] == grad_H_ad[i] [ FAILED ] mathMixScalFun.lcdf_derivatives (1 ms) [----------] 1 test from mathMixScalFun (1 ms total) I get relative diff = -2 in each case, but I'm not sure what it means? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

Cool. Thanks Bob. @PhilClemson I'll dig through this tomorrow and get back to you.

(stat_comp_benchmarks/benchmarks/gp_pois_regr/gp_pois_regr.stan, 0.99)
(stat_comp_benchmarks/benchmarks/low_dim_corr_gauss/low_dim_corr_gauss.stan, 0.96)
(stat_comp_benchmarks/benchmarks/irt_2pl/irt_2pl.stan, 1.0)
(stat_comp_benchmarks/benchmarks/pkpd/one_comp_mm_elim_abs.stan, 1.02)
(stat_comp_benchmarks/benchmarks/eight_schools/eight_schools.stan, 0.97)
(stat_comp_benchmarks/benchmarks/gp_regr/gp_regr.stan, 1.0)
(stat_comp_benchmarks/benchmarks/arK/arK.stan, 0.75)
(performance.compilation, 1.03)
(stat_comp_benchmarks/benchmarks/low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan, 1.0)
(stat_comp_benchmarks/benchmarks/low_dim_gauss_mix/low_dim_gauss_mix.stan, 1.01)
(stat_comp_benchmarks/benchmarks/sir/sir.stan, 0.95)
(stat_comp_benchmarks/benchmarks/pkpd/sim_one_comp_mm_elim_abs.stan, 0.98)
(stat_comp_benchmarks/benchmarks/garch/garch.stan, 0.89)
(stat_comp_benchmarks/benchmarks/gp_regr/gen_gp_data.stan, 1.0)
(stat_comp_benchmarks/benchmarks/arma/arma.stan, 1.0)
Result: 0.97019854218
Commit hash: c329aca

Code looks good! I tried to break this and it seems really tough (which I think was the goal). I left three easy comments, and then a question regarding the test framework in test/prob. Apologies for asking for the mix tests -- looks like this test framework does similar things lol. Keep the new tests though. I like them and I assume we'll move in that direction anyway.

Oh yeah is there a particular thing you need a CDF to be super robust for? Like is this for truncated distributions? Or for a link function? Might be worth throwing a comment in there saying something like these expansions look complicated but they're here because of X.

@PhilClemson question -- @nhuurre pointed out this code might translate to other functions. This isn't something that needs done in this pull req, but do any of these expansions generalize? Clearly it's a lot of work to get them right and so it'd be nice if it could be boxed up at all.

@bob-carpenter does your autodiff testing stuff replace the distribution autodiff tests in test/prob? The mix tests are here: https://github.com/stan-dev/math/blob/c329aca4247ff1c1d40b96d5530ab18e568b189a/test/unit/math/mix/scal/prob/normal_cdf_log_test.cpp . Also for Bob, why are there separate implementations of every CDF in the tests (cdf_log vs. cdf_log_function): https://github.com/stan-dev/math/blob/c329aca4247ff1c1d40b96d5530ab18e568b189a/test/prob/normal/normal_cdf_log_test.hpp? Is that just an artifact of history or is there something important about that?

I understand the reasoning for doing relative error checks away from zero and absolute error checks near zero. I'm not familiar with this log asymptotic thing though. Is it related?

Also I get our tolerances are bad here, but 10% off for the relative error seems high? Was that just an arbitrary choice or does it need to stay like that?

I'm not super concerned about these tests, but this changes touches other code yada yada so gotta be careful. I trust Bob's autodiff stuff and I'm pretty sure it's covering these tests (didn't realize this finite difference stuff was being checked already!). I'll ask Bob about this separately in this thread.

The asymptotic thing is just from the fact that as a function asymptotically converges towards 0, the log of that function goes to -inf. When you've got a function that converges asymptotically then the absolue error is (hopefully) going to shrink the closer it gets to that value. For the other functions it's fine to use a single absolue error check over the whole function, but for log(x) functions you have to switch to a relative error check as x->0.

10% was a completely arbitrary choice. I've changed it ot 1% now (just checked and it's still passing all the tests).

I guess these are the same sorts of tests as above. Just dropping a comment here so I remember them.

I've changed this to a 1% relative error check as well

I'm familiar with how relative differences near zero fail and we switch to absolute errors, but I don't follow the log of functions thing here. Is it similar? Is this something special with gradients?

Also 10% seems like a really coarse error? I know things in this test framework are pretty rough but 0.1 seems really large.

Oh yeah is there a particular thing you need a CDF to be super robust for? Like is this for truncated distributions? Or for a link function? Might be worth throwing a comment in there saying something like these expansions look complicated but they're here because of X.

Yes the reason we came across this issue is because we were investigating the case of rare base-rate logistic regression problems in stan. We're proposing an alternative link function for dealing with these problems, which is able to more closely fit the analytic posteriors.

@PhilClemson question -- @nhuurre pointed out this code might translate to other functions. This isn't something that needs done in this pull req, but do any of these expansions generalize? Clearly it's a lot of work to get them right and so it'd be nice if it could be boxed up at all.

Other than normal lccdf (which could be implemented as a simple y-axis reflection of normal lcdf) I'm not sure what else would be generalizable for use in other functions. The gradient approximations are definitely very specific and that's where most of the work has gone in.

ordered_probit_lpmf uses log(Phi(...)) directly so it can borrow this without any problem.
It's more complicated for skew_normal and exp_mod_normal. Naive substitution where they use log(erfc(...)) would result in catastrophic cancellation in the derivatives (in a subset of the region where they currently fail anyway...)

@nhuurre are there problems where order_probit suffers from numeric breakdowns like happened here? If so I guess we make a separate issue out of this.

…ng the other suggested changes.

…al approximations used.

Should be Abramowitz. :)

why are there separate implementations of every CDF in the tests (cdf_log vs. cdf_log_function)

No idea. Maybe to have a baseline against which to test.

does your autodiff testing stuff replace the distribution autodiff tests

No.

(stat_comp_benchmarks/benchmarks/gp_pois_regr/gp_pois_regr.stan, 0.99)
(stat_comp_benchmarks/benchmarks/low_dim_corr_gauss/low_dim_corr_gauss.stan, 0.98)
(stat_comp_benchmarks/benchmarks/irt_2pl/irt_2pl.stan, 1.0)
(stat_comp_benchmarks/benchmarks/pkpd/one_comp_mm_elim_abs.stan, 0.99)
(stat_comp_benchmarks/benchmarks/eight_schools/eight_schools.stan, 0.96)
(stat_comp_benchmarks/benchmarks/gp_regr/gp_regr.stan, 1.01)
(stat_comp_benchmarks/benchmarks/arK/arK.stan, 0.75)
(performance.compilation, 1.01)
(stat_comp_benchmarks/benchmarks/low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan, 1.0)
(stat_comp_benchmarks/benchmarks/low_dim_gauss_mix/low_dim_gauss_mix.stan, 1.01)
(stat_comp_benchmarks/benchmarks/sir/sir.stan, 0.95)
(stat_comp_benchmarks/benchmarks/pkpd/sim_one_comp_mm_elim_abs.stan, 1.01)
(stat_comp_benchmarks/benchmarks/garch/garch.stan, 0.9)
(stat_comp_benchmarks/benchmarks/gp_regr/gen_gp_data.stan, 0.96)
(stat_comp_benchmarks/benchmarks/arma/arma.stan, 1.0)
Result: 0.96712407902
Commit hash: 59171bf

I think revert the relative error check change to the probability test framework and remove the test that triggers it. I left details in the comments. I think it's unncessary there and we already covered it and I don't wanna go touching that old test framework.

The tests in rev are fine (though they pass without the relative checks too I think). They're self contained (and the stuff in mix covers them anyway).

I'll merge afterwards! This looks great, thanks!

Just delete this extra code and the test that triggers it (the y = -50 example). I don't wanna worry about changing this test framework.

Weird to say delete a test, but lemme just justify that by saying:

I had a closer look today and I agree this totally looks like a shortcoming in the distribution test framework, not the code here.

Also, you have the test in test/unit/math/mix/scal/prob/normal_cdf_log_test.cpp that covers this exact case with Bob's autodiff stuff.

Okay I've reverted back to the original file. I guess the relative error checks in the rev test should at least keep a record of the issue in the negative tail. Let me know if you need anything else!

…ks). Also corrected a couple of typos.

(stat_comp_benchmarks/benchmarks/gp_pois_regr/gp_pois_regr.stan, 1.0)
(stat_comp_benchmarks/benchmarks/low_dim_corr_gauss/low_dim_corr_gauss.stan, 0.96)
(stat_comp_benchmarks/benchmarks/irt_2pl/irt_2pl.stan, 1.0)
(stat_comp_benchmarks/benchmarks/pkpd/one_comp_mm_elim_abs.stan, 0.98)
(stat_comp_benchmarks/benchmarks/eight_schools/eight_schools.stan, 0.94)
(stat_comp_benchmarks/benchmarks/gp_regr/gp_regr.stan, 0.96)
(stat_comp_benchmarks/benchmarks/arK/arK.stan, 0.75)
(performance.compilation, 1.01)
(stat_comp_benchmarks/benchmarks/low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan, 1.0)
(stat_comp_benchmarks/benchmarks/low_dim_gauss_mix/low_dim_gauss_mix.stan, 1.0)
(stat_comp_benchmarks/benchmarks/sir/sir.stan, 0.97)
(stat_comp_benchmarks/benchmarks/pkpd/sim_one_comp_mm_elim_abs.stan, 0.95)
(stat_comp_benchmarks/benchmarks/garch/garch.stan, 0.89)
(stat_comp_benchmarks/benchmarks/gp_regr/gen_gp_data.stan, 1.0)
(stat_comp_benchmarks/benchmarks/arma/arma.stan, 0.98)
Result: 0.95944473214
Commit hash: b148f2b

In we go! Thanks Phil!

does your autodiff testing stuff replace the distribution autodiff tests

no.

why are there separate implementations of every CDF in the tests (cdf_log vs. cdf_log_function)

I think one is used as a reference implementation to test the other, but I'm not 100% sure. It prettymuch has to be that, though, if it's being used and defined only in the tests.

The two functions had 102 near-identical lines of doc block each: the same A&S 7.1.26 / Cody (1969) / DLMF 7.12.1 references, the same R pnorm and SciPy log_ndtr cross-references, the same #1411 provenance and quote, and the same two cutoff tables. Keeping both in sync by hand is a losing proposition -- the copy in std_normal_lcdf.hpp had already drifted, referring to `scaled_diff` in a function whose variable is `scaled_y`. Documented once on normal_lcdf. std_normal_lcdf now carries a 20-line pointer plus the only two things that genuinely differ: the scaled variable is `scaled_y = y * INV_SQRT_TWO` rather than `scaled_diff` (they coincide at mu=0, sigma=1, so every cutoff transfers unchanged), and its enforcing test is mix/prob/std_normal_cdf_log_test.cpp. The short two-line markers at each changed branch stay in both files -- a reader in the middle of the function should not have to open another file to see why the expression is shaped the way it is. Comments only; no behaviour change. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

PR #3363 rewrites normal_lccdf as normal_lcdf(-y, -mu, sigma), which routes lccdf into argument regions of lcdf that lccdf never touched before and that lcdf's own mix test barely probed -- it stopped at y = +10. Three defects live in those regions, all pre-existing on develop. Verified with git diff that normal_lcdf.hpp was untouched by #3363: the PR reaches these, it does not introduce them. Notation: s = scaled_diff, z = (y - mu) / sigma = s * sqrt(2). P1 exp(x2) sat in a denominator in the s > 2.9 derivative branch. Autodiff of 1/w needs w^(k+1) representable, so the k-th derivative dies at x2 > 709.78/(k+1). Measured at mu=0 sigma=1: d2 silently 0 for y in [26.6, 37.5] then NaN from 37.6; d3 silently 0 at y = 20, NaN from 22. Abramowitz & Stegun (1964) 7.1.26 publishes exp(-x^2) as a NUMERATOR factor; Stan had algebraically inverted it. Restoring A&S's arrangement is an exact identity -- max 6 ulp over a 200,001-point scan of s in [2.9, 26.64] -- and lets the exponential underflow instead of overflowing. R's pnorm does the same, forming the Gaussian factor only in the numerator (do_del) and guarding underflow, never overflow. P2 erfc(|s|)^2 underflowed inside log(erfc(-s)), giving a NaN third-order mixed derivative. Reachable ONLY in fvar<fvar<var>>, i.e. expect_ad's grad_hessian; fvar<var> and fvar<fvar<double>> are both clean. Window bisects to s in (-20.0000000000, -19.2103473480). Moved the crossover to R pnorm's own cutoff: it switches to the Cody tail form at y > M_SQRT_32, i.e. |x| > sqrt(32), which in scaled_diff is exactly sqrt(32)/sqrt(2) = 4. Our Cody coefficient set differs from R's p[]/q[] (ours is in s, theirs in x), so its range was measured separately: <= 9.593e-17 relative down to s = -4, degrading past -3.52. P3 The s < -29 residual correction added 0.0015065154280332 * x2, growing quadratically, but d/ds log Phi grows linearly, to 2|s| (DLMF 7.12.1). A quadratic fit cannot track a linear asymptote. Relative gradient error was 1.0e-4 at y = -54.46, 8.2e-3 at -100, and 68% at s = -1000. This corrupted plain var gradients, i.e. HMC transitions. Replaced with the truncated asymptotic Mills ratio, which needs no exponential at all: 1.3e-11 at s = -29 improving to 7.8e-18 at s = -1000. All three land at four sites together -- prim normal_lcdf, prim std_normal_lcdf, and the two OpenCL kernels, which are line-by-line transliterations whose results are compared against CPU in opencl/rev. Also in this change: - OpenCL normal_lccdf and std_normal_lccdf now reflect through lcdf like prim does. They had been left on the old erf algorithm, whose select(scaled_diff > 8.25 * INV_SQRT_TWO, 0.0, ...) makes log(0) = -inf. That is the reported CI failure: OpenCL -inf against CPU -136.16378053699881 at z = (2, 16, -0.01). - A const char* template parameter carries the caller's name into the error checks, so normal_lccdf(1, 0, -1) reports normal_lccdf rather than normal_lcdf. This breaks normal_lcdf<double, double, double>(...), an idiom that should not compile anyway; the two in-repo users were the deprecated *_cdf_log wrappers. - normal_lcdf's doc block records what each cutoff rests on: A&S 7.1.26, Cody (1969) and DLMF 7.12.1 with links, the R pnorm and SciPy log_ndtr cross-references, and the #1411 origin of the interior Taylor cutoffs, which are original and undocumented -- the PR describes them as "original Taylor expansions ... to bridge the gap where the numerical approximations are unstable". std_normal_lcdf points at it rather than duplicating; the duplicate had already drifted. - Doxygen completed on both: \ingroup prob_dists, @return, and removal of @tparam T_loc, @tparam T_scale, @PARAM mu and @PARAM sigma from std_normal_lcdf, which does not take them. Tests, all in prim and mix: - prim: 16 rows of log Phi against references carried to 40 significant digits in comments beside the correctly-rounded doubles, spanning all three value branches out to s = -212. Mutation-checked: a 0.15% perturbation of temp_p's leading coefficient fires all six Cody rows. One row is asserted loosely and documented -- at s = 11.31 this platform's libm erfc is itself 8.0e-06 out, which passes straight through. - mix: expect_ad at every branch cutoff, bracketed rather than landing on the seam, since finite_diff_grad_hessian_auto's stencil would straddle it. - mix: per-branch derivative accuracy at each branch's own measured error, because expect_ad's blanket 1e-4 leaves the (0.8, 1.5] branch only 1.6x. - mix: order-2 and order-3 tail derivatives against 60-digit references, plus a finiteness sweep over y in [-120, 120]. These are asserted directly rather than through expect_ad because expect_ad cannot see either failure: the silent zeros are ~1e-158, invisible to a finite-difference comparison whose relative tolerance floors its denominator at 1, and test_ad.hpp stops at fvar<fvar<var>>, so fvar<fvar<fvar<double>>> is never instantiated. OpenCL cannot be executed on the development machine -- the only local device has no fp64 -- so those two kernels were verified by transliterating them back to C++ and diffing against prim over 16,000 points in y in [-400, 400]: value bit-identical, gradient worst 3.6e-16. Jenkins remains the real check. Known and deliberately untouched: the s > 2.9 branch is ~6e-5 relatively inaccurate because it uses P(t)/t where A&S 7.1.26 requires P(t)/2. This change preserves that bit-for-bit and makes no accuracy claim. Fixing it means replacing all five piecewise Taylor branches with INV_SQRT_PI * exp(-x2) / (1 - 0.5 * erfc(scaled_diff)), reusing the erfc the value branch already computes, which would shift many gradient expectations. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

PhilClemson and others added 4 commits October 15, 2019 15:25

Added two valid_values with y far from mu/sigma to demonstrate the is…

d44c0e6

…sue (causes test failures).

fixes stan-dev#1284

54945c9

Merge commit 'fffe604a7ead4525be2551eb81578c5f351e5c87' into HEAD

fd701d4

[Jenkins] auto-formatting by clang-format version 6.0.0 (tags/google/…

3673fe7

…stable/2017-11-14)

bbbales2 requested changes Oct 31, 2019

View reviewed changes

SteveBronder reviewed Oct 31, 2019

View reviewed changes

PhilClemson commented Nov 1, 2019

View reviewed changes

PhilClemson added 2 commits November 14, 2019 11:47

Added the proposed changes and autodiff tests, with some additional i…

055c640

…mprovements to the numerical approximations.

Merge branch 'bugfix/issue-#1284-Numerical-precision-of-normal_lcdf' of

9345e21

https://github.com/PhilClemson/math into bugfix/issue-#1284-Numerical-precision-of-normal_lcdf

Autodiff tests moved into new mix test file

c329aca

bbbales2 requested changes Nov 16, 2019

View reviewed changes

Adjusted relative error check to make it more strict, as well as addi…

68ad6f6

…ng the other suggested changes.

Added brief comment explaining reasoning behind the extensive numeric…

59171bf

…al approximations used.

mcol reviewed Nov 18, 2019

View reviewed changes

bbbales2 requested changes Nov 19, 2019

View reviewed changes

Reverted test framework file to original (without relative error chec…

b148f2b

…ks). Also corrected a couple of typos.

bbbales2 approved these changes Nov 20, 2019

View reviewed changes

bbbales2 merged commit 858dd45 into stan-dev:develop Nov 20, 2019

PhilClemson mentioned this pull request Jul 27, 2020

Accuracy of normal_lccdf much worse than normal_lcdf #1985

Open

rok-cesnovar mentioned this pull request Apr 13, 2021

Phi_approx is not as robust in the tails as normal_lcdf stan-dev/docs#352

Open

spinkney mentioned this pull request Apr 13, 2021

Make normal_lcdf(x | 0, 1) equivalent to Phi() #2470

Open

+                  if (scaled_diff > 0.0) {
+                    // CDF(x) = 1/2 + 1/2erf(x) = 1 - 1/2erfc(x)
+                    cdf_log += log1p(-0.5 * erfc(scaled_diff));
+                    if (isnan(cdf_log)) {

+                    // based on analytic form given by:
+                    // dln(CDF)/dx = exp(-x^2)/(sqrt(pi)*(1/2+erf(x)/2)
+                    T_partials_return dncdf_log_dbl = 0.0;
+                    T_partials_return t = 0.0, t2 = 0.0, t4 = 0.0;

+                      // use fact that erf(x)=-erf(-x)
+                      // technically only true for -inf<x<0 but seems to be accurate
+                      // for -inf<x<0.1
+                      t = 1.0 / (1.0 - 0.3275911 * scaled_diff);

+                        if (scaled_diff < -4.0) {
+                          // need to add correction term (from cubic fit of residuals)
+                          dncdf_log_dbl += 0.00024073 * x2 * scaled_diff + 0.020656 * x2
+                                           + 0.10349 * scaled_diff + 0.97431;

+                        t = scaled_diff - 1.85;
+                        t2 = square(t);
+                        t4 = pow(t, 4);
+                        dncdf_log_dbl = 0.01849212058 - 0.06876280470 * t + 0.1099906382 * t2

+                    // calculate using piecewise funciton
+                    // (due to instability / inaccuracy in the various approximations)
+                    if (scaled_diff < 0.1) {

+                      // need to use direct numerical approximation of cdf_log instead
+                      // the following based on W. J. Cody, Math. Comp. 23(107):631-638 (1969)
+                      // CDF(x) = 1/2erfc(-x)
+                      typename stan::return_type<T_y, T_loc, T_scale>::type temp_p = 0.0;

-                        << i << ": " << parameters << std::endl
-                        << "  finite diffs: " << finite_dif << std::endl
-                        << "  grads:        " << gradients;
+                    // use relative error check in first case as logs of functions do not

-                    EXPECT_NEAR(finite_diffs[0], gradients[0], 1e-2);
-                  else
+                  if (!is_nan(finite_diffs[0])) {
+                    if (abs((finite_diffs[0] - gradients[0]) / finite_diffs[0]) > 0.1) {

Bugfix/issue #1284 numerical precision of normal lcdf by PhilClemson · Pull Request #1411 · stan-dev/math · GitHub

Navigation Menu

Uh oh!

Uh oh!

Bugfix/issue #1284 numerical precision of normal lcdf - #1411

Conversation

PhilClemson commented Oct 20, 2019

Summary

Tests

Side Effects

Checklist

Uh oh!

bob-carpenter commented Oct 20, 2019

Uh oh!

nhuurre commented Oct 21, 2019

Uh oh!

serban-nicusor-toptal commented Oct 21, 2019

Uh oh!

PhilClemson commented Oct 21, 2019

Uh oh!

bob-carpenter commented Oct 21, 2019 via email

Uh oh!

bbbales2 left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

SteveBronder left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

PhilClemson left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

PhilClemson commented Nov 1, 2019

Uh oh!

bbbales2 commented Nov 1, 2019

Uh oh!

bbbales2 commented Nov 1, 2019

PhilClemson commented Nov 14, 2019 •
edited
Loading

PhilClemson commented Nov 14, 2019 •
edited
Loading