首页 首页 大数据 查看内容

标准正态分布函数的快速计算方法

木马童年 2020-1-31 16:54 71 0

标准正态分布的分布函数 Φ(x)Φ(x) 可以说是统计计算中非常重要的一个函数,基本上有正态分布的地方都或多或少会用上它。在一些特定的问题中,我们需要大量多次地计算这个函数的取值,比如我经常需要算正态分布与另 ...

标准正态分布函数的快速计算方法

标准正态分布的分布函数 Φ(x)Φ(x) 可以说是统计计算中非常重要的一个函数,基本上有正态分布的地方都或多或少会用上它。在一些特定的问题中,我们需要大量多次地计算这个函数的取值,比如我经常需要算正态分布与另一个随机变量之和的分布,这时候就需要用到数值积分,而被积函数就包含 Φ(x)Φ(x)。如果 Z?N(0,1),X?f(x)Z?N(0,1),X?f(x),ff 是 XX 的密度函数,那么 Z+XZ+X 的分布函数就是

P(Z+X≤t)=∫+∞?∞Φ(t?x)f(x)dxP(Z+X≤t)=∫?∞+∞Φ(t?x)f(x)dx

我们知道,Φ(x)Φ(x) 没有简单的显式表达式,所以它需要用一定的数值方法进行计算。在大部分的科学计算软件中,计算的精度往往是第一位的,因此其算法一般会比较复杂。当这个函数需要被计算成千上万次的时候,速度可能就成为了一个瓶颈。

当然有问题就会有对策,一种常见的做法是略微放弃一些精度,以换取更简单的计算。在大部分实际应用中,一个合理的误差大小,例如 10?710?7,一般就足够了。在这篇文章中,给大家介绍两种简单的方法,它们都比R中自带的 pnorm() 更快,且误差都控制在 10?710?7 的级别。

第一种办法来自于经典参考书 Abramowitz and Stegun: Handbook of Mathematical Functions

的公式 26.2.17。其基本思想是把 Φ(x)Φ(x) 表达成正态密度函数 ?(x)?(x) 和一个有理函数的乘积。这种办法可以保证误差小于 7.5×10?87.5×10?8,一段C++实现可以在这里找到。(代码中的常数与书中的略有区别,是因为代码是针对误差函数 erf(x)erf(x) 编写的,它与 Φ(x)Φ(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() 的三倍,且误差也在预定的范围里,是对计算效率的一次巨大提升。

那么还有没有可能更快呢?答案是肯定的,而且你其实已经多次使用过这种方法了。怎么,不相信?看看下面这张图,你就明白了。

标准正态分布函数的快速计算方法

没错,这种更快的方法其实就是两个字:查表。它的基本想法是,我们预先计算好一系列的函数取值 (xi,Φ(xi))(xi,Φ(xi)),然后当我们需要计算某个点 x0x0 时,就找到离它最近的两个点 xkxk 和 xk+1xk+1,再用线性插值的方法得到 Φ(x0)Φ(x0) 的近似取值:

Φ(x0)≈x0?xkxk+1?xkΦ(xi+1)+xk+1?x0xk+1?xkΦ(xi)Φ(x0)≈x0?xkxk+1?xkΦ(xi+1)+xk+1?x0xk+1?xkΦ(xi)

什么?觉得这个方法太简单了?先别急,这里面还有不少学问。之前我们说了,我们需要保证这种方法的误差不超过 ?=10?7?=10?7,因此就需要合理地选择预先计算的点。由于 Φ(?x)=1?Φ(x)Φ(?x)=1?Φ(x),我们暂且只需要考虑 xx 为正的情况。如果让 xi=ih,i=0,1,…,Nxi=ih,i=0,1,…,N,那么对函数 ff 进行线性插值的误差将不超过(来源)

E(x)≤18∥f′′∥∞h2E(x)≤18∥f″∥∞h2

其中 ∥f′′∥∞∥f″∥∞ 是函数二阶导绝对值的最大值。对于正态分布函数来说,它等于 ?(1)≈0.242?(1)≈0.242。于是令 E(x)=10?7E(x)=10?7,我们就可以解出 h≈0.001818h≈0.001818。最后,只要 xN>5.199xN>5.199,即 N≥2860N≥2860 并另所有 x>xNx>xN 的取值等于1,就可以保证整个实数域上 Φ(x)Φ(x) 的近似误差都不超过 10?710?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倍!

我们似乎一直以为,在计算机和统计软件普及以后,一些传统的做法就会慢慢被淘汰,例如现在除了考试,或许大部分的时间我们都是在用软件而不是正态概率表。从教学与实际应用的角度来看,这种做法是应该进行推广和普及的,但这也不妨碍我们从一些“旧知识”中汲取营养。关于这种大巧若拙的做法的故事还有很多,比如广为流传的这一则。在计算资源匮乏的年代,科学家们想出了各种巧妙的办法来解决他们遇到的各种问题。现如今计算机的性能已经远不是当年可以媲迹,但前人的很多智慧却依然穿透了时间来为现在的我们提供帮助,不得不说这也是一种缘分吧。

来自:统计之都

本文出处:,链接:,采用「CC BY-SA 4.0 CN」协议转载学习交流,内容版权归原作者所有,如涉作品、版权和其他问题请联系「我们」处理。

在不久的将来,多智时代一定会彻底走入我们的生活,有兴趣入行未来前沿产业的朋友,可以收藏多智时代,及时获取人工智能、大数据、云计算和物联网的前沿资讯和基础知识,让我们一起携手,引领人工智能的未来!

科学计算 实际应用 Python
0
为您推荐
大数据技术改变城市的运作方式,智慧城市呼之欲出

大数据技术改变城市的运作方式,智慧城市呼

纽奥良虽像大多数城市一样有火灾侦测器安装计划,但直到最近还是要由市民主动申装。纽…...

大数据分析面临生死边缘,未来之路怎么走?

大数据分析面临生死边缘,未来之路怎么走?

大数据分析开始朝着营销落地,尤其像数果智能这类服务于企业的大数据分析供应商,不仅…...

什么是工业大数据,要通过3B和3C来理解?

什么是工业大数据,要通过3B和3C来理解?

核心提示:工业视角的转变如果说前三次工业革命分别从机械化、规模化、标准化、和自动…...

大数据普及为什么说肥了芯片厂商?

大数据普及为什么说肥了芯片厂商?

科技界默默无闻的存在,芯片行业年规模增长到了3520亿美元。半导体给无人驾驶汽车带来…...

大数据技术有哪些,为什么说云计算能力是大数据的根本!

大数据技术有哪些,为什么说云计算能力是大

历史规律告诉我们,任何一次大型技术革命,早期人们总是高估它的影响,会有一轮一轮的…...

个人征信牌照推迟落地,大数据 重新定义个人信用!!

个人征信牌照推迟落地,大数据 重新定义个

为金融学的基础正日益坚实。通过互联网大数据精准记录海量个人行为,进而形成分析结论…...