如何在Python中基于Poisson变量正确应用Skellam分布？

你之前的核心错误 #

Skellam分布描述的是两个独立泊松随机变量的差值的概率分布，即 ( Z = X - Y )（其中 ( X \sim \text{Poisson}(\mu_X) )，( Y \sim \text{Poisson}(\mu_Y) )），不是把两个泊松事件的概率值直接相减。你之前计算的21.75%是 ( P(X=3) ) 的概率，7.58%是 ( P(Y=2) ) 的概率，直接相减没有任何统计意义。

用Scipy正确实现Skellam分布 #

Scipy的scipy.stats模块已经内置了Skellam分布的实现，直接调用即可完成概率计算、样本生成等操作。

步骤1：导入依赖库

from scipy.stats import skellam
import numpy as np

步骤2：初始化Skellam分布

传入两个泊松分布的均值参数 ( \mu_X ) 和 ( \mu_Y )：

mu_X = 2.6
mu_Y = 0.5
# 初始化Skellam分布对象
skellam_dist = skellam(mu_X, mu_Y)

步骤3：计算概率质量函数（PMF）

PMF用于计算Z取某个具体整数值的概率，比如你之前计算的 ( X=3 )、( Y=2 )，对应 ( Z=1 )，计算 ( P(Z=1) )：

# 计算Z=1时的概率
pmf_z1 = skellam_dist.pmf(1)
print(f"P(Z=1) = {pmf_z1 * 100:.2f}%")

如果需要批量计算多个Z值的概率：

# 计算Z从-2到5的概率
z_values = np.arange(-2, 6)
pmf_values = skellam_dist.pmf(z_values)
for z, pmf in zip(z_values, pmf_values):
    print(f"P(Z={z}) = {pmf * 100:.2f}%")

步骤4：计算累积分布函数（CDF）

CDF用于计算 ( Z \leq \text{某个值} ) 的累积概率，比如 ( P(Z \leq 2) )：

cdf_z2 = skellam_dist.cdf(2)
print(f"P(Z ≤ 2) = {cdf_z2 * 100:.2f}%")

步骤5：生成Skellam分布的样本

如果需要生成符合该分布的随机样本：

# 生成10个随机样本
samples = skellam_dist.rvs(size=10)
print("生成的Skellam分布样本：", samples)

手动验证（可选） #

如果想理解Skellam分布的计算逻辑，可以手动累加所有满足 ( x-y=z ) 的 ( (x,y) ) 组合的概率乘积，以 ( Z=1 ) 为例：

manual_pmf = 0.0
x = 1  # y = x-1 ≥ 0 → x≥1
while True:
    # 计算P(X=x)
    poisson_x = (mu_X ** x * np.exp(-mu_X)) / np.math.factorial(x)
    # 计算P(Y=x-1)
    poisson_y = (mu_Y ** (x-1) * np.exp(-mu_Y)) / np.math.factorial(x-1)
    term = poisson_x * poisson_y
    manual_pmf += term
    # 当项足够小时停止计算
    if term < 1e-10:
        break
    x += 1
print(f"手动计算P(Z=1) = {manual_pmf * 100:.2f}%")

这个结果会和Scipy计算的PMF值几乎一致。

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