如何用np.cumsum复现scipy.stats.expon.cdf的输出？

我需要编写一个兼容scipy.stats.kstest的自定义分布函数，因此想了解如何用np.cumsum复现scipy.stats.expon.cdf的输出。目前找不到数值计算CDF的优质指南，大多资源只指向内置方法。我已经实现了custom_exponential_cdf函数，但在部分场景下失效（如下方示例），请问哪里出错了？该如何修改以匹配scipy.stats.expon.cdf的输出？

代码示例 #

import numpy as np
from matplotlib import pyplot as plt
from scipy.stats import kstest, expon

def custom_exponential_cdf(x, lamb):
    x = x.copy()
    x[x < 0] = 0.0
    pdf = lamb * np.exp(-lamb * x)
    cdf = np.cumsum(pdf * np.diff(np.concatenate(([0], x))))
    return cdf

unique_values = np.array([0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.20, 0.21, 0.22, 0.23, 0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 0.3, 0.31, 0.32, 0.33, 0.34, 0.35, 0.36, 0.37, 0.38, 0.39, 0.4, 0.41, 0.42, 0.43, 0.44, 0.45, 0.46, 0.47, 0.48, 0.49, 0.5, 0.51, 0.52, 0.53, 0.54, 0.55, 0.56, 0.57, 0.58, 0.59, 0.6, 0.61, 0.62, 0.63, 0.64, 0.65, 0.66, 0.67, 0.68, 0.69, 0.71, 0.72, 0.73, 0.74, 0.75, 0.76, 0.77, 0.78, 0.79, 0.85, 0.87, 0.89])
counts = np.array([1597, 1525, 1438, 1471, 1311, 1303, 1202, 1147, 1002, 918, 893, 801, 713, 680, 599, 578, 478, 430, 409, 353, 350, 292, 245, 211, 224, 182, 171, 151, 125, 111, 94, 85, 72, 73, 57, 36, 53, 35, 35, 27, 19, 20, 15, 10, 20, 12, 10, 13, 11, 10, 17, 15, 8, 3, 3, 3, 5, 6, 6, 2, 3, 3, 4, 6, 1, 1, 3, 1, 2, 1, 3, 1, 1, 2, 2, 2, 2, 2, 1, 1, 2])
x = np.repeat(unique_values, counts)
lamb = 9.23

fig, ax = plt.subplots()
ax.plot(x, expon.cdf(x, 0.0, 1.0 / lamb))
ax.plot(x, custom_exponential_cdf(x, lamb))

print(kstest(x, custom_exponential_cdf, (lamb,)))
print(kstest(x, "expon", (0.0, 1.0 / lamb)))

打印输出 #

KstestResult(statistic=0.08740741955472273, pvalue=6.857709296861777e-145, statistic_location=0.02, statistic_sign=-1)
KstestResult(statistic=0.0988550670723162, pvalue=2.7098163860110364e-185, statistic_location=0.04, statistic_sign=-1)

绘图输出 #

错误原因 #

你的custom_exponential_cdf函数存在两个核心问题：

• 输入数据未排序：np.cumsum依赖有序的输入区间来累积概率，但你传入的x是重复后的原始数据，并非严格递增序列，导致区间计算混乱。
• 区间积分方法错误：用当前点的PDF值乘以区间宽度的方式，仅适用于均匀采样的微小区间；对于非均匀或重复的点，会重复计算相同区间的概率，导致累积结果偏离真实CDF。

修改思路 #

要复现指数分布的CDF，需：

1. 先对输入的x进行排序去重，确保区间递增；
2. 用更准确的积分方式计算每个区间的概率贡献；
3. 将计算后的CDF值映射回原始输入的所有元素，保证输出长度与输入一致。

修正后的数值计算版本 #

如果是为了学习数值计算CDF的方法，可使用以下代码：

import numpy as np
from matplotlib import pyplot as plt
from scipy.stats import kstest, expon

def custom_exponential_cdf(x, lamb):
    x = np.asarray(x)
    # 处理负数输入
    x_clamped = np.maximum(x, 0.0)
    # 获取唯一值并排序，同时记录原始值对应的索引
    unique_x, idx = np.unique(x_clamped, return_inverse=True)
    # 计算区间宽度：第一个区间从0到第一个unique_x，后续为相邻值的差
    dx = np.diff(np.concatenate(([0], unique_x)))
    # 用区间中点的PDF值计算积分，提升准确性
    mid_points = np.concatenate(([unique_x[0]/2], (unique_x[1:] + unique_x[:-1])/2))
    pdf = lamb * np.exp(-lamb * mid_points)
    # 累积积分得到每个unique_x对应的CDF
    cdf_unique = np.cumsum(pdf * dx)
    # 将CDF映射回原始输入的每个元素
    return cdf_unique[idx]

# 后续代码保持不变
unique_values = np.array([0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.20, 0.21, 0.22, 0.23, 0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 0.3, 0.31, 0.32, 0.33, 0.34, 0.35, 0.36, 0.37, 0.38, 0.39, 0.4, 0.41, 0.42, 0.43, 0.44, 0.45, 0.46, 0.47, 0.48, 0.49, 0.5, 0.51, 0.52, 0.53, 0.54, 0.55, 0.56, 0.57, 0.58, 0.59, 0.6, 0.61, 0.62, 0.63, 0.64, 0.65, 0.66, 0.67, 0.68, 0.69, 0.71, 0.72, 0.73, 0.74, 0.75, 0.76, 0.77, 0.78, 0.79, 0.85, 0.87, 0.89])
counts = np.array([1597, 1525, 1438, 1471, 1311, 1303, 1202, 1147, 1002, 918, 893, 801, 713, 680, 599, 578, 478, 430, 409, 353, 350, 292, 245, 211, 224, 182, 171, 151, 125, 111, 94, 85, 72, 73, 57, 36, 53, 35, 35, 27, 19, 20, 15, 10, 20, 12, 10, 13, 11, 10, 17, 15, 8, 3, 3, 3, 5, 6, 6, 2, 3, 3, 4, 6, 1, 1, 3, 1, 2, 1, 3, 1, 1, 2, 2, 2, 2, 2, 1, 1, 2])
x = np.repeat(unique_values, counts)
lamb = 9.23

fig, ax = plt.subplots()
ax.plot(x, expon.cdf(x, 0.0, 1.0 / lamb), label='scipy expon.cdf')
ax.plot(x, custom_exponential_cdf(x, lamb), label='custom cdf')
ax.legend()

print(kstest(x, custom_exponential_cdf, (lamb,)))
print(kstest(x, "expon", (0.0, 1.0 / lamb)))

高效解析版本 #

如果只是需要兼容kstest的自定义CDF，直接用指数分布的解析公式会更高效准确，完全等价于scipy.stats.expon.cdf：

def custom_exponential_cdf(x, lamb):
    x = np.asarray(x)
    return np.where(x < 0, 0.0, 1 - np.exp(-lamb * x))

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