 |
Stan Math Library
5.3.0
Automatic Differentiation
|
|
|
Stan Math Library
5.3.0
Automatic Differentiation
|
|
|
inline |
Calculates the log of the cdf of the normal distribution.
Tail branching follows three published results, and matches what other libraries do:
erf(x) = 1 - P(t) exp(-x^2), t = 1/(1 + p x), with exp(-x^2) as a numerator factor: https://archive.org/details/handbookofmathem1964abra/page/298/mode/2upR's pnorm only ever forms the Gaussian factor in the numerator (the do_del macro), and switches to the Cody tail form at y > M_SQRT_32, i.e. |x| > sqrt(32) ~= 5.657. Since scaled_diff = x / sqrt(2), that same crossover is |scaled_diff| > sqrt(32)/sqrt(2) = 4 exactly, which is the cutoff used here for the cdf value. Measured by evaluating the temp_p/temp_q expression below against log(erfc(-scaled_diff)/2) at 60 significant digits, our Cody set holds to 8.6e-17 relative down to scaled_diff = -4 and degrades past about -3.52, so 4 sits just inside its range. This is the same as R's impl. https://github.com/wch/r-source/blob/trunk/src/nmath/pnorm.c SciPy's log_ndtr uses the identical log1p(-erfc(t)/2) upper branch with the same t = x/sqrt(2), and needs no rational approximation or Taylor patches at all below x = -1, because erfcx never forms exp(+t^2): https://github.com/scipy/xsf/blob/main/include/xsf/stats.h
The interior cutoffs for the gradient are not from the literature. They were derived in stan-dev/math#1411 (Phil Clemson, Univ. of Liverpool, Nov 2019), fixing stan-dev/math#1284, which describes them as "original Taylor expansions that have been derived to bridge the gap where the numerical approximations are unstable". The author's account of how they were placed, from the review thread, was: "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." https://github.com/stan-dev/math/pull/1411 https://github.com/stan-dev/math/issues/1284
Each row of the chart below is one Taylor expansion around scaled_diff. interval is the range of scaled_diff it covers. centre is the point the series is expanded about.
| interval | centre | half-width |
|---|---|---|
| (2.5, 2.9] | 2.7 | 0.200 |
| (2.1, 2.5] | 2.3 | 0.200 |
| (1.5, 2.1] | 1.85 | 0.300 |
| (0.8, 1.5] | 1.15 | 0.350 |
| (0.1, 0.8] | 0.45 | 0.350 |
The cut points are not arbitrary. Each interval of scaled_diff is centred on its own Taylor point, which minimises the largest |t| the series has to cover.
The second table below shows why each interval ends where it does: the series is accurate across its own range and falls apart just outside it, so the intervals cannot be widened to use fewer of them.
Worst in-range comes from evaluating each Taylor polynomial as written against (2/sqrt(pi)) * exp(-scaled_diff^2) / erfc(-scaled_diff), the exact derivative, at 60 significant digits.
Columns below: worst in-range: The largest relative error of that branch's series over its own interval; 0.2 below lo: The relative error the same series would give at lo - 0.2; 0.2 above hi: The relative error the same series would give at hi + 0.2, i.e. what widening the interval either way would cost.
| interval | worst in-range | 0.2 below lo | 0.2 above hi |
|---|---|---|---|
| (2.5, 2.9] | 3.27e-05 | 2.32e-05 | 1.61e-02 |
| (2.1, 2.5] | 2.80e-05 | 3.25e-04 | 9.15e-03 |
| (1.5, 2.1] | 2.30e-05 | 8.11e-05 | 3.77e-03 |
| (0.8, 1.5] | 6.09e-05 | 9.29e-05 | 2.93e-03 |
| (0.1, 0.8] | 7.61e-06 | 3.62e-05 | 2.82e-04 |
The worst in-range error is uniformly 1e-5 to 6e-5, just inside the 1e-4 relative gradient tolerance expect_ad applies by default (gradient_grad_ in test/unit/math/ad_tolerances.hpp; the 1e-8 there is gradient_val_, which bounds the value rather than the gradient), while widening an interval by 0.2 costs one to three orders of magnitude. That uniformity, not any published result, is the placement criterion.
The worst in-range column is enforced: the branch_accuracy test in mix/prob/normal_cdf_log_test.cpp asserts every branch against 60-digit references at that error plus a small margin.
The negative-tail residual corrections at scaled_diff of -2.1, -3.9, -7 and -17 are cubic fits of residuals from the same PR and likewise have no external source. The -29 cutoff below is justified by DLMF 7.12.1
| func | name reported by the error checks. Reflected distributions such as normal_lccdf delegate here and pass their own name so that exceptions name the function the user actually called. |
| T_y | A vector or scalar type for the random variable. |
| T_loc | A vector or scalar type for the location parameter. |
| T_scale | A vector or scalar type for the scale parameter. |
| y | (Sequence of) scalar(s). |
| mu | (Sequence of) scalar(s). |
| sigma | (Sequence of) scalar(s). |
Definition at line 143 of file normal_lcdf.hpp.