标准正态分布函数的快速计算方法
标准正态分布的分布函数 $\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倍!
我们似乎一直以为,在计算机和统计软件普及以后,一些传统的做法就会慢慢被淘汰,例如现在除了考试,或许大部分的时间我们都是在用软件而不是正态概率表。从教学与实际应用的角度来看,这种做法是 应该进行推广和普及的 ,但这也不妨碍我们从一些“旧知识”中汲取营养。关于这种大巧若拙的做法的故事还有很多,比如广为流传的 这一则 。在计算资源匮乏的年代,数据科学家"数据分析师"们想出了各种巧妙的办法来解决他们遇到的各种问题。现如今计算机的性能已经远不是当年可以媲迹,但前人的很多智慧却依然穿透了时间来为现在的我们提供帮助,不得不说这也是一种缘分吧。
数据分析咨询请扫描二维码
CDA数据分析师在中国航信高科技产业园进行了面向测试度量的数据分析培训课程,培训人数近2 ...
2024-05-01CDA数据分析师走进深圳迈瑞生物医疗电子股份有限公司,在迈瑞总部展开了为期两天的培训,本次课程参训人员线上及线下近百人, ...
2024-05-01CDA数据分析师在合肥市对合肥阳光新能源科技有限公司开展了为期8天的企业内训。 合肥阳光新能源科技 ...
2024-05-01CDA数据分析师走进海尔大学,进行了《数据治理与数据中台建设的道与术》专题培训,培训现场爆满,近百人参加了此次培训。 ...
2024-05-01在中国银行苏州分行培训中心开始数据分析师培训,此次培训课程共10天内容,包括Excel、MySQL、概率论与数理统计、SPSS等内容, ...
2024-05-01从实际的业务需求出发,结合行业的典型应用特点,围绕实际的商业问题,探讨数据挖掘、机器学习模型在金融领域的应用,包括获客、信用评分、细分画像、交叉销售、反欺诈、违规识别、时序预测、运筹优化、流程挖掘九个方面,形成 ...
2024-05-01本次培训课程为线上+线下的模式,由于学员编程能力不一、部分学员没有编程基础,故提供统计学、python基 ...
2024-05-01华夏银行信用卡中心-机器学习培训 1、课程亮点 取材于业界一流企业和顶级咨询公司的行业实践;已经被证明是人人 ...
2024-05-01主 题:数据中台建设及数据分析应用主题分享 1. 数据中台市场洞察 2. 主流数据中台产品比较 3. 某企业数据中 ...
2024-05-01围绕“数据驱动”战略,全力打造我行 300 人数字化人才梯队,着力培养数字化管理人才、大数据专业团队 ...
2024-05-01在当今数据驱动的商业环境中,数据分析成为了企业决策的重要依据。通过对大量数据的收集、处理和分析,企业能够更好地理解市场 ...
2024-04-29在人工智能(AI)的世界里,提示词(Prompt)是一种强大的工具,它能够引导AI按照用户的需求产生特定的输出。本文将深入探讨AI ...
2024-04-29CDA立足未来职场,拓展前沿视野——对外经贸大学保险学院举办“三全育人大讲堂”分享行业最新动态。 ...
2024-04-294月2日,CDA数据分析师创始发起人兼协会理事长赵坚毅博士受邀在浙江万里学院举办了一场以“数字化能力在职场中的作用” ...
2024-04-29随机森林(Random Forests)现在机器学习中比较火的一个算法,是一种基于Bagging的集成学习方法,能够很好地处理分类和回归的问 ...
2022-12-23方差分析是数据分析中常用的一种统计分析方法,接下来让我们简单了解一下方差分析的基本思想和原理吧。 方差分析(Analysis ...
2022-12-23来源:关于数据分析与可视化 关于streamlit-aggrid 数据排序 表格样式的调整 数据 ...
2022-08-03作者:麦叔 定义 「把上面晦涩的概念汇成一句话就是:」 ❝ 回调函数就是一个被作为参 ...
2022-08-03现今,高学历人群日益增多,物以稀为贵的高学历光环淡去。无论本科生还是研究生,甚至博士生,求职竞争力都大不如前,就业压力越来越大。
2022-06-01某家企业10个人面试,有9个本科生……如何脱颖而出,除得体的举止和良好的沟通力外,证书成重要筹码,这也是很多人考证的关键所在。
2022-04-14