累积正态分布函数在C/ c++中。

I was wondering if there were statistics functions built into math libraries that are part of the standard C++ libraries like cmath. If not, can you guys recommend a good stats library that would have a cumulative normal distribution function? Thanks in advance.

我想知道在数学库中是否有一些统计功能，它们是标准c++库(如cmath)的一部分。如果没有，你们能推荐一个有累积正态分布函数的好统计库吗?提前谢谢。

More specifically, I am looking to use/create a cumulative distribution function.

更具体地说，我希望使用/创建一个累积分布函数。

6 个解决方案

#1

Here's a stand-alone C++ implementation of the cumulative normal distribution in 14 lines of code.

这是一个独立的c++实现，在14行代码中累积正态分布。

http://www.johndcook.com/cpp_phi.html

#include <cmath>

double phi(double x)
{
    // constants
    double a1 =  0.254829592;
    double a2 = -0.284496736;
    double a3 =  1.421413741;
    double a4 = -1.453152027;
    double a5 =  1.061405429;
    double p  =  0.3275911;

    // Save the sign of x
    int sign = 1;
    if (x < 0)
        sign = -1;
    x = fabs(x)/sqrt(2.0);

    // A&S formula 7.1.26
    double t = 1.0/(1.0 + p*x);
    double y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*exp(-x*x);

    return 0.5*(1.0 + sign*y);
}

void testPhi()
{
    // Select a few input values
    double x[] = 
    {
        -3, 
        -1, 
        0.0, 
        0.5, 
        2.1 
    };

    // Output computed by Mathematica
    // y = Phi[x]
    double y[] = 
    { 
        0.00134989803163, 
        0.158655253931, 
        0.5, 
        0.691462461274, 
        0.982135579437 
    };

        int numTests = sizeof(x)/sizeof(double);

    double maxError = 0.0;
    for (int i = 0; i < numTests; ++i)
    {
        double error = fabs(y[i] - phi(x[i]));
        if (error > maxError)
            maxError = error;
    }

        std::cout << "Maximum error: " << maxError << "\n";
}

#2

Theres is no straight function. But since the gaussian error function and its complementary function is related to the normal cumulative distribution function (see here) we can use the implemented c-function erfc:

它不是直线函数。但由于高斯误差函数及其补充函数与正态累积分布函数有关(见此处)，我们可以使用实现的c函数erfc:

double normalCFD(double value)
{
   return 0.5 * erfc(-value * M_SQRT1_2);
}

I use it for statistical calculations and it works great. No need for using coefficients.

我用它来进行统计计算，效果很好。不需要使用系数。

#3

Boost is as good as the standard :D here you go: boost maths/statistical.

Boost和标准一样好:D在这里:提高数学/统计学。

#4

I figured out how to do it using gsl, at the suggestion of the folks who answered before me, but then found a non-library solution (hopefully this helps many people out there who are looking for it like I was):

我找到了如何使用gsl的方法，根据在我之前回答的人的建议，然后找到了一个非库的解决方案(希望这能帮助许多正在寻找它的人，像我一样):

#ifndef Pi 
#define Pi 3.141592653589793238462643 
#endif 

double cnd_manual(double x)
{
  double L, K, w ;
  /* constants */
  double const a1 = 0.31938153, a2 = -0.356563782, a3 = 1.781477937;
  double const a4 = -1.821255978, a5 = 1.330274429;

  L = fabs(x);
  K = 1.0 / (1.0 + 0.2316419 * L);
  w = 1.0 - 1.0 / sqrt(2 * Pi) * exp(-L *L / 2) * (a1 * K + a2 * K *K + a3 * pow(K,3) + a4 * pow(K,4) + a5 * pow(K,5));

  if (x < 0 ){
    w= 1.0 - w;
  }
  return w;
}

#5

The implementations of the normal CDF given here are single precision approximations that have had float replaced with double and hence are only accurate to 7 or 8 significant (decimal) figures.
For a VB implementation of Hart's double precision approximation, see figure 2 of West's Better approximations to cumulative normal functions.

这里给出的普通CDF的实现是单精度逼近，它的浮点数替换为double，因此只能精确到7或8个重要(十进制)数字。对于一个VB实现的Hart的双精度逼近，参见图2的西方更好的近似于累积的正常函数。

Edit: My translation of West's implementation into C++:

编辑:我把西方的实现翻译成c++:

double
phi(double x)
{
  static const double RT2PI = sqrt(4.0*acos(0.0));

  static const double SPLIT = 7.07106781186547;

  static const double N0 = 220.206867912376;
  static const double N1 = 221.213596169931;
  static const double N2 = 112.079291497871;
  static const double N3 = 33.912866078383;
  static const double N4 = 6.37396220353165;
  static const double N5 = 0.700383064443688;
  static const double N6 = 3.52624965998911e-02;
  static const double M0 = 440.413735824752;
  static const double M1 = 793.826512519948;
  static const double M2 = 637.333633378831;
  static const double M3 = 296.564248779674;
  static const double M4 = 86.7807322029461;
  static const double M5 = 16.064177579207;
  static const double M6 = 1.75566716318264;
  static const double M7 = 8.83883476483184e-02;

  const double z = fabs(x);
  double c = 0.0;

  if(z<=37.0)
  {
    const double e = exp(-z*z/2.0);
    if(z<SPLIT)
    {
      const double n = (((((N6*z + N5)*z + N4)*z + N3)*z + N2)*z + N1)*z + N0;
      const double d = ((((((M7*z + M6)*z + M5)*z + M4)*z + M3)*z + M2)*z + M1)*z + M0;
      c = e*n/d;
    }
    else
    {
      const double f = z + 1.0/(z + 2.0/(z + 3.0/(z + 4.0/(z + 13.0/20.0))));
      c = e/(RT2PI*f);
    }
  }
  return x<=0.0 ? c : 1-c;
}

Note that I have rearranged expressions into the more familiar forms for series and continued fraction approximations. The last magic number in West's code is the square root of 2π, which I've deferred to the compiler on the first line by exploiting the identity acos(0) = ½ π.
I've triple checked the magic numbers, but there's always the chance that I've mistyped something. If you spot a typo, please comment!

注意，我已经将表达式重新排列成更常见的级数和连续的分数近似。在西方代码中，最后一个神奇的数字是√2，我通过使用identity acos(0) =，将其递延到第一行的编译器上。我已经三次检查了这些神奇的数字，但总有可能我把一些东西搞错了。如果你发现一个错误，请评论!

The results for the test data John Cook used in his answer are

约翰·库克在回答中使用的测试数据的结果是。

 x               phi                Mathematica
-3     1.3498980316301150e-003    0.00134989803163
-1     1.5865525393145702e-001    0.158655253931
 0     5.0000000000000000e-001    0.5
0.5    6.9146246127401301e-001    0.691462461274
2.1    9.8213557943718344e-001    0.982135579437

I take some small comfort from the fact that they agree to all of the digits given for the Mathematica results.

我从他们同意所有为Mathematica结果给出的数字中得到一些小小的安慰。

#6

From NVIDIA CUDA samples:

从NVIDIA CUDA样本:

static double CND(double d)
{
    const double       A1 = 0.31938153;
    const double       A2 = -0.356563782;
    const double       A3 = 1.781477937;
    const double       A4 = -1.821255978;
    const double       A5 = 1.330274429;
    const double RSQRT2PI = 0.39894228040143267793994605993438;

    double
    K = 1.0 / (1.0 + 0.2316419 * fabs(d));

    double
    cnd = RSQRT2PI * exp(- 0.5 * d * d) *
          (K * (A1 + K * (A2 + K * (A3 + K * (A4 + K * A5)))));

    if (d > 0)
        cnd = 1.0 - cnd;

    return cnd;
}

#1