C/C++ 中的累积正态分布函数
声明：本页面是StackOverFlow热门问题的中英对照翻译，遵循CC BY-SA 4.0协议，如果您需要使用它，必须同样遵循CC BY-SA许可，注明原文地址和作者信息，同时你必须将它归于原作者（不是我）：StackOverFlow
原文地址: http://stackoverflow.com/questions/2328258/
Warning: these are provided under cc-by-sa 4.0 license. You are free to use/share it, But you must attribute it to the original authors (not me):
StackOverFlow
Cumulative Normal Distribution Function in C/C++
提问by Tyler Brock
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.
更具体地说，我希望使用/创建一个累积分布函数。
采纳答案by Tyler Brock
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;
}
回答by JFS
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):
没有直接的功能。但由于高斯误差函数及其互补函数与正态累积分布函数有关（参见此处或此处），我们可以使用实现的 c 函数erfc（互补误差函数）：
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.
其中考虑erfc(x) = 1-erf(x)与M_SQRT1_2= √0,5的关系。
I use it for statistical calculations and it works great. No need for using coefficients.
我用它进行统计计算，效果很好。无需使用系数。
回答by John D. Cook
Here's a stand-alone C++ implementation of the cumulative normal distribution in 14 lines of code.
这是一个用 14 行代码实现的累积正态分布的独立 C++ 实现。
http://www.johndcook.com/cpp_phi.html
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";
}
回答by Hassan Syed
Boost is as good as the standard :D here you go: boost maths/statistical.
Boost 和标准一样好 :D 给你：boost maths/statistical。
回答by thus spake a.k.
The implementations of the normal CDF given here are single precisionapproximations that have had floatreplaced with doubleand hence are only accurate to 7 or 8 significant (decimal) figures.
For a VB implementation of Hart's double precisionapproximation, see figure 2 of West's Better approximations to cumulative normal functions.
这里给出的普通 CDF 的实现是单精度近似值，已float被替换，double因此只能精确到 7 或 8 位有效（十进制）数字。
有关 Hart双精度逼近的 VB 实现，请参见 West对累积正态函数的Better approximations 的图 2 。
Edit: My translation of West's implementation into C++:
编辑：我将 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;
}
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!
请注意，我已将表达式重新排列为更熟悉的级数和连分数近似形式。West 代码中的最后一个神奇数字是 2π 的平方根，我在第一行通过利用等式 acos(0) = ½ π 将其推迟给编译器。
我已经对幻数进行了三次检查，但总有可能我打错了一些东西。如果你发现错别字，请评论！
The results for the test data John Cook used in his answer are
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
I take some small comfort from the fact that they agree to all of the digits given for the Mathematica results.
我从他们同意 Mathematica 结果给出的所有数字这一事实中得到一些小小的安慰。
回答by serbaut
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;
}
Copyright 1993-2012 NVIDIA Corporation. All rights reserved.
回答by Manohar Reddy Poreddy
From https://en.cppreference.com/w/cpp/numeric/math/erfc
来自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; }
正常 CDF 可以计算如下：
#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.
在分母中使用 2.0 而不是 2 有助于获得小数而不是整数。
Hope that helps.
希望有帮助。
