当前位置：文江博客话题详情

C/C++ 中的累积正态分布函数

发布于 2024-08-23 00:39:39 字数 112 浏览 5 评论 0原文

我想知道数学库中是否内置了统计函数，这些数学库是标准 C++ 库（如 cmath）的一部分。如果没有，你们能推荐一个具有累积正态分布函数的好的统计库吗？

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

原文

分享到QQ

分享到微博

如果你对这篇内容有疑问，欢迎到本站社区发帖提问参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

发布评论

需要登录才能够评论，你可以免费注册一个本站的账号。

听不够的曲调 2024-08-30 00:39:39

没有直接函数。但由于高斯误差函数及其补函数与正态累积分布函数有关（请参见此处，或此处）我们可以使用实现的c函数erfc（互补误差函数）：

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

考虑 erfc(x) = 1-erf(x) 与 M_SQRT1_2 = √0,5 的关系。

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

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, or here) we can use the implemented c-function erfc (complementary error function):

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

Which considers the relation of erfc(x) = 1-erf(x) with M_SQRT1_2 = √0,5.

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

回复收藏 0 原文

猫腻 2024-08-30 00:39:39

下面是用 14 行代码实现的累积正态分布的独立 C++ 实现。

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";
}

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

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";
}

回复收藏 0 原文

秋凉 2024-08-30 00:39:39

根据在我之前回答的人的建议，我想出了如何使用 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;
}

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):

#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;
}

回复收藏 0 原文

浅唱ヾ落雨殇 2024-08-30 00:39:39

Boost 与标准一样好：D 在这里：提升数学/统计。

回复收藏 0 原文

年华零落成诗 2024-08-30 00:39:39

这里给出的普通 CDF 的实现是单精度近似值，已将浮点型替换为双精度型，因此仅精确到 7 或 8 位有效值（十进制）数字。
有关 Hart 双精度近似的 VB 实现，请参见 West 的图 2 Better累积正态函数的近似。

编辑：我将 West 的实现翻译成 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;
}

请注意，我已将表达式重新排列为更熟悉的级数和连分数近似形式。 West 代码中的最后一个幻数是 2π 的平方根，我通过利用恒等式 acos(0) = ½ π 将其推迟到第一行的编译器。
我已经三次检查了魔法数字，但总有可能我输错了一些东西。如果发现错别字，请评论！

John Cook 在他的答案中使用的测试数据的结果让

 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

我感到一些小小的安慰，因为他们同意 Mathematica 结果给出的所有数字。

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.

Edit: My translation of West's implementation into 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!

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.

回复收藏 0 原文

拥抱没勇气 2024-08-30 00:39:39

来自 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;
}

From NVIDIA CUDA samples:

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;
}

回复收藏 0 原文

难理解 2024-08-30 00:39:39

来自 https://en.cppreference.com/w/cpp/numeric/math /erfc

正态 CDF 的计算公式如下：

#include ;
#include ;
#include ;
使用命名空间 std；

double normalCDF(double x) // Phi(-∞, x) 又名 N(x)
{
    返回 erfc(-x / sqrt(2))/2;
}

在分母中使用 2.0 而不是 2 有助于得到小数而不是整数。

希望有帮助。

From https://en.cppreference.com/w/cpp/numeric/math/erfc

Normal CDF can be calculated as below:

#include <iostream>
#include <cmath>
#include <iomanip>
using namespace std;

double normalCDF(double x) // Phi(-∞, x) aka N(x)
{
    return erfc(-x / sqrt(2))/2;
}

Using 2.0 instead of 2 in the denominator helps in getting decimals instead of integers.

Hope that helps.

回复收藏 0 原文

小…楫夜泊 2024-08-30 00:39:39

由于这个问题是在大约 13 年前提出的，而且答案已经相当过时了，到目前为止，为了计算正态分布的 cdf，我们可以使用 boost 库，可以从这里下载 https://www.boost.org/。安装最新版本后，#include 任何发行版，例如#include“boost/math/distributions/normal.hpp”，您可以直接使用 cdf。请记住使用命名空间 boost::math。
您可以参考此链接以获取更多参考： https:// www.boost.org/doc/libs/1_80_0/boost/math/distributions.hpp

回复收藏 0 原文

~没有更多了~