使用没有数值下溢/上溢的 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²)，即使这样也会下溢：对于float64，PDF 已经太小了 z=±39。