SciPy binom.logsf仅为便捷函数？精度问题及配置问询

我原本认为scipy.stats.binom.log*系列函数因为全程在对数空间计算，应该能返回更大范围的数值。但测试显示，binom.logsf和binom.sf在相同输入下都会因精度问题失效，看起来log*类函数只是log(binom.sf())的便捷实现。

问题：这是预期行为吗？还是我遗漏了某个配置步骤？

测试代码 #

for n in range(n_min_max, n_min_max + 20):        
    p_trial_log = binom.logsf(n, N, p_let)
    p_trial = binom.sf(n, N, p_let)
    print(f"args = {n}, {N}, {p_let}")
    print(f"p_trial_log={p_trial_log:.1e}")
    print(f"p_trial=    {p_trial:.4e}")
    print(f"p_trial_exp={exp(p_trial_log):.4e}")
    print("")

测试输出 #

args = 42, 100000, 4e-10
p_trial_log=-5.6e+02
p_trial=    1.2691e-242
p_trial_exp=1.2691e-242

args = 43, 100000, 4e-10
p_trial_log=-5.7e+02
p_trial=    1.1532e-248
p_trial_exp=1.1532e-248

args = 44, 100000, 4e-10
p_trial_log=-5.8e+02
p_trial=    1.0246e-254
p_trial_exp=1.0246e-254

args = 45, 100000, 4e-10
p_trial_log=-6.0e+02
p_trial=    8.9059e-261
p_trial_exp=8.9059e-261

args = 46, 100000, 4e-10
p_trial_log=-6.1e+02
p_trial=    7.5760e-267
p_trial_exp=7.5760e-267

args = 47, 100000, 4e-10
p_trial_log=-6.3e+02
p_trial=    6.3104e-273
p_trial_exp=6.3104e-273

args = 48, 100000, 4e-10
p_trial_log=-6.4e+02
p_trial=    5.1488e-279
p_trial_exp=5.1488e-279

args = 49, 100000, 4e-10
p_trial_log=-6.5e+02
p_trial=    4.1171e-285
p_trial_exp=4.1171e-285

args = 50, 100000, 4e-10
p_trial_log=-6.7e+02
p_trial=    3.2274e-291
p_trial_exp=3.2274e-291

args = 51, 100000, 4e-10
p_trial_log=-6.8e+02
p_trial=    2.4814e-297
p_trial_exp=2.4814e-297

args = 52, 100000, 4e-10
p_trial_log=-7.0e+02
p_trial=    1.8718e-303
p_trial_exp=1.8718e-303

args = 53, 100000, 4e-10
p_trial_log=-7.1e+02
p_trial=    1.3858e-309
p_trial_exp=1.3858e-309

args = 54, 100000, 4e-10
p_trial_log=-7.3e+02
p_trial=    1.0073e-315
p_trial_exp=1.0073e-315

args = 55, 100000, 4e-10
p_trial_log=-7.4e+02
p_trial=    7.2134e-322
p_trial_exp=7.2134e-322

args = 56, 100000, 4e-10
p_trial_log=-inf
p_trial=    0.0000e+00
p_trial_exp=0.0000e+00

args = 57, 100000, 4e-10
p_trial_log=-inf
p_trial=    0.0000e+00
p_trial_exp=0.0000e+00

args = 58, 100000, 4e-10
p_trial_log=-inf
p_trial=    0.0000e+00
p_trial_exp=0.0000e+00

args = 59, 100000, 4e-10
p_trial_log=-inf
p_trial=    0.0000e+00
p_trial_exp=0.0000e+00

args = 60, 100000, 4e-10
p_trial_log=-inf
p_trial=    0.0000e+00
p_trial_exp=0.0000e+00

args = 61, 100000, 4e-10
p_trial_log=-inf
p_trial=    0.0000e+00
p_trial_exp=0.0000e+00

解答 #

这是预期行为，并非遗漏了配置步骤，核心原因是双精度浮点数的数值范围限制：

1. 双精度浮点数能表示的最小正数值约为2.2e-308，当概率值小于这个阈值时，会直接下溢为0；对应的对数形式则会变成-inf。
2. binom.logsf的设计初衷是在概率值未达到下溢阈值时，避免直接计算sf再取对数带来的精度损失——全程在对数空间计算能保留更多有效数字。但当概率小到超出浮点数的表示极限时，不管是sf还是logsf都无法绕过硬件的数值限制，最终都会出现下溢。
3. 你观察到logsf和log(sf)结果一致，是因为当数值还未下溢时，logsf的计算精度确实和直接取对数相当，但这并不代表它只是简单的封装；而当下溢发生后，两者自然都会返回-inf和0。

如果需要处理这种极端参数下的二项分布尾概率，可以考虑：

• 使用对数形式的近似公式（比如泊松近似的对数版本，因为当N很大、p极小时，二项分布可近似为泊松分布Poisson(lambda=N*p)）
• 借助专门处理极端数值的符号计算库，或者使用任意精度浮点数进行计算。

内容的提问来源于stack exchange，提问作者Bob