使用无数值下溢/溢出的cdf计算概率（Python）

Question

考虑以下任务：对于任意值x和正数S，计算正态分布随机变量落在以x为中心的长度S区间内的概率。p>

原则上这很容易做到：

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

正如您所料，问题在于，尽管对数CDF不相等，但此时对数CDF差异的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

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

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对于z=&；下午；三十九,