標準正態分布函數的快速計算方法

標準正態分布的分布函數 $\Phi(x)$ 可以說是"數據分析師"統計計算中非常重要的一個函數，基本上有正態分布的地方都或多或少會用上它。在一些特定的問題中，我們"數據分析師"需要大量多次地計算這個函數的取值，比如我經常需要算正態分布與另一個隨機變量之和的分布，這時候就需要用到數值積分，而被積函數就包含 $\Phi(x)$。如果 $Z\sim N(0,1), X\sim f(x)$，$f$ 是 $X$ 的密度函數，那么 $Z+X$ 的分布函數就是

我們"數據分析師"知道，$\Phi(x)$ 沒有簡單的顯式表達式，所以它需要用一定的數值方法進行計算。在大部分的科學計算軟件中，計算的精度往往是第一位的，因此其算法一般會比較復雜。當這個函數需要被計算成千上萬次的時候，速度可能就成為了一個瓶頸。

當然有問題就會有對策，一種常見的做法是略微放棄一些精度，以換取更簡單的計算。在大部分實際應用中，一個合理的誤差大小，例如 $10^{-7}$，一般就足夠了。在這篇文章中，給大家介紹兩種簡單的方法，它們都比R中自帶的 pnorm() 更快，且誤差都控制在 $10^{-7}$ 的級別。

第一種辦法來自于經典參考書 Abramowitz and Stegun: Handbook of Mathematical Functions 的 公式 26.2.17 。其基本思想是把 $\Phi(x)$ 表達成正態密度函數 $\phi(x)$ 和一個有理函數的乘積。這種辦法可以保證誤差小于 $7.5\times 10^{-8}$，一段C++實現可以在 這里 找到。（代碼中的常數與書中的略有區別，是因為代碼是針對誤差函數 $\mathrm{erf}(x)$ 編寫的，它與 $\Phi(x)$ 相差一些常數）

我們來對比一下這種方法與R中 pnorm() 的速度，并驗證其精度。

library(Rcpp)
sourceCpp("test_as26217.cpp") x = seq(-6, 6, by = 1e-6) system.time(y <- pnorm(x)) ## user  system elapsed ## 1.049   0.000   1.048 system.time(asy <- r_as26217ncdf(x)) ## user  system elapsed ## 0.293   0.019   0.311 max(abs(y - asy)) ## [1] 6.968772e-08

可以看出，A&S 26.2.17 的速度大約是 pnorm() 的三倍，且誤差也在預定的范圍里，是對計算效率的一次巨大提升。

那么還有沒有可能更快呢？答案是肯定的，而且你其實已經多次使用過這種方法了。怎么，不相信？看看下面這張圖，你就明白了。

沒錯，這種更快的方法其實就是兩個字：查表。它的基本想法是，我們預先計算好一系列的函數取值 $(x_i,\Phi(x_i))$，然后當我們需要計算某個點 $x_0$ 時，就找到離它最近的兩個點 $x_k$ 和 $x_{k+1}$，再用線性插值的方法得到 $\Phi(x_0)$ 的近似取值：

什么？覺得這個方法太簡單了？先別急，這里面還有不少學問。之前我們"數據分析師"說了，我們需要保證這種方法的誤差不超過 $\epsilon=10^{-7}$，因此就需要合理地選擇預先計算的點。由于 $\Phi(-x)=1-\Phi(x)$，我們暫且只需要考慮 $x$ 為正的情況。如果讓 $x_i = ih,i=0,1,\ldots,N$，那么對函數 $f$ 進行線性插值的誤差將不超過（ 來源 ）

其中 $\Vert f’’ \Vert_{\infty}$ 是函數二階導絕對值的最大值。對于正態分布函數來說，它等于 $\phi(1)\approx 0.242$。于是令 $E(x)=10^{-7}$，我們就可以解出 $h\approx 0.001818$。最后，只要 $x_N>5.199$，即 $N\ge 2860$ 并另所有 $x>x_N$ 的取值等于1，就可以保證整個實數域上 $\Phi(x)$ 的近似誤差都不超過 $10^{-7}$。

這種簡單方法的實現我放在了 Github 上 ，源程序和測試代碼也可以在文章最后找到。下面給出它的表現：

library(Rcpp)
sourceCpp("test_fastncdf.cpp") x = seq(-6, 6, by = 1e-6) system.time(fasty <- r_fastncdf(x)) ## user  system elapsed ## 0.043   0.024   0.066 max(abs(y - fasty)) ## [1] 9.99999e-08

與之前的結果相比，相當于速度是 pnorm() 的15倍！

我們似乎一直以為，在計算機和統計軟件普及以后，一些傳統的做法就會慢慢被淘汰，例如現在除了考試，或許大部分的時間我們都是在用軟件而不是正態概率表。從教學與實際應用的角度來看，這種做法是 應該進行推廣和普及的 ，但這也不妨礙我們從一些“舊知識”中汲取營養。關于這種大巧若拙的做法的故事還有很多，比如廣為流傳的 這一則 。在計算資源匱乏的年代，數據科學家"數據分析師"們想出了各種巧妙的辦法來解決他們遇到的各種問題?，F如今計算機的性能已經遠不是當年可以媲跡，但前人的很多智慧卻依然穿透了時間來為現在的我們提供幫助，不得不說這也是一種緣分吧。