考虑以下任务：对于任意值 x 和正数 s，计算正态分布随机变量落在以 x 为中心的长度 s 区间的概率。
原则上，这很容易做到：
def normal_inverval_prob(y, s, mean, sd):
return norm.cdf(x=y+s/2.0, loc=mean, scale=sd) - norm.cdf(x=y-s/2.0, loc=mean, scale=sd)
normal_inverval_prob(-3, .2, 1, 1)#2.7438837105055897e-05
normal_inverval_prob(-3, .2, 1, .1)# 0.0
我的问题是最后一行：对于某些值，我得到的概率为零，尽管实际概率是大于零的小数字。这会导致我稍后在我的代码中出现除以零的问题。
事实证明，我可以处理对数概率，所以我重新设计了函数，仅使用对数 cdf 为我提供对数概率：
def normal_inverval_logprob(y, s, mean, sd):
p1 = norm.logcdf(x=y+s/2.0, loc=mean, scale=sd)
p0 = norm.logcdf(x=y-s/2.0, loc=mean, scale=sd)
return p1 + np.log1p(-np.exp(p0 - p1))
np.exp(normal_inverval_logprob(-3, .2, 1, 1))#2.7438837105055897e-05
normal_inverval_logprob(-3, .2, 1, .1)#-765.0831565643776
对于其他值，此对数概率函数会遇到问题：
normal_inverval_logprob(3, .2, 1, .1)
/home/keith/.local/lib/python3.6/site-packages/ipykernel_launcher.py:4: RuntimeWarning: divide by zero encountered in log1p
after removing the cwd from sys.path.
-inf
正如您所料，问题在于，尽管对数 cdfs 不相等，但此时对数 cdfs 差异的 exp 计算结果为 1(一种不同类型的数值下溢问题(：
np.exp(norm.logcdf(2.9, 1, .1) - norm.logcdf(3.1, 1, .1))#1.0
norm.logcdf(3.1, 1, .1) > norm.logcdf(2.9, 1, .1)#True
np.allclose(norm.logcdf(3.1, 1, .1), norm.logcdf(2.9, 1, .1))#True
我不确定如何解决这个问题(或者是否有一些完全不同的方法来实现我的目标(。
一种简单的方法是使用expm1而不是log1p：
return p1 + np.log(-np.expm1(p0 - p1))
如果即使失败了，你也可以用黎曼和(这里只有一个项(近似
：def normal_inverval_prob(y, s, mean, sd):
return norm.pdf(x=y, loc=mean, scale=sd) * s
这将低估尾部;您可以平均区间端点处的值以获得那里的上限。 当然，随着exp(-x2(最终甚至会溢出：PDF已经太小了，无法通过z=±39float64。