我想知道标准C++库（例如cmath）中是否内置了统计函数。如果没有，你们能否推荐一个好的统计库，其中包括累积正态分布函数？
更具体地说，我正在寻找使用/创建累积分布函数。
在C/C++中的累积正态分布函数
59
- Tyler Brock
1
3如果您只需要正态分布的累积分布函数(CDF)，为什么不自己实现呢？它并没有什么神奇的地方，因此实现起来很简单。 - Hannes Ovrén
8个回答
39
这是一个14行代码的独立C++实现累积正态分布的函数。具体内容可参考以下链接：http://www.johndcook.com/cpp_phi.html。请注意保留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";
}
- John D. Cook
3
请注意，这是一个单精度近似值，但下面有一个双精度实现，由 thus spake a.k. 提供。https://dev59.com/p3E95IYBdhLWcg3wY9D6#23119456 - David
1该函数计算
[pnorm](https://cosmosweb.champlain.edu/people/stevens/webtech/R/Chapter-6-R.pdf)，其中mean=0，sd=1。我能否指定自定义的mean和sd？如果是这样，如何修改该函数？ - Suman Khanal1@suman-khanal：你需要进行z变换。只需在函数头中添加参数mean和sd，并在行int sign = 1之前将x =（x - mean）/ sd即可。 - jwdietrich
13
在得到之前回答者的建议后，我使用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;
}
- Tyler Brock
5
7哎呀……不要使用
pow 函数，使用霍纳法则。在这个问题被修正之前我会投反对票（请通知我）。 - Alexandre C.5此代码将失去精度。霍纳法则更稳定（并且更快）。 - Alexandre C.
1为什么不直接使用
double pK3 = K*K*K 等等呢？ - Daniel Bonetti据我所知，但
pow被定义为宏，我认为这是为了允许良好的实现优化常见的幂，比如2和3。因此不要太快放弃pow！ - Aaron McDaid感谢提供代码！我猜想这是正态分布的泰勒展开式，然后对得到的多项式进行积分。总之，我在z范围为-3.8到+3.8，增量为0.01的情况下使用Boost库测试了该代码，并且绝对差值之和abs（boost-cnd_manul）的数量级约为10^-6。 - macroland
11
这里给出的正态分布函数的实现是单精度近似值，将
有关Hart的双精度近似值的VB实现，请参见West的更好的累积正常函数近似第2图。
float替换为double，因此只能精确到7或8个有效数字（十进制）。有关Hart的双精度近似值的VB实现，请参见West的更好的累积正常函数近似第2图。
编辑：我将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结果所给出的所有数字。
- thus spake a.k.
5
这与erfc有何区别？ - Johan Lundberg
这将取决于erfc的精度保证。参数和一半平方根的乘积肯定会略微舍入，这可能会传播到最终值。据称Hart算法对于每个参数都可以精确到双精度，尽管我没有独立验证过。无论如何，两者都比单精度近似要好得多，其中float被替换为double！ - thus spake a.k.
@tantrev，恐怕没有简单的方法将正态分布函数的级数和连分数逼近转换为$t$分布的级数和连分数逼近。 - thus spake a.k.
谢谢@thusspakea.k！我想快速近似计算相同N大小的约1E6个p值的希望是愚蠢的。 - tantrev
11
Boost就像标准库一样好:D，以下是链接：boost数学/统计。
- Hassan Syed
5
有一个标准的内置吗？ - Tyler Brock
不，标准库目前还没有。 - dirkgently
正态分布，是的我这么认为。或者你是在谈论内置于标准库中的情况吗？如果是后者，那么不是。 - Hassan Syed
使用boost的3个原因：1）您可以使用boost附带的实用程序复制所选库。2）这些库通常只是头文件。3）大多数库非常好用，并且在使用C++编译器的任何地方都能正常工作。 - Hassan Syed
2不，Boost文档页面永远不会消失。 - ABCD
9
来自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;
}
版权所有 ©1993-2012 NVIDIA公司。保留所有权利。
这段文字是关于版权信息的，来自于NVIDIA公司的一个网页链接。
- serbaut
1
来自 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; }
将分母中的2改为2.0有助于获得小数，而不是整数。
希望这有所帮助。
- Manohar Reddy Poreddy
1
自从13年前提出这个问题以来，现在的答案已经过时了。要计算正态分布的cdf，我们可以使用boost库，该库可以从https://www.boost.org/下载。一旦您安装了最新版本，就可以包含任何分布，例如#include "boost/math/distributions/normal.hpp"，然后可以直接使用cdf。记得使用命名空间boost::math。
您可以参考此链接进行进一步参考：https://www.boost.org/doc/libs/1_80_0/boost/math/distributions.hpp。
- Khanh
1
网页内容由stack overflow 提供, 点击上面的可以查看英文原文，
原文链接
原文链接