用C标准数学库精确计算标准正态分布的CDF [英] Accurate computation of CDF of standard normal distribution using C standard math library

问题描述

标准的C数学库不提供一个函数来计算标准正态分布的CDF， normcdf（）。但是，它提供了密切相关的函数：错误函数 erf（）和补充错误函数 erfc（）。计算CDF的最快方法通常是通过错误函数，使用预定义常量M_SQRT1_2来表示√½：

  double normcdf （double a）
 {
 return 0.5 + 0.5 * erf（M_SQRT1_2 * a）; 
 
 
 
 
 
 $ b显然，这在负半平面上遭受了大量的减法消除，不适合大多数应用。由于使用 erfc（）很容易避免取消问题，但是通常性能比 erf（），最常推荐的计算是：
 
 pre $ $ $ $ $ $ $ $ 0.5 * erfc（-M_SQRT1_2 * a）; 
 
 $ / code> 

有些测试显示负半平面上的最大ulp误差是还是比较大的。使用 erfc（）精确到0.51 ulps的双精度实现，可以在 normcdf中观察到
的错误，高达1705.44 ulps ）。这里的问题是，对 erfc（）的输入的计算错误被 erfc（）（见回答
的解释错误放大倍数）。

下面的文章展示了当浮点操作数乘以任意精度常数（如在本文中，我们将讨论如何使用任意精度常数来正确舍入乘法运算。IEEE Transactions on Computers，第2卷， 57，第2号，2008年2月，pp.165-174

本文提出的方法依赖于融合的乘加操作，这个操作在最近的实现所有常见的处理器体系结构，并通过标准的数学函数 fma（）在C中公开。这导致以下版本：

pre $ 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））;

$ / code>

测试结果显示，与前一版本相比， 。像以前一样使用高精度的 erfc（）实现，观察到的最大错误是842.71个ulps。这与通常的提供基本数学函数的目标相距甚远，最多只有几个ulps。



是否有一个有效的方法可以进行精确的计算 normcdf（），并且只使用标准C数学库中可用的函数？

解决方案

这个问题中概述的方法的精度限制的一个方法是有限的使用双重计算。这包括将 -sqrt（0.5）* a 作为一对 double 变量 h / code>和 l 以头/尾的方式。产品的高阶部分 h 传递给 erfc（），而低阶部分<$然后用c $ c> l 来插入 erfc（）结果，这个结果是根据互补误差函数在<$ c
$ b $ erfc（x）的导数是-2 * exp（-x * x）/√π。但是，人们希望避免exp（-x * x）的相当昂贵的计算。这是已知对于x> 0，erfc（x）〜= 2 * exp（-x * x）/（√π*（x + sqrt（x * x + 4 /π））因此，渐近地，
erfc'（x）〜= -2 * x * erfc那么对于| l |«| h |，erfc（h + l）〜= erfc（h） - 2 * h * l * erfc（h）可以很容易地推导出后一项的否定。 code> 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 
 
 / *钳制输入为normcdf（x）如果（fabs（a）> 38.625）a =（a <0.0）-38.625：38.625; 
 
h = fma（-SQRT_HALF_HI，a ，-SQRT_HALF_LO * a）; 
l = fm a（SQRT_HALF_LO，a，fma（SQRT_HALF_HI，a，h））; 
 r = erfc（h）;如果（h> 0.0）r = fma（2.0 * h * l，r，r）则
; 
 return 0.5 * r; 
 
 
 
 $ b 
使用相同的 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;如果（fabsf（a）> 14.171875f）a =（a），则输入为normcdf（x）或者为0或1渐近* / 
 < 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）;如果（h> 0.0f）r = fmaf（2.0f * h * l，r，r）则
; 
 return 0.5f * r; 
} 
  
 
The standard C math library does not provide a function to compute the CDF of the standard normal distribution, normcdf(). It does, however, offer closely related functions: the error function, erf() and the complementary error function, erfc(). The fastest way to compute the CDF is often via the error function, using the predefined constant M_SQRT1_2 to represent √½:
double normcdf (double a) 
{
    return 0.5 + 0.5 * erf (M_SQRT1_2 * a);
}
Clearly, this suffers from massive subtractive cancellation in the negative half-plane and is not  suitable for the majority of applications. Since the cancellation issue is easily avoided by use of erfc(), which however is typically of somewhat lower performance than erf(), the most frequently recommended computation is:
double normcdf (double a) 
{
    return 0.5 * erfc (-M_SQRT1_2 * a);
}
Some testing shows that the maximum ulp error incurred in the negative half-plane is still rather large however. Using a double-precision implementation of erfc() accurate to 0.51 ulps, one can observe errors of 
up to 1705.44 ulps in normcdf(). The issue here is that the computational error in the input to erfc() is magnified by the exponential scaling that is inherent to erfc() (See this answer 
for an explanation of error magnification caused by exponentiation).


The following paper shows how one can achieve (nearly) correctly rounded, products when multiplying floating-point operands with arbitrary-precision constants such as √½:

Nicolas Brisebarre and Jean-Michel Muller, "Correctly rounded multiplication by arbitrary precision constants", IEEE Transactions on Computers, Vol. 57, No. 2, February 2008, pp. 165-174

The method advocated by the paper relies on the fused multiply-add operation, which is available on recent implementations of all common processor architectures, and is exposed in C via the standard math function fma(). This leads to the following version:
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));
}
Tests show that this cuts the maximum error roughly in half compared to the previous version. Using the same highly-accurate erfc() implementation as before, the maximum observed error is 842.71 ulps. This is still far away from the usual goal of providing basic mathematical functions with an error of at most a few ulps.

Is there an efficient method that allows for the accurate computation of normcdf(), and that uses only functions available in the standard C math library?
解决方案
One way around the accuracy limitations of the approaches outlined in the question is the limited use of double-double computation. This involves the computation of -sqrt (0.5) * a as a pair of double variables h and l in head/tail fashion. The high-order portion h of the product is passed to erfc(), while the low-order portion l is then used to interpolate the erfc() result, based on the local slope of the complementary error function at h.

The derivative of erfc(x) is -2 * exp (-x * x) / √π. However, one would like to avoid the fairly expensive computation of exp(-x * x). It is known that for x > 0, erfc(x) ~= 2 * exp (-x * x) / (√π * (x + sqrt (x* x + 4/π)). Therefore, asymptotically, 
erfc'(x) ~= -2 * x * erfc(x), and it follows that for |l| ≪|h|, erfc (h+l) ~= erfc (h) - 2 * h * l * erfc(h). The negation of the latter term can easily be pulled into the computation of l. One arrives at the following implementation for double precision (using 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;
}
The maximum observed error, using the same erfc() implementation as used before, is 1.96 ulps. The corresponding single-precision implementation (using IEEE-754 binary32) is:
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;
}


                        
这篇关于用C标准数学库精确计算标准正态分布的CDF的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！