如何以数值稳定方式计算log(1 - normal_cdf(x))？

这个问题我太熟悉了！当x是很大的正数时，直接执行np.log(1 - norm.cdf(x))确实会因为数值下溢返回-inf，还触发警告——本质原因是norm.cdf(x)太接近1，导致1 - norm.cdf(x)在浮点数精度下直接变成0，自然对数对0就会报错。下面给你几个实用的解决办法，都是数值稳定的：

方法1：用Scipy内置的logsf函数（最推荐） #

Scipy的正态分布模块里专门提供了生存函数的对数——norm.logsf(x)，而生存函数的定义就是1 - CDF(x)，这个函数内部已经做了数值稳定处理，完全不用我们手动计算差值，直接用就行：

import scipy
from scipy.stats import norm
import numpy as np

# 计算log(1 - norm.cdf(10))
result = norm.logsf(10)
print(result)  # 输出约为 -50.00083331943853

这个方法最简单，而且精度有保障，优先用它就对了。

方法2：借助互补误差函数手动推导 #

标准正态分布的CDF可以用误差函数表示：

norm.cdf(x) = 0.5 * (1 + erf(x / sqrt(2)))

由此可以推导出1 - norm.cdf(x) = 0.5 * erfc(x / sqrt(2))，其中erfc是互补误差函数。取对数后就是：

log(1 - norm.cdf(x)) = log(0.5) + log(erfc(x / sqrt(2)))

Scipy的scipy.special模块里有log_erfc函数，专门计算互补误差函数的对数，这样计算也能避免数值下溢：

from scipy.special import log_erfc
import numpy as np

x = 10
result = np.log(0.5) + log_erfc(x / np.sqrt(2))
print(result)  # 和logsf的结果完全一致

方法3：极端大x的近似公式（补充） #

如果x特别大（比如x>5），还可以用近似公式快速计算，精度也足够：
当x趋近于正无穷时，log(1 - norm.cdf(x)) ≈ -x²/2 - log(x) - log(sqrt(2π))
代码示例：

import numpy as np

x = 10
approx_result = -0.5 * x**2 - np.log(x) - 0.5 * np.log(2 * np.pi)
print(approx_result)  # 输出约为 -50.00079663272972，和精确值非常接近

总结一下：优先用norm.logsf(x)，这是最省心且可靠的方案，不需要自己推导公式，Scipy已经帮我们搞定了数值稳定性的问题。

内容的提问来源于stack exchange，提问作者p-value