Scipy stats中二项式PMF输入极端参数时概率和大于1的解决方法

问题描述 #

计算二项式概率质量函数（PMF）时，传入极端参数（极大的n和极小的p）会导致计算结果总和超过1，甚至出现明显不合理的值。示例代码如下：

from scipy.stats import binom
n = 1541096362225563.0
p = 1.0477878413173978e-18
vals = [binom.pmf(n=n, p=p, k=i) for i in range(3)]
print(vals)  # -> [0.9983865609638467, 1.0, 0.060730098812886574]

过滤掉明显错误的1.0后，剩余值总和仍大于1，误差远超可接受的数值精度范围。

原因分析 #

出现这种问题的核心是直接计算二项式PMF在极端参数下的数值不稳定性：

• 当n极大、p极小时，直接计算组合数C(n,k)或高次幂(1-p)^(n-k)会因浮点数精度限制（如双精度浮点数无法精确表示极小的与1的差值）、数值溢出或截断导致结果失真。
• 本例中np ≈ 1.615，符合泊松近似的条件（n→∞，p→0，np=λ为常数），二项分布会收敛于参数为λ=np的泊松分布。

解决办法 #

方法1：使用泊松分布近似二项分布 #

这是最简便且高效的方案，完全适配本例的参数场景：

from scipy.stats import poisson

n = 1541096362225563.0
p = 1.0477878413173978e-18
λ = n * p

# 计算前3个k值的PMF
vals_poisson = [poisson.pmf(k=i, mu=λ) for i in range(3)]
print(vals_poisson)
# 输出示例：[0.1992, 0.322, 0.261]（实际计算值会更精确）
print(sum(vals_poisson))  # 总和接近1，符合概率分布要求

方法2：手动实现对数形式的稳定二项式PMF计算 #

如果需要严格使用二项式分布的精确计算（而非近似），可以通过对数形式避免数值溢出和精度损失：
二项式PMF的对数表达式为：

log(P(k)) = log(C(n,k)) + k·log(p) + (n-k)·log(1-p)

其中组合数的对数可以通过伽马函数的对数（gammaln）计算：log(C(n,k)) = gammaln(n+1) - gammaln(k+1) - gammaln(n-k+1)

代码实现：

import numpy as np
from scipy.special import gammaln

def stable_binom_pmf(n, p, k):
    if k < 0 or k > n:
        return 0.0
    # 计算组合数的对数
    log_comb = gammaln(n + 1) - gammaln(k + 1) - gammaln(n - k + 1)
    # 计算PMF的对数
    log_pmf = log_comb + k * np.log(p) + (n - k) * np.log(1 - p)
    # 转换为原始概率值
    return np.exp(log_pmf)

n = 1541096362225563.0
p = 1.0477878413173978e-18
vals_stable = [stable_binom_pmf(n, p, i) for i in range(3)]
print(vals_stable)
print(sum(vals_stable))  # 总和接近1，精度符合要求

说明 #

• 泊松近似在np ≤ 10时精度极高，完全满足绝大多数场景需求；
• 对数形式计算适用于必须保留二项式分布定义的场景，避免了直接计算大数值带来的精度问题。

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