Tyl*_*ock 47 c c++ math statistics distribution
我想知道是否有数学库内置的统计函数,它们是标准C++库的一部分,如cmath.如果没有,你们可以推荐一个具有累积正态分布函数的好的统计库吗？提前致谢.
更具体地说,我希望使用/创建累积分布函数.
JFS*_*JFS 39
Theres没有直接的功能.但由于高斯误差函数及其互补函数与正态累积分布函数有关(见这里或这里),我们可以使用实现的c函数erfc(互补误差函数):
double normalCDF(double value)
{
return 0.5 * erfc(-value * M_SQRT1_2);
}
考虑erfc(x) = 1-erf(x)与M_SQRT1_2=√0,5 的关系.
我用它进行统计计算,效果很好.无需使用系数.
Joh*_*ook 32
这是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";
}
根据在我之前回答的人的建议,我想出了如何使用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;
}
这里给出的普通CDF的实现是已经替换的单精度近似,因此仅精确到7或8个有效(十进制)数字.
对于Hart的双精度近似的VB实现,请参见West 对累积正态函数的更好近似的图2 .floatdouble
编辑:我将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π的平方根,我通过利用identity 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结果的所有数字这一事实中得到一些小小的安慰.
来自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 Corporation。版权所有。