以数学形式，
$$ \ begin {align *}
\ textrm {norm_cdf}（z）= 1 -0.3989423 e ^ \ frac {-z ^ 2} {2}（＆1.330274429 t ^ 5-1.821255978 t ^ 4 \\
＆+ 1.781477937 t ^ 3-0.356563782 t ^ 2 + 0.319381530 t）
\ end {align *} $$
其中
$$ t = \ frac {1} {1 + 0.2316419 z} $$
#1 楼
让我引用《 C ++数值食谱》一书（但也适用）：我们假设您已经知道永远不会用这种方式求出多项式：
p=c[0]+c[1]*x+c[2]*x*x+c[3]*x*x*x+c[4]*x*x*x*x;
p=c[0]+c[1]*x+c[2]*pow(x,2.0)+c[3]*pow(x,3.0)+c[4]*pow(x,4.0);
您当然可以方便地拆分，因为括号内的换行符被忽略。记住PEP：“绕开二进制运算符的首选位置是在运算符之后，而不是在运算符之前。”
-t * (0.319381530 + t * (-0.356563782 + t * (1.781477937 + t * (-1.821255978 + 1.330274429 * t))))
另一种选择是将系数保存在列表中：
-t * (0.319381530 + t * (-0.356563782 +
t * (1.781477937 + t * (-1.821255978 + 1.330274429 * t))))
我在进行适当的操作以避免分配和取消分配内存，但这仅在
x是NumPy数组的情况下才有意义。如果要对单个数字求值，则可以改用更好的表达式：coeff = [0, 0.319381530, -0.356563782, 1.781477937, -1.821255978, 1.330274429]
poly = coeff[-1]
for c in coeff[-2::-1]:
poly *= x
poly += c
但是我会坚持使用第一个数字，因为它更通用。必须乘以前置因子：
poly = poly * x + coeff[i]
或者，如果您要就位：
return 1 - 0.3989423*math.exp(-z*z/2) * poly
奖金轨道：对于大型数组（超过一千个数字），您可能会受益于使用numexpr中的第一个技术： />
评论
\ $ \ begingroup \ $
我可以看到以任何不鼓励的方式实施它都会带来严重的性能损失，但没有准确性损失。为什么对准确性不利？我在SO上问了它：stackoverflow.com/q/25203873/174365
\ $ \ endgroup \ $
– Emilio M Bumachar
2014年8月8日12:46
\ $ \ begingroup \ $
也许我错过了一些上下文，但是看起来您出于速度的考虑正在编写更长，更丑陋的代码，这似乎根本不是OP所关心的。
\ $ \ endgroup \ $
–raptortech97
2014年8月8日13:20
\ $ \ begingroup \ $
@EmilioMBumachar取决于多项式。但是对于实践中使用的许多多项式，幼稚的方法将导致减去较大的相似大小的值，这会增加绝对误差，因为较大的数字使用浮点数会产生较大的绝对表示误差。
\ $ \ endgroup \ $
– CodesInChaos
14年8月8日在16:21
\ $ \ begingroup \ $
这种评估多项式的方法称为“霍纳方案”（或取决于您居住的地方的方法）。
\ $ \ endgroup \ $
–艾米莉·L。
2014年8月9日在17:03
\ $ \ begingroup \ $
而不是使用xrange（len（coeff）-2，-1，-1），然后通过i索引coeff，我将遍历coeff [-2 ::-1]
\ $ \ endgroup \ $
–弗朗西斯科·库佐（Francisco Couzo）
16-10-19在6:06
#2 楼
事实证明，最近在Math.SE上提出了类似的问题。无需重新设计轮子，而要利用Python的内置功能。 frac {1} {\ sqrt {2 \ pi}} \ int _ {-\ infty} ^ {z} e ^ {-t ^ 2/2} \ dt= \ int _ {-\ infty} ^ { z} Z（t）\ dt
= \ frac {1} {2} \ left（1 + \ mathrm {erf} \ left（\ frac {z} {\ sqrt {2}} \ right）\ right）
= \ frac {1} {2} \ mathrm {erfc} \ left（-\，\ frac {z} {\ sqrt {2}} \ right）$$
因此，您可以仅使用
norm_cdf(z)（自python 2.7起可用）并获得更准确的结果（尤其是对于\ $ z \ $的非常负的值）。 >更好的是，只需使用math.erfc()！评论
\ $ \ begingroup \ $
导入scipy.stats.norm.cdf（）时，我有一个带有try-except的完整代码，以防万一用户没有安装scipy，但math.erfc会很好。
\ $ \ endgroup \ $
–alvas
2014年8月7日在19:08
\ $ \ begingroup \ $
但我想这会淘汰<= py 2.6用户...
\ $ \ endgroup \ $
–alvas
2014年8月7日19:10
\ $ \ begingroup \ $
如果在python≤2.6上运行，为什么不仅仅要求或重新分发SciPy？
\ $ \ endgroup \ $
– 200_success
2014年8月7日19:12
\ $ \ begingroup \ $
或者，只需窃取此erf（）的公共实现即可。
\ $ \ endgroup \ $
– 200_success
2014年8月7日在19:16
\ $ \ begingroup \ $
@ 200_success：A＆S 26.2.17明确表示有效输入为非负。不应使用负值来调用代码，并且在可能发生无效调用的环境中，可能应包括检查。错误界限列为7.5 * 10 ** -8。是否以某种方式为math.erfc指定了精度？我并不是说算法是好是坏，我只是说应该在解决问题的背景下进行权衡。
\ $ \ endgroup \ $
–花岗岩罗伯特
2014年8月8日在13:05
#3 楼
我只讲几位审稿人提到的所谓的“幻数”。有时，当您在从事纯数学工作时，乍一看似乎是一个“幻数”真的不是。数字本身可能只是问题陈述的一部分。我认为问题可以归结为：您能提出一个比数字更具描述性的名称吗？如果有个好名字，您应该使用它。
乍一看，我认为您的电话号码是问题的固有部分。但是当我查看Abramowitz和Stegun时，我看到引用的公式已经为您看起来丑陋的常量命名。名称是
p（您在注释中提到），以及b1至b5。您应该在代码中使用这些名称，因为它们创建了指向原始公式定义的非常清晰的链接。该号码应命名。 （并且一旦代码显示p=0.2316419，则应删除该注释。）顺便说一句，在注释中包含确切的Abramowitz和Stegun参考对您非常有用。
#4 楼
代替使用math.pow，请使用内置的**运算符。您无需在EOL处使用\，因为表达式周围的括号使它可以隐式跨越多行。因此，在完成所有这些更改之后，我得出以下结论：def norm_cdf(z):
""" Use the norm distribution functions as of Gale-Church (1993) srcfile. """
# Equation 26.2.17 from Abramowitz and Stegun (1964:p.932)
t = 1.0 / (1 + 0.2316419*z) # t = 1/(1+pz) , p=0.2316419
probdist = 1.0 - ( (0.319381530 * t)
+ (-0.356563782 * t**2)
+ (1.781477937 * t**3)
+ (-1.821255978 * t**4)
+ (1.330274429 * t**5)) * 0.3989423*math.exp(-z*z/2)
return probdist
我还重新排列了乘法之一，使优先级更加明显和可读。 />
经过这些更改之后，我仍然看到的唯一问题是所有魔术数字。我不知道这些常量是如何得出的，但是如果可以给这些常量赋予有意义的名称，则可能有助于提高可读性。有时候，公式实际上并没有太多有意义的名称。
#5 楼
不确定这是否有帮助，但是您可以轻松定义一个函数来计算给定位置处的多项式的值def evaluate_polynom(pol, x):
return sum(a * math.pow(x, i) for i, a in enumerate(pol))
然后
(0.319381530 * t) + (-0.356563782* math.pow(t,2)) + (1.781477937 * math.pow(t,3)) + (-1.821255978* math.pow(t,4)) + (1.330274429 * math.pow(t,5))
成为：
evaluate_polynom([0, 0.319381530, -0.356563782, 1.781477937, -1.821255978, 1.330274429], t)
评论
\ $ \ begingroup \ $
如果是这样，那么您可以定义另一个函数来调用您描述的代码片段，作为某种辅助方法。
\ $ \ endgroup \ $
– Pimgd
2014年8月7日13:37
\ $ \ begingroup \ $
很显然，我没有足够的知识将其放入具有有意义名称的函数中，但是我什至不确定是否需要这样做。我认为它带来了使它看起来像所涉及的数学对象所需的抽象级别，将其移动到其他位置可能会破坏数学公式。
\ $ \ endgroup \ $
– SylvainD
2014年8月7日13:42
\ $ \ begingroup \ $
我唯一拥有的nitpick可能是将pol重命名为coeff或系数以使其更加清晰，或者可能使用如下所述的内置函数？
\ $ \ endgroup \ $
– C.B.
2014年8月7日15:53
#6 楼
冒着（进一步）激怒版主的风险，我建议您不要遵循任何源自任何数字食谱的数字建议1。万灵药很长的路要走。它没有解决求和值时的舍入问题，而是试图避免该问题。不幸的是，这样做只是部分成功。对于大多数多项式，仍然会有结果相对较差的输入，即使充其量也是如此。如果多项式项生成的值在求和时可能导致数值不稳定，则可以考虑单独生成每个项，然后使用类似Kahan求和的方式对这些项求和。如果您愿意的话，这还会给您带来误差余量以及总和本身。至少在上次我仔细看时，这似乎是评估多项式的最新技术。与使用Horner方案使用双精度浮点数的结果保持大致相同的精度（例如，使用64位double，它提供的精度与使用128位Quad-精度浮点数，但几乎没有通常会带来的速度损失）。像Kahan求和一样，它支持对结果上的误差范围进行计算。
1。作为参考，请参见以下注释：http：//www.stat.uchicago.edu/~lekheng/courses/302/wnnr/nr.htmlhttp：//www.lysator.liu.se/c/num-recipes-in- c.html
我认为比“精彩”更准确的总结是：“我发现数字食谱为一个人提供了足够的信息来使自己陷入困境，因为阅读NR后，人们认为人们可以理解发生了什么事。“
#7 楼
相反：a*x**3 + b*x**2 + c*x = ((a*x + b)*x + c)*x
h / t @ Emily L.，代表“霍纳”：
https://en.wikipedia.org/ wiki / Horner％27s_method
和h / t @Davidmh指出此方法的计算速度/精度有所改进
gal-church在1990年这样引用：
import math
def pnorm(z):
t = 1 / (1 + 0.2316419 * z)
pd = (1 - 0.3989423 *
math.exp(-z * z / 2) *
((((1.330274429 * t - 1.821255978) * t
+ 1.781477937) * t - 0.356563782) * t + 0.319381530) * t)
return pd
该方法可以方便地避免出现t ^ n问题。
引用：
来源：
http://www.aclweb.org/anthology/J93-1004
第28页中的第21页pdf
《计算语言学》第19卷第95页第1页
我可能会“美化”为：
def pnorm(z):
t = 1 / (1 + 0.2316419 * z)
pd = (1 - 0.3989423 * math.exp(-z * z / 2) *
((((1.330274429 * t -
1.821255978) * t +
1.781477937) * t -
0.356563782) * t +
0.319381530) * t )
return pd
如果您查看
Abromowitz and Stegun，数学函数手册
第932页等式26.2.17
引用：
http://people.math.sfu.ca/~cbm/aands/page_932.htm
您会发现：
我们可以从中创建表格给我们：
def pnorm(z):
p = 0.2316419
b1 = 0.319381530
b2 = -0.356563782
b3 = 1.781477937
b4 = -1.821255978
b5 = 1.330274429
t = 1 / (1 + p * z)
pd = (1 - 0.3989423 * math.exp(-z * z / 2) *
((((b5 * t + b4) * t + b3) * t + b2) * t + b1) * t )
return pd
然后从上一页开始； 931您会找到：
Zx = (1/√(2* π))*e(-z*z/2)
在python中： >我们发现（1 /√（2 *π））= 0.3989423
我也喜欢这样：
Zx = (1/math.sqrt(2* math.pi))*math.exp(-z*z/2)
优于：
t * (b1 + t * (b2 + t * (b3 + t * (b4 + t * b5))))
所以，最后，
(((b5 * t + b4) * t + b3) * t + b2) * t + b1) * t
对照操作员检查我的工作
import math
def pnorm(z):
p = 0.2316419
b1 = 0.319381530
b2 = -0.356563782
b3 = 1.781477937
b4 = -1.821255978
b5 = 1.330274429
t = 1 / (1 + p * z)
Zx = (1 / math.sqrt(2 * math.pi)) * math.exp(-z * z / 2)
pd = Zx * t * (b1 + t * (b2 - t * (b3 + t * (b4 - t * b5))))
return (1 - pd)
我期望Horner方法会更快，所以我进行了时间测试，替换为：
import matplotlib.pyplot as plt
import numpy as np
import math
def norm_cdf(z):
""" Use the norm distribution functions as of Gale-Church (1993) srcfile. """
# Equation 26.2.17 from Abramowitz and Stegun (1964:p.932)
t = 1.0 / (1+0.2316419*z) # t = 1/(1+pz) , p=0.2316419
probdist = 1 - 0.3989423*math.exp(-z*z/2) * ((0.319381530 * t)+ \
(-0.356563782* math.pow(t,2))+ \
(1.781477937 * math.pow(t,3)) + \
(-1.821255978* math.pow(t,4)) + \
(1.330274429 * math.pow(t,5)))
return probdist
for z in np.arange (-3,3,0.01):
zf = pnorm(z)
plt.plot(z,zf, c='red', marker = '.', ms=1)
for z in np.arange (-3,3,0.01):
zf = norm_cdf(z)+0.1 #offset 0.1
plt.plot(z,zf, c='blue', marker = '.', ms=1)
plt.show()
plt.pause(0.1)
使其公平，并将np.arrange分辨率提高到0.0001：
#Zx = (1.0 / math.sqrt(2.0 * math.pi)) * math.exp(-z * z / 2.0)
Zx = 0.3989423* math.exp(-z * z / 2.0)
结果，旋转四核AMD 7950 FM2 +带16G内存的硬盘相当困难（尽管正在运行其他多个应用程序）...违背了我的期望：
评论
\ $ \ begingroup \ $
欢迎使用StackExchange代码复习！很酷的历史，感谢分享。
\ $ \ endgroup \ $
–斯蒂芬·劳赫（Stephen Rauch）
17 Mar 20 '17在2:44
\ $ \ begingroup \ $
我也建议您在此处使用MathJax改进公式。如果您不知道该怎么做，我稍后可以进行处理。
\ $ \ endgroup \ $
– Jamal♦
17 Mar 20 '17 at 14:32
\ $ \ begingroup \ $
嗯，我以前从未与Mathjax一起工作过，所以有兴趣看到Jamal的结果；感谢您的欢迎Stephen Rauch
\ $ \ endgroup \ $
– litepresence
17 Mar 20 '17 at 14:41
评论
@ 200_success喜欢公式图像！python中的Jeez注释很难看...另一方面，您可以为您的数学常数做类似C样式的typedef吗？它将使代码变得更加易读。诸如typdef 0.319381530 NORMALITY之类的东西。但您知道，请选择更好的名称和名称。对于Python，您可以将所有常量放在方法的顶部（出于可读性考虑），例如NORMALITY = 0.319381530。无论如何，这是组织数学函数的首选方法。
@ChrisCirefice在Python中没有typedef这样的东西。将系数保存在列表中是一种更具可读性的方法。
@Davidmh我知道Python中没有typedef，这就是为什么我建议只定义一个“常数”，例如NORMALITY = 0.319381530。我个人不喜欢@Davidmh的答案中建议的列表方法：coeff = [1，0.319381530，-0.356563782，1.781477937，-1.821255978，1.330274429]如果有很多常量，请使用coeff [0]，coeff [ 1]等开始在长函数中变得不可读，尤其是多次使用常量的函数。再说一次，只是我的意见。
@ChrisCirefice列表在循环中很有用，系数不具有任何特殊含义，只能在从多项式派生的表达式（导数，积分等）中重用。