Repository navigation
Make normal_lcdf(x | 0, 1) equivalent to Phi() #2470
Description
Since #1411 normal_lcdf is more robust in the tails than Phi or even Phi_approx. Update Phi() to be as robust as normal_lcdf and add log_Phi() as discussed in the linked PR above (should this be a separate issue?)
In fact, Phi_approx() is less robust than Phi(). Adding @PhilClemson @nhuurre @bbbales2
generated quantities {
vector[N] b;
vector[N] c;
vector[N] d;
for (i in 1:N) {
b[i] = normal_lcdf(i | 0, 1);
c[i] = log(Phi_approx(i));
d[i] = log(Phi(i));
}
}> data.table(mean_b = test_b, mean_c = test_c, mean_d = test_d)
mean_b mean_c mean_d
1: -1.72754e-01 -1.72771e-01 -1.72754e-01
2: -2.30129e-02 -2.30241e-02 -2.30129e-02
3: -1.35081e-03 -1.23271e-03 -1.35081e-03
4: -3.16717e-05 -1.83432e-05 -3.16717e-05
5: -2.86652e-07 -5.01624e-08 -2.86652e-07
6: -9.86588e-10 -1.65181e-11 -9.86588e-10
7: -1.27981e-12 -4.44089e-16 -1.27987e-12
8: -6.22096e-16 0.00000e+00 -6.66134e-16
9: -1.12859e-19 0.00000e+00 0.00000e+00
10: -7.61985e-24 0.00000e+00 0.00000e+00
11: -1.91066e-28 0.00000e+00 0.00000e+00
12: -1.77648e-33 0.00000e+00 0.00000e+00
13: -6.11716e-39 0.00000e+00 0.00000e+00
14: -7.79354e-45 0.00000e+00 0.00000e+00
15: -3.67097e-51 0.00000e+00 0.00000e+00
16: -6.38875e-58 0.00000e+00 0.00000e+00
17: -4.10600e-65 0.00000e+00 0.00000e+00
18: -9.74095e-73 0.00000e+00 0.00000e+00
19: -8.52722e-81 0.00000e+00 0.00000e+00
20: -2.75362e-89 0.00000e+00 0.00000e+00
21: -3.27928e-98 0.00000e+00 0.00000e+00
22: -1.43989e-107 0.00000e+00 0.00000e+00
23: -2.33064e-117 0.00000e+00 0.00000e+00
24: -1.39039e-127 0.00000e+00 0.00000e+00
25: -3.05670e-138 0.00000e+00 0.00000e+00
26: -2.47606e-149 0.00000e+00 0.00000e+00
27: -7.38948e-161 0.00000e+00 0.00000e+00
28: -8.12387e-173 0.00000e+00 0.00000e+00
29: -3.28979e-185 0.00000e+00 0.00000e+00
30: -4.90671e-198 0.00000e+00 0.00000e+00
31: -2.69525e-211 0.00000e+00 0.00000e+00
32: -5.45208e-225 0.00000e+00 0.00000e+00
33: -4.06119e-239 0.00000e+00 0.00000e+00
34: -1.11390e-253 0.00000e+00 0.00000e+00
35: -1.12491e-268 0.00000e+00 0.00000e+00
36: -4.18262e-284 0.00000e+00 0.00000e+00
37: -5.72557e-300 0.00000e+00 0.00000e+00
38: -2.88543e-316 0.00000e+00 0.00000e+00
39: 0.00000e+00 0.00000e+00 0.00000e+00
40: 0.00000e+00 0.00000e+00 0.00000e+00
Activity
andrjohns commentedon Apr 16, 2021
It looks while
log(Phi(x))does not perform as well asstd_normal_lcdf(x), there is only a marginal difference betweenPhi(x)andexp(std_normal_lcdf(x)). Also,log(Phi_approx(x))performs better (but not equivalently) whenlog_inv_logitis used.Testing
Model functions
log(Phi(x))Results:Phi(x)Results:bbbales2 commentedon Apr 16, 2021
Also normal_lccdf should be updated too -- there might be a separate issue for that. Tnx for looking at this.
PhilClemson commentedon Apr 19, 2021
I'm not sure if there's a way to get
log(Phi(x))to replicate the precision ofnormal_lcdf(x|0,1). The precision is achieved through analytic approximations of the log function. Usinglog()severely reduces the precision and leads to faster overflow to 0. Alog_Phi()function definitely makes sense though, as does updatingnormal_lccdf().Also, I just happened to check the function reference doc and it still quotes the old values for underflow / overflow (-37.5 and 8.25 respectively) so this should probably be updated too.
spinkney commentedon Apr 19, 2021
Here is the docs issue :) stan-dev/docs#352
avehtari commentedon Oct 23, 2026
I checked the three requests in this issue against mpmath, on
developwith PR #3406.Phias robust asnormal_lcdf. With Expanded numerical stability and performance forstd_normal_lcdf, reused across distributions & functions #3406,Phihas no cutoffs. Before, it returned 0 below -37.5 and 1 above 8.25. Its relative error for z in [-37.5, 8] went from 4e-11 (at z = -5) to 4e-16. On a grid from -40 to 40 with step 1/64, there is now no point where one ofPhi(z)andexp(std_normal_lcdf(z))is 0 and the other is not. Before there were 63, all between -38.5 and -37.5. In the deep lower tailPhiis the more accurate of the two, becauseexpturns the error of a log value near -600 into a relative error of about 2e-13.The upper tail cannot be fixed in
Phi. Above z = 8.30,Phi(z)rounds to 1 in double precision, solog(Phi(z))is 0 there with any implementation ofPhi. Onlystd_normal_lcdfkeeps that information, or the complementstd_normal_lccdfandPhi(-z).log_Phi(). This is still missing.std_normal_lcdf(x)gives log Phi for a scalar, but for a vector it returns the sum, so Stan has no elementwise log Phi. The copula use case in Phi_approx is not as robust in the tails as normal_lcdf docs#352 would benefit from one. The inverse already exists as the elementwisestd_normal_log_qf. Followinglog_inv_logitandlog1m_inv_logit, the new functions would belog_Phi(T x)andlog1m_Phi(T x). After Expanded numerical stability and performance forstd_normal_lcdf, reused across distributions & functions #3406 the Stan Math part is small, becauseinternal::std_normal_lcdf_value_gradalready computes the elementwise values and slopes, and the OpenCL device functions exist. The functions also need a stanc3 signature and a docs entry. I can make the Stan Math PR if you want these functions.normal_lccdf. Aligned normal_lccdf and std_normal_lccdf with the lcdf functions #3363 fixed this (merged 2026-08-30) and closed Accuracy of normal_lccdf much worse than normal_lcdf #1985.The functions reference still gave the old limits and said that
Phi_approxis more robust in the tails. stan-dev/docs#990 corrects both.