以数值稳定的方式评估 log(1 - normal_cdf(x))答案

【问题标题】：Evaluate log(1 - normal_cdf(x)) in a numerically stable way以数值稳定的方式评估 log(1 - normal_cdf(x))
【发布时间】：2018-05-10 23:18:18
【问题描述】：

如何以数值稳定的方式评估log(1 - normal_cdf(x))？这里normal_cdf是标准正态分布的累积分布函数。

例如，在 Python 中：

import scipy 
from scipy.stats import norm

np.log(1 - norm.cdf(10))

给出-inf 和RuntimeWarning: divide by zero encountered in log，因为norm.cdf(10) 几乎等于1。有没有像logsumexp 这样可以避免数值下溢的函数？

【问题讨论】：

因为它不是产品的对数，所以不能利用对数函数的属性。检查此答案以了解一些有趣的事实math.stackexchange.com/questions/75074/…
也许norm.cdf(10)=1-norm.cdf(-10) 所以np.log(1 - norm.cdf(10))=np.log(norm.cdf(-10))

标签： python scipy numerical softmax underflow

【解决方案1】：

由于正态分布关于 0 对称，我们有

1 - F(x) = P(X > x)
         = P(X < -x)
         = F(-x)

因此

np.log(1 - norm.cdf(10)) = np.log(norm.cdf(-10))
                         = norm.logcdf(-10)

【讨论】：

【解决方案2】：

@HongOoi 利用对称性的建议很棒。但是对于scipy.stats（包括norm）中的任意分布，您可以使用logsf 方法进行精确计算。 sf代表survival function，也就是函数名1 - cdf(x)。

例如，

In [25]: import numpy as np

In [26]: from scipy.stats import norm, gamma

这是norm.logsf 的示例：

In [27]: norm.logsf(3, loc=1, scale=1.5)
Out[27]: -2.3945773661586434

In [28]: np.log(1 - norm.cdf(3, loc=1, scale=1.5))
Out[28]: -2.3945773661586434

下面是gamma.logsf 的示例：

In [29]: gamma.logsf(1.2345, a=2, scale=1.8)
Out[29]: -0.16357333194167956

In [30]: np.log(1 - gamma.cdf(1.2345, a=2, scale=1.8))
Out[30]: -0.16357333194167956

这说明了为什么要使用logsf(x) 而不是log(1 - cdf(x))：

In [35]: norm.logsf(50, loc=1, scale=1.5)
Out[35]: -537.96178420294677

In [36]: np.log(1 - norm.cdf(50, loc=1, scale=1.5))
/Users/warren/miniconda3scipy/bin/ipython:1: RuntimeWarning: divide by zero encountered in log
  #!/Users/warren/miniconda3scipy/bin/python
Out[36]: -inf

【讨论】：