使用无数值 underflow/overflow 的 cdf 计算概率（在 Python 中）

Question

考虑以下任务：对于任意值 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

如您所料，问题是此时 log cdfs 差异的 exp 求值为 1（另一种数值下溢问题），尽管 log cdfs 不相等：

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

我不确定如何解决这个问题（或者是否有一些完全不同的方法来实现我的目标）。

Answer 1

一个简单的方法是使用 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(-x²) 最终甚至会下溢：PDF 对于 float64 by z[=24= 已经太小了]=±39.

使用无数值 underflow/overflow 的 cdf 计算概率（在 Python 中）

calculating probabilities using cdf without numerical underflow/overflow (in Python)

python

precision

probability