正态累积分布函数的近似计算方法

1. 正态分布的提出

一些学者将发现正态分布的功劳归于德·莫伊夫(de Moivre)，他于1738年在其著作《The Doctrine of Chances》第二版中发表了对二项式 $(a+b)^n$ 展开系数的研究。De Moivre证明，在这个展开的中项有近似的大小 $2^{n}/{\sqrt {2\pi n}}$ ，如果 $m$ 或 $n/2$ 是充分大，那么与中间项间隔 $\ell$ 的项的幅度与中间项幅度之比的对数为 $-2\ell^2/n$ 。虽然这个结论可以被解释为正态概率分布的第一个模糊表达，但Stigler指出，de Moivre本人并没有把他的结果解释为二项式系数的近似规则以外的任何东西，特别是de Moivre缺乏概率密度函数的概念。

1823年，高斯出版了他的专著《Theoria combinationis observationum erroribus minimis obnoxiae》，在其中，他介绍了几个重要的统计概念，如最小二乘方法、最大似然方法和正态分布。高斯用 $M,M',M'',\dots$ 来表示某个未知量 $V$ 的测量值，并寻求这个量的“最有可能”的估计值：使 $\varphi (M-V)\cdot\varphi (M'-V)\cdot\varphi (M''-V)\cdot\dots$ 最大化，得到观测的实验结果。在他的符号中 $\varphi \Delta$ 表示测量误差的概率密度函数。由于不知道函数 $\varphi$ 的表达式，高斯将该方法简化为众所周知的形式：测量值的算术平均值。从该原则出发，高斯证明可将算术平均值合理地作为参数估计的唯一分布是误差服从正态分布：

$\varphi {\mathit {\Delta }}={\frac {h}{\sqrt{ \pi} }}\,e^{-\mathrm {h}^2 \Delta ^2},$

其中 $h$ 为“观测精度的度量”。用这个正态分布作为实验中误差的一般模型，高斯公式现在被称为非线性加权最小二乘法。

虽然高斯是第一个提出正态分布的人，但Pierre-Simon Laplace做出的贡献显著。1774年，Laplace首先提出了几个观测结果的聚合问题，尽管他自己的解决方案导致了Laplace分布。1782年，Laplace首次计算出 $\int e^{-t^2}dt = \sqrt{\pi}$ 的值，为正态分布提供了归一化常数。最后，Laplace在1810年证明了中心极限定理，强调了正态分布的理论重要性。

Pearson首次采用现代符号标准差 $\sigma$ 来表示正态分布，并将该分布正式命名为正态分布。不久之后，1915年，Fisher在正态分布公式中加入了位置参数，表示成现代的方式：

$df= \frac{1}{\sqrt{2\sigma^2\pi}}e^{-\frac{(x-m)^2}{2\sigma^2}}dx$

“标准正态”是指的均值为零、方差为1的正态分布，于20世纪50年代左右开始普遍使用。

正态分布在统计学中的重要性在于：

实际中很多的分布可近似为正态分布。中心极限定理说明，充分多独立随机变量的和近似服从正态分布。
在具有相同方差的所有可能的概率分布中，正态分布在实数上具有最大的不确定性。

2. 正态分布

若随机变量 $X$ 服从均值为 $\mu$ 、方差为 $\sigma^2$ 的正态分布，可以记为 $X\sim N(\mu, \sigma^2)$ ，其概率密度函数(probability density function，pdf)为

$f(x)= \frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(x-m)^2}{2\sigma^2}}$

相应的累积分布函数(cumulative distribution function，CDF)为

$F(x) =\frac{1}{\sqrt{2\pi}\sigma}\int_{-\infty}^x e^{-\frac{(x-m)^2}{2\sigma^2}}dx$

显然，正态CDF的求解决十分困难。

3. 近似求解方法

正态分布的CDF求解决可以采用查表的方式解决，但对于编程而言，查表的方式过于复杂且精度不高，不太适合。下面将介绍Zelent Severo于1964年在“Handbook of mathematical functions with formulas, and mathematical tables”的第26章“Probability Functions”中给出多项式的计算方法，该方法的计算绝对精度可达 $7.5\times 10^{-8}$ 。该方法在张善杰、金建铭著的《特殊函数计算手册》第16章也进行了相应的介绍。

令 $\Phi(x)$ 为标准正态分布的CDF，对于 $x>0$

$\Phi(x) = 1 - \varphi(x)(b_1t + b_2t^2 + b_3t^3 + b_4t^4 + b_5t^5) + \epsilon (x),\qquad t = \frac{1}{1+b_0x}$

式中， $\varphi(x)$ 是标准正态分布的pdf， $b_0 = 0.2316419$ ， $b_1 = 0.319381530$ ， $b_2 = -0.356563782$ ， $b_3 = 1.781477937$ ， $b_4 = -1.821255978$ 和 $b_5 = 1.330274429$ 。

若 $x<0$ ，则相应的CDF为 $1 - \Phi(|x|)$ 。

对于一般意义上的正态分布，可将其转化为标准正态分布后，再行计算。

4. 程序代码

public double cdfGaussian(double x) {double t =  1.0/(1+0.2316419 * Math.abs(x));double b1 = 0.31938153;double b2 = - 0.356563782;double b3 = 1.781477937;double b4 = - 1.821255978;double b5 = 1.330274429;double temp = 1 - Math.exp(- x*x/2)/Math.sqrt(2*Math.PI) * (b1*t + b2 * Math.pow(t, 2) + b3 * Math.pow(t, 3) + b4 * Math.pow(t, 42) + b5 * Math.pow(t, 5));if (x < 0) {temp = 1 - temp;}return temp;
}