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

Question

如何以数值稳定的方式评估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这样的函数可以避免数值under-flow？

Answer 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)

Answer 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

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

Evaluate log(1 - normal_cdf(x)) in a numerically stable way

python

numerical

scipy

underflow

softmax