互补误差函数erfc的近似估计，累积分布函数CDF的计算

        最近由于需要计算累积分布函数(Cumulative Distribution Function/CDF) ，根据公式转化需要计算互补误差函数（erfc函数），在C++11标准库中是有这个函数的，但是比如用VS2010就不能用，木得办法，老师要求的，手动编程erfc函数，计算方式有很多种，第一种就是积分，这个积分不会编写；另一种方法就是近似计算，不同近似方法有不同的近似精度，本次将给出一种双精度的近似方法。

• 概率密度函数

•     pdf转化为关于计算erfc函数的公式

•     erfc函数近似计算

 

• 概率密度函数

• pdf转化为关于计算erfc函数的公式

        关于概率密度函数的计算，一般不直接计算，因为需要计算积分（个人观点），所以转化为计算erfc函数。计算公式可以参考erf函数-中文维基百科，链接没有转到误差函数的可以在里面搜索误差函数，这个链接不用翻墙也可以用。我使用的公式如下图：

转化后为CDF(x)=0.5*erfc(-x / sqrt(2));其他的不再做介绍，可以自己看维基百科相关内容，比百度百科也就好了几百倍。

• erfc函数近似计算

在erf函数-中文维基百科中的数值近似模块有关于erf函数近似计算的内容，精度还可以，我找到了对应的参考文献，发现里面有可以得到双精度的近似计算.h文件。注意，erfc(x)=1-erf(x)，很简单。

具体大家可以找The Art of Scientific Computing第二版、第三版pdf无水印高清版及相应源代码  ,我是在这里下载的，在6.2.2节。书一定要新版本的，老版本的如维基百科中引用的，给出的公式精度并不高。新版本的可以达到双精度，书中没有说精度是多少，我实验发现可以达到1e-16精度。没币下载的我已经截图如下：

 

具体代码，书中代码简写还是怎样，并不能直接运行，接下来给出erf.h的代码（仅仅包括计算erfc的部分，求逆部分书中SQR是求平方根还是求平方我也不知道，可能需要好好研究研究）：

#include <iostream>

struct Erf_approximate {
	//Object for error function and related functions.
	static const int ncof=28;
	static const double cof[28]; //Initialization at end of struct.
	inline double erf_a(double x) {
	//Return erf.x/ for any x.
		if (x >=0.) return 1.0 - erfccheb(x);
		else return erfccheb(-x) - 1.0;
	}
	inline double erfc_a(double x) {
	//Return erfc.x/ for any x.
		if (x >= 0.) return erfccheb(x);
		else return 2.0 - erfccheb(-x);
}
double erfccheb(double z){
	//Evaluate equation (6.2.16) using stored Chebyshev coefficients. User should not call directly.
		int j;
		double t,ty,tmp,d=0.,dd=0.;
		if (z < 0.) throw("erfccheb requires nonnegative argument");
		t = 2./(2.+z);
		ty = 4.*t - 2.;
		for (j=ncof-1;j>0;j--) {
		tmp = d;
		d = ty*d - dd + cof[j];
		dd = tmp;
		}
		return t*exp(-z*z + 0.5*(cof[0] + ty*d) - dd);
	}
};
const double Erf_approximate::cof[28] = {-1.3026537197817094, 6.4196979235649026e-1,
1.9476473204185836e-2,-9.561514786808631e-3,-9.46595344482036e-4,
3.66839497852761e-4,4.2523324806907e-5,-2.0278578112534e-5,
-1.624290004647e-6,1.303655835580e-6,1.5626441722e-8,-8.5238095915e-8,
6.529054439e-9,5.059343495e-9,-9.91364156e-10,-2.27365122e-10,
9.6467911e-11, 2.394038e-12,-6.886027e-12,8.94487e-13, 3.13092e-13,
-1.12708e-13,3.81e-16,7.106e-15,-1.523e-15,-9.4e-17,1.21e-16,-2.8e-17};

接下来我是用VS2019编译的，标准库中有erfc函数，我们对我们写的和标准库中的进行对比，

主函数中的函数如下：

double normalcdf(double value){
	double param = -value / sqrt(double(2));
	cout.precision(16);
	Erf_approximate erf_approximate;
	cout << endl << endl<< 'x'<<param << endl <<"erf_approximate"<< erf_approximate.erfc_a(param) << endl << "erfc" << erfc(param) << endl;
	return 0.5 * erfc(param);
}

结果如下：

可以看到第16位输出不一样，效果还是刚刚的，够用了。

注意，这个x不是我指定的，是我在计算过程中得到的x，所以你使用我的x计算结果可能会有一点差异，你可以自己取x验证。

若有错误，欢迎指正。