
標(biāo)準(zhǔn)正態(tài)分布函數(shù)的快速計(jì)算方法
標(biāo)準(zhǔn)正態(tài)分布的分布函數(shù) $\Phi(x)$ 可以說是"數(shù)據(jù)分析師"統(tǒng)計(jì)計(jì)算中非常重要的一個(gè)函數(shù),基本上有正態(tài)分布的地方都或多或少會(huì)用上它。在一些特定的問題中,我們"數(shù)據(jù)分析師"需要大量多次地計(jì)算這個(gè)函數(shù)的取值,比如我經(jīng)常需要算正態(tài)分布與另一個(gè)隨機(jī)變量之和的分布,這時(shí)候就需要用到數(shù)值積分,而被積函數(shù)就包含 $\Phi(x)$。如果 $Z\sim N(0,1), X\sim f(x)$,$f$ 是 $X$ 的密度函數(shù),那么 $Z+X$ 的分布函數(shù)就是
我們"數(shù)據(jù)分析師"知道,$\Phi(x)$ 沒有簡(jiǎn)單的顯式表達(dá)式,所以它需要用一定的數(shù)值方法進(jìn)行計(jì)算。在大部分的科學(xué)計(jì)算軟件中,計(jì)算的精度往往是第一位的,因此其算法一般會(huì)比較復(fù)雜。當(dāng)這個(gè)函數(shù)需要被計(jì)算成千上萬次的時(shí)候,速度可能就成為了一個(gè)瓶頸。
當(dāng)然有問題就會(huì)有對(duì)策,一種常見的做法是略微放棄一些精度,以換取更簡(jiǎn)單的計(jì)算。在大部分實(shí)際應(yīng)用中,一個(gè)合理的誤差大小,例如 $10^{-7}$,一般就足夠了。在這篇文章中,給大家介紹兩種簡(jiǎn)單的方法,它們都比R中自帶的 pnorm() 更快,且誤差都控制在 $10^{-7}$ 的級(jí)別。
第一種辦法來自于經(jīng)典參考書 Abramowitz and Stegun: Handbook of Mathematical Functions 的 公式 26.2.17 。其基本思想是把 $\Phi(x)$ 表達(dá)成正態(tài)密度函數(shù) $\phi(x)$ 和一個(gè)有理函數(shù)的乘積。這種辦法可以保證誤差小于 $7.5\times 10^{-8}$,一段C++實(shí)現(xiàn)可以在 這里 找到。(代碼中的常數(shù)與書中的略有區(qū)別,是因?yàn)榇a是針對(duì)誤差函數(shù) $\mathrm{erf}(x)$ 編寫的,它與 $\Phi(x)$ 相差一些常數(shù))
我們來對(duì)比一下這種方法與R中 pnorm() 的速度,并驗(yàn)證其精度。
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() 的三倍,且誤差也在預(yù)定的范圍里,是對(duì)計(jì)算效率的一次巨大提升。
那么還有沒有可能更快呢?答案是肯定的,而且你其實(shí)已經(jīng)多次使用過這種方法了。怎么,不相信?看看下面這張圖,你就明白了。
沒錯(cuò),這種更快的方法其實(shí)就是兩個(gè)字:查表。它的基本想法是,我們預(yù)先計(jì)算好一系列的函數(shù)取值 $(x_i,\Phi(x_i))$,然后當(dāng)我們需要計(jì)算某個(gè)點(diǎn) $x_0$ 時(shí),就找到離它最近的兩個(gè)點(diǎn) $x_k$ 和 $x_{k+1}$,再用線性插值的方法得到 $\Phi(x_0)$ 的近似取值:
什么?覺得這個(gè)方法太簡(jiǎn)單了?先別急,這里面還有不少學(xué)問。之前我們"數(shù)據(jù)分析師"說了,我們需要保證這種方法的誤差不超過 $\epsilon=10^{-7}$,因此就需要合理地選擇預(yù)先計(jì)算的點(diǎn)。由于 $\Phi(-x)=1-\Phi(x)$,我們暫且只需要考慮 $x$ 為正的情況。如果讓 $x_i = ih,i=0,1,\ldots,N$,那么對(duì)函數(shù) $f$ 進(jìn)行線性插值的誤差將不超過( 來源 )
其中 $\Vert f’’ \Vert_{\infty}$ 是函數(shù)二階導(dǎo)絕對(duì)值的最大值。對(duì)于正態(tài)分布函數(shù)來說,它等于 $\phi(1)\approx 0.242$。于是令 $E(x)=10^{-7}$,我們就可以解出 $h\approx 0.001818$。最后,只要 $x_N>5.199$,即 $N\ge 2860$ 并另所有 $x>x_N$ 的取值等于1,就可以保證整個(gè)實(shí)數(shù)域上 $\Phi(x)$ 的近似誤差都不超過 $10^{-7}$。
這種簡(jiǎn)單方法的實(shí)現(xiàn)我放在了 Github 上 ,源程序和測(cè)試代碼也可以在文章最后找到。下面給出它的表現(xiàn):
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
與之前的結(jié)果相比,相當(dāng)于速度是 pnorm() 的15倍!
我們似乎一直以為,在計(jì)算機(jī)和統(tǒng)計(jì)軟件普及以后,一些傳統(tǒng)的做法就會(huì)慢慢被淘汰,例如現(xiàn)在除了考試,或許大部分的時(shí)間我們都是在用軟件而不是正態(tài)概率表。從教學(xué)與實(shí)際應(yīng)用的角度來看,這種做法是 應(yīng)該進(jìn)行推廣和普及的 ,但這也不妨礙我們從一些“舊知識(shí)”中汲取營(yíng)養(yǎng)。關(guān)于這種大巧若拙的做法的故事還有很多,比如廣為流傳的 這一則 。在計(jì)算資源匱乏的年代,數(shù)據(jù)科學(xué)家"數(shù)據(jù)分析師"們想出了各種巧妙的辦法來解決他們遇到的各種問題?,F(xiàn)如今計(jì)算機(jī)的性能已經(jīng)遠(yuǎn)不是當(dāng)年可以媲跡,但前人的很多智慧卻依然穿透了時(shí)間來為現(xiàn)在的我們提供幫助,不得不說這也是一種緣分吧。
數(shù)據(jù)分析咨詢請(qǐng)掃描二維碼
若不方便掃碼,搜微信號(hào):CDAshujufenxi
LSTM 模型輸入長(zhǎng)度選擇技巧:提升序列建模效能的關(guān)鍵? 在循環(huán)神經(jīng)網(wǎng)絡(luò)(RNN)家族中,長(zhǎng)短期記憶網(wǎng)絡(luò)(LSTM)憑借其解決長(zhǎng)序列 ...
2025-07-11CDA 數(shù)據(jù)分析師報(bào)考條件詳解與準(zhǔn)備指南? ? 在數(shù)據(jù)驅(qū)動(dòng)決策的時(shí)代浪潮下,CDA 數(shù)據(jù)分析師認(rèn)證愈發(fā)受到矚目,成為眾多有志投身數(shù) ...
2025-07-11數(shù)據(jù)透視表中兩列相乘合計(jì)的實(shí)用指南? 在數(shù)據(jù)分析的日常工作中,數(shù)據(jù)透視表憑借其強(qiáng)大的數(shù)據(jù)匯總和分析功能,成為了 Excel 用戶 ...
2025-07-11尊敬的考生: 您好! 我們誠摯通知您,CDA Level I和 Level II考試大綱將于 2025年7月25日 實(shí)施重大更新。 此次更新旨在確保認(rèn) ...
2025-07-10BI 大數(shù)據(jù)分析師:連接數(shù)據(jù)與業(yè)務(wù)的價(jià)值轉(zhuǎn)化者? ? 在大數(shù)據(jù)與商業(yè)智能(Business Intelligence,簡(jiǎn)稱 BI)深度融合的時(shí)代,BI ...
2025-07-10SQL 在預(yù)測(cè)分析中的應(yīng)用:從數(shù)據(jù)查詢到趨勢(shì)預(yù)判? ? 在數(shù)據(jù)驅(qū)動(dòng)決策的時(shí)代,預(yù)測(cè)分析作為挖掘數(shù)據(jù)潛在價(jià)值的核心手段,正被廣泛 ...
2025-07-10數(shù)據(jù)查詢結(jié)束后:分析師的收尾工作與價(jià)值深化? ? 在數(shù)據(jù)分析的全流程中,“query end”(查詢結(jié)束)并非工作的終點(diǎn),而是將數(shù) ...
2025-07-10CDA 數(shù)據(jù)分析師考試:從報(bào)考到取證的全攻略? 在數(shù)字經(jīng)濟(jì)蓬勃發(fā)展的今天,數(shù)據(jù)分析師已成為各行業(yè)爭(zhēng)搶的核心人才,而 CDA(Certi ...
2025-07-09【CDA干貨】單樣本趨勢(shì)性檢驗(yàn):捕捉數(shù)據(jù)背后的時(shí)間軌跡? 在數(shù)據(jù)分析的版圖中,單樣本趨勢(shì)性檢驗(yàn)如同一位耐心的偵探,專注于從單 ...
2025-07-09year_month數(shù)據(jù)類型:時(shí)間維度的精準(zhǔn)切片? ? 在數(shù)據(jù)的世界里,時(shí)間是最不可或缺的維度之一,而year_month數(shù)據(jù)類型就像一把精準(zhǔn) ...
2025-07-09CDA 備考干貨:Python 在數(shù)據(jù)分析中的核心應(yīng)用與實(shí)戰(zhàn)技巧? ? 在 CDA 數(shù)據(jù)分析師認(rèn)證考試中,Python 作為數(shù)據(jù)處理與分析的核心 ...
2025-07-08SPSS 中的 Mann-Kendall 檢驗(yàn):數(shù)據(jù)趨勢(shì)與突變分析的有力工具? ? ? 在數(shù)據(jù)分析的廣袤領(lǐng)域中,準(zhǔn)確捕捉數(shù)據(jù)的趨勢(shì)變化以及識(shí)別 ...
2025-07-08備戰(zhàn) CDA 數(shù)據(jù)分析師考試:需要多久?如何規(guī)劃? CDA(Certified Data Analyst)數(shù)據(jù)分析師認(rèn)證作為國內(nèi)權(quán)威的數(shù)據(jù)分析能力認(rèn)證 ...
2025-07-08LSTM 輸出不確定的成因、影響與應(yīng)對(duì)策略? 長(zhǎng)短期記憶網(wǎng)絡(luò)(LSTM)作為循環(huán)神經(jīng)網(wǎng)絡(luò)(RNN)的一種變體,憑借獨(dú)特的門控機(jī)制,在 ...
2025-07-07統(tǒng)計(jì)學(xué)方法在市場(chǎng)調(diào)研數(shù)據(jù)中的深度應(yīng)用? 市場(chǎng)調(diào)研是企業(yè)洞察市場(chǎng)動(dòng)態(tài)、了解消費(fèi)者需求的重要途徑,而統(tǒng)計(jì)學(xué)方法則是市場(chǎng)調(diào)研數(shù) ...
2025-07-07CDA數(shù)據(jù)分析師證書考試全攻略? 在數(shù)字化浪潮席卷全球的當(dāng)下,數(shù)據(jù)已成為企業(yè)決策、行業(yè)發(fā)展的核心驅(qū)動(dòng)力,數(shù)據(jù)分析師也因此成為 ...
2025-07-07剖析 CDA 數(shù)據(jù)分析師考試題型:解鎖高效備考與答題策略? CDA(Certified Data Analyst)數(shù)據(jù)分析師考試作為衡量數(shù)據(jù)專業(yè)能力的 ...
2025-07-04SQL Server 字符串截取轉(zhuǎn)日期:解鎖數(shù)據(jù)處理的關(guān)鍵技能? 在數(shù)據(jù)處理與分析工作中,數(shù)據(jù)格式的規(guī)范性是保證后續(xù)分析準(zhǔn)確性的基礎(chǔ) ...
2025-07-04CDA 數(shù)據(jù)分析師視角:從數(shù)據(jù)迷霧中探尋商業(yè)真相? 在數(shù)字化浪潮席卷全球的今天,數(shù)據(jù)已成為企業(yè)決策的核心驅(qū)動(dòng)力,CDA(Certifie ...
2025-07-04CDA 數(shù)據(jù)分析師:開啟數(shù)據(jù)職業(yè)發(fā)展新征程? ? 在數(shù)據(jù)成為核心生產(chǎn)要素的今天,數(shù)據(jù)分析師的職業(yè)價(jià)值愈發(fā)凸顯。CDA(Certified D ...
2025-07-03