python - 以数值稳定的方式评估 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)1logsumexp

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