高效准确地计算标准正态分布的米尔斯比

Question

Mills ratio M(x) 由 John Mills 引入，用于表达分布的累积分布函数与其概率密度函数之间的关系：

J。 P.米尔斯，"Table of the ratio: Area to bounding ordinate, for any portion of normal curve"。 Biometrika，卷。 18，第 3/4 期（1926 年 11 月），第 395-400 页。 (online)

米尔斯比的definition为(1 - D(x)) / P(x)，其中D表示分布函数，P(x)为概率密度函数。在标准正态分布的特定情况下，我们有 M(x) = (1 - Φ(x)) / ϕ(x) = Φ(-x) / ϕ(x)，或者当通过互补误差表示时函数，M(x) = e^x² √(π/2) erfc (x/√2) = √(π/2) erfcx (x/√2).

之前的问题涉及 R and Matlab 等数学环境中米尔斯比的计算，但这些环境的复杂计算设施在 C 中没有等效项。如何计算标准的米尔斯比仅使用 C 标准数学库准确高效地进行正态分布？

Answer 1

在之前的回答中，我展示了如何使用 C 标准数学库高效准确地计算 , normpdf(), the , normcdf(), and the erfcx()。基于这三种实现，可以通过以下两种方式之一以直接的方式轻松编写米尔斯比的计算代码：

double my_mills_ratio_1 (double a)
{
    return my_normcdf (-a) / my_normpdf (a);
}

double my_mills_ratio_2 (double a)
{
    const double SQRT_HALF_HI = 0x1.6a09e667f3bccp-01; // 1/sqrt(2), msbs;
    const double SQRT_HALF_LO = 0x1.21165f626cdd5p-54; // 1/sqrt(2), lsbs;
    const double SQRT_PIO2_HI = 0x1.40d931ff62705p+00; // sqrt(pi/2), msbs;
    const double SQRT_PIO2_LO = 0x1.2caf9483f5ce4p-53; // sqrt(pi/2), lsbs;
    double r;

    a = fma (SQRT_HALF_HI, a, SQRT_HALF_LO * a);
    r = my_erfcx (a);
    return fma (SQRT_PIO2_HI, r, SQRT_PIO2_LO * r);
}

但是，这两种方法在数值上都存在缺陷。对于 mills_ratio_1()，PDF 项和 CDF 项都在正半平面中随着自变量的增加而迅速消失。在 IEEE-754 中，双精度都在 a = 38 附近变为零，由于零被零除而导致 NaN 结果。关于 my_mills_ratio_2()，负半平面的指数增长导致误差放大，因此导致大的 ulp 误差。解决这个问题的一种方法是简单地组合两个近似值中每一个的表现良好的部分：

double my_mills_ratio_3 (double a)
{
    return (a < 0) ? my_mills_ratio_1 (a) : my_mills_ratio_2 (a);
}

这相当有效。使用英特尔编译器版本 13.1.3.198 构建我在之前的答案中提供的代码，使用 40 亿个测试向量，在正半平面中观察到的最大误差为 2.79346 ulps，而在正半平面中观察到的最大误差为 6.81248 ulps负半平面。对于接近溢出的大结果，负半平面中的误差稍大一些，因为此时 PDF 的值非常小，以至于它们被表示为精度降低的次正规双精度数。

另一种解决方案是解决影响负半平面中 my_mills_ratio_2() 的误差放大问题。可以通过将 erfcx() 的参数计算为优于双精度，并使用该参数的低位位对 erfcx() 结果进行线性插值来实现。

为此，还需要erfcx(x)的斜率，即erfcx'(x) = 2x erfcx(x) - 2/√π。通过 C 标准数学函数 fma() 的 FMA（融合乘加）运算的可用性提供了这种准双双计算的有效实现。可以通过局部重新缩放来避免中间计算期间斜率大小溢出的风险。

最终实现在整个输入域中的误差小于 4 ulps：

/* Compute Mills ratio of the standard normal distribution:
 *
 * M(x) = normcdf(-x)/normpdf(x) = sqrt(pi/2) * erfcx(x/sqrt(2)) 
 *
 * maximum ulp error in positive half-plane: 2.79346
 * maximum ulp error in negative half-plane: 3.90753
 */ 
double my_mills_ratio (double a)
{
    double s, t, r, h, l;

    const double SQRT_HALF_HI = 0x1.6a09e667f3bccp-01; // 1/sqrt(2), msbs
    const double SQRT_HALF_LO = 0x1.21165f626cdd5p-54; // 1/sqrt(2), lsbs
    const double SQRT_PIO2_HI = 0x1.40d931ff62705p+00; // sqrt(pi/2), msbs
    const double SQRT_PIO2_LO = 0x1.2caf9483f5ce4p-53; // sqrt(pi/2), lsbs
    const double TWO_RSQRT_PI = 0x1.20dd750429b6dp+00; // 2/sqrt(pi)
    const double MAX_IEEE_DBL = 0x1.fffffffffffffp+1023;
    const double SCALE_DOWN = 0.03125; // prevent ovrfl in intermed. computation
    const double SCALE_UP = 1.0 / SCALE_DOWN;

    // Compute argument a/sqrt(2) as a head-tail pair of doubles h:l
    h = fma (SQRT_HALF_HI, a, SQRT_HALF_LO * a);
    l = fma (-SQRT_HALF_LO, a, fma (-SQRT_HALF_HI, a, h));
    // Compute scaled complementary error function for argument "head"
    t = my_erfcx (h);
    // Enhance accuracy if in negative half-plane, if result has not overflowed
    if ((a < -1.0)  && (t <= MAX_IEEE_DBL)) {
        // Compute slope: erfcx'(x) = 2x * erfcx(x) - 2/sqrt(pi)
        s = fma (h, t * SCALE_DOWN, -TWO_RSQRT_PI * SCALE_DOWN); // slope
        // Linearly interpolate result based on derivative and argument "tail"
        t = fma (s, -2.0 * SCALE_UP * l, t);
    }
    // Scale by sqrt(pi/2) for final result
    r = fma (SQRT_PIO2_HI, t, SQRT_PIO2_LO * t);
    return r;
}

单精度实现看起来几乎相同，除了涉及的常量：

/* Compute Mills ratio of the standard normal distribution:
 *
 * M(x) = normcdf(-x)/normpdf(x) = sqrt(pi/2) * erfcx(x/sqrt(2)) 
 *
 * maximum ulp error in positive half-plane: 2.41987
 * maximum ulp error in negative half-plane: 3.39521
 */ 
float my_mills_ratio_f (float a)
{
    float h, l, r, s, t;

    const float SQRT_HALF_HI = 0x1.6a09e6p-01f; // sqrt(1/2), msbs
    const float SQRT_HALF_LO = 0x1.9fcef4p-27f; // sqrt(1/2), lsbs
    const float SQRT_PIO2_HI = 0x1.40d930p+00f; // sqrt(pi/2), msbs
    const float SQRT_PIO2_LO = 0x1.ff6270p-24f; // sqrt(pi/2), lsbs
    const float TWO_RSQRT_PI = 0x1.20dd76p+00f; // 2/sqrt(pi)
    const float MAX_IEEE_FLT = 0x1.fffffep+127f;
    const float SCALE_DOWN = 0.0625f; // prevent ovrfl in intermed. computation
    const float SCALE_UP = 1.0f / SCALE_DOWN;

    // Compute argument a/sqrt(2) as a head-tail pair of floats h:l
    h = fmaf (SQRT_HALF_HI, a, SQRT_HALF_LO * a);
    l = fmaf (-SQRT_HALF_LO, a, fmaf (-SQRT_HALF_HI, a, h));
    // Compute scaled complementary error function for argument "head"
    t = my_erfcxf (h);
    // Enhance accuracy if in negative half-plane, if result has not overflowed
    if ((a < -1.0f)  && (t <= MAX_IEEE_FLT)) {
        // Compute slope: erfcx'(x) = 2x * erfcx(x) - 2/sqrt(pi)
        s = fmaf (h, t * SCALE_DOWN, -TWO_RSQRT_PI * SCALE_DOWN);
        // Linearly interpolate result based on derivative and argument "tail"
        t = fmaf (s, -2.0f * SCALE_UP * l, t);
    }
    // Scale by sqrt(pi/2) for final result
    r = fmaf (SQRT_PIO2_HI, t, SQRT_PIO2_LO * t);
    return r;
}

高效准确地计算标准正态分布的米尔斯比

Efficient and accurate computation of the Mills ratio of the standard normal distribution

c

math

floating-point

floating-accuracy