c - 使用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() double 实现，可以观察到误差 在 normcdf() 中高达 1705.44 ulps。这里的问题是 erfc() 输入中的计算误差被 erfc() 固有的指数缩放放大了(参见这个 answer 求幂引起的误差放大的解释)。

以下论文展示了如何在将浮点操作数与任意精度常量(例如 √½)相乘时实现(几乎)正确舍入的乘积:

Nicolas Brisebarre 和 Jean-Michel Muller，“任意精度常数的正确舍入乘法”，IEEE 计算机交易，卷。 57, No. 2, February 2008, pp. 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。这与提供误差至多几个 ulp 的基本数学函数的通常目标相去甚远。

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

Answer 1

解决问题中概述的方法的准确性限制的一种方法是限制使用 double-double 计算。这涉及计算 -sqrt (0.5) * a 作为一对 double 变量 h 和 l头/尾时尚。乘积的高阶部分 h 被传递给 erfc()，而低阶部分 l 则用于插值erfc() 结果，基于 h 处互补误差函数的局部斜率。

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 的计算中。一个到达以下 double 实现(使用 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;
}