使用 C 标准数学库精确计算标准正态分布的 CDF

Question

标准 C 数学库不提供计算标准正态分布的 CDF 的函数，normcdf()。然而，它确实提供了密切相关的函数：误差函数erf() 和互补误差函数erfc()。计算 CDF 最快的方法通常是通过误差函数，使用预定义的常数 M_SQRT1_2 来表示 √½：

double normcdf (double a) 
{
    return 0.5 + 0.5 * erf (M_SQRT1_2 * a);
}

显然，这会在负半平面中受到大量减法抵消，不适合大多数应用。由于使用erfc() 可以轻松避免取消问题，但其性能通常比erf() 稍低，因此最常用的推荐计算是：

double normcdf (double a) 
{
    return 0.5 * erfc (-M_SQRT1_2 * a);
}

一些测试表明，在负半平面产生的最大 ulp 误差仍然相当大。使用精确到 0.51 ulps 的erfc() 的双精度实现，可以观察到 normcdf() 中高达 1705.44 ulps。这里的问题是erfc() 的输入中的计算错误被erfc() 固有的指数缩放放大（参见answer 求幂导致的误差放大的解释）。

以下论文展示了如何在将浮点操作数与任意精度常量（例如 √½）相乘时获得（几乎）正确舍入的乘积：

Nicolas Brisebarre 和 Jean-Michel Muller，“任意精度常数的正确舍入乘法”，IEEE Transactions on Computers，卷。 57，第 2 期，2008 年 2 月，第 165-174 页

论文所提倡的方法依赖于融合乘加运算，该运算可用于所有常见处理器架构的最新实现，并通过标准数学函数fma() 在 C 中公开。这导致以下版本：

double normcdf (double a) 
{
    double SQRT_HALF_HI =  0x1.6a09e667f3bcd0p-01; //  7.0710678118654757e-01
    double SQRT_HALF_LO = -0x1.bdd3413b264560p-55; // -4.8336466567264567e-17

    return 0.5 * erfc (fma (-SQRT_HALF_HI, a, -SQRT_HALF_LO * a));
}

测试表明，与以前的版本相比，这将最大错误减少了大约一半。使用与以前相同的高精度erfc() 实现，观察到的最大误差为 842.71 ulps。这与提供误差最多为几个 ulps 的基本数学函数的通常目标相去甚远。

是否有一种有效的方法可以精确计算 normcdf()，并且只使用标准 C 数学库中可用的函数？

Answer 1

解决问题中概述的方法的准确性限制的一种方法是有限使用双双计算。这涉及将-sqrt (0.5) * a 计算为一对double 变量h 和l 以头/尾方式。乘积的高阶部分h 被传递给erfc()，而低阶部分l 然后用于根据@ 处互补误差函数的局部斜率对erfc() 结果进行插值987654330@.

erfc(x) 的导数是 -2 * exp (-x * x) / √π。然而，人们希望避免 exp(-x * x) 的相当昂贵的计算。 known 对于 x > 0，erfc(x) ~= 2 * exp (-x * x) / (√π * (x + sqrt (x* x + 4/π))。因此，渐近， erfc'(x) ~= -2 * x * erfc(x)，因此对于 |l| ≪|h|, erfc (h+l) ~= erfc (h) - 2 * h * l * erfc(h)。后一项的否定很容易被拉入l的计算中。一个到达以下双精度实现（使用 IEEE-754 binary64）：

double my_normcdf (double a)
{
    double h, l, r;
    const double SQRT_HALF_HI =  0x1.6a09e667f3bcd0p-01; //  7.0710678118654757e-01
    const double SQRT_HALF_LO = -0x1.bdd3413b264560p-55; // -4.8336466567264567e-17

    /* clamp input as normcdf(x) is either 0 or 1 asymptotically */
    if (fabs (a) > 38.625) a = (a < 0.0) ? -38.625 : 38.625;

    h = fma (-SQRT_HALF_HI, a, -SQRT_HALF_LO * a);
    l = fma (SQRT_HALF_LO, a, fma (SQRT_HALF_HI, a, h));
    r = erfc (h);
    if (h > 0.0) r = fma (2.0 * h * l, r, r);
    return 0.5 * r;
}

使用与之前相同的erfc() 实现，观察到的最大错误为 1.96 ulps。对应的单精度实现（使用 IEEE-754 binary32）是：

float my_normcdff (float a)
{
    float h, l, r;
    const float SQRT_HALF_HI = 0x1.6a09e6p-01f; // 7.07106769e-1
    const float SQRT_HALF_LO = 0x1.9fcef4p-27f; // 1.21016175e-8

    /* clamp input as normcdf(x) is either 0 or 1 asymptotically */
    if (fabsf (a) > 14.171875f) a = (a < 0.0f) ? -14.171875f : 14.171875f;

    h = fmaf (-SQRT_HALF_HI, a, -SQRT_HALF_LO * a);
    l = fmaf (SQRT_HALF_LO, a, fmaf (SQRT_HALF_HI, a, h));
    r = erfcf (h);
    if (h > 0.0f) r = fmaf (2.0f * h * l, r, r);
    return 0.5f * r;
}