在做数据模拟或者蒙特卡洛实验时，经常遇到这样一个需求：手里只有一组观测样本，没有明确的解析分布形式，但希望生成服从这组样本分布的新随机数。这时候经验累积分布函数（Empirical CDF）就派上用场了。直接对经验CDF做逆变换抽样是最直观的思路，不过它只能抽到原始样本点本身；如果想要在样本点之间也能取值，让分布更平滑，就需要插值或者核密度估计的介入。本文把这两种路线都完整梳理一遍，配上可运行的代码，方便直接套用到自己的项目里。
一、直接法：基于经验CDF的逆变换抽样
逆变换抽样的核心思想很简单：如果U是[0,1]上的均匀随机数，那么F⁻¹(U)就服从以F为累积分布函数的分布。对于经验CDF来说，F是一个阶梯函数，它的逆就是把均匀随机数映射到阶梯对应的样本值上。落到工程实现上，最省事的方式是numpy.random.choice配合归一化的权重，或者用numpy.searchsorted在排序后的样本里查找分位点。
先看第一种写法，直接用choice按样本频率抽样：
import numpy as np # 假设这是观测到的样本数据 samples = np.array([1.2, 2.5, 2.5, 3.1, 4.7, 4.7, 4.7, 6.3, 8.9]) # 统计每个唯一值出现的次数，作为抽样权重 unique_values, counts = np.unique(samples, return_counts=True) probabilities = counts / counts.sum() # 从经验分布中抽取10000个新样本 rng = np.random.default_rng(seed=42) new_samples = rng.choice(unique_values, size=10000, p=probabilities) print(new_samples[:10])
这种写法的特点是简单直接，计算开销小，但它有一个明显局限：抽出来的值永远只能是原始数据中出现过的那些点，重复样本越多，被抽中的概率就越大。这本质上是一个离散分布的抽样。如果观测数据本身是连续型的，只是被采样成了有限个点，这种做法会丢失分布的连续性特征。
再看searchsorted的做法，它是逆变换思想的忠实实现：
import numpy as np
def sample_from_empirical_cdf(samples, n, rng):
sorted_samples = np.sort(samples)
u = rng.uniform(size=n)
# 查找每个均匀随机数落在经验CDF的哪个台阶上
indices = np.searchsorted(sorted_samples, u, side='left')
indices = np.clip(indices, 0, len(sorted_samples) - 1)
return sorted_samples[indices]
samples = np.array([1.2, 2.5, 3.1, 4.7, 6.3, 8.9])
rng = np.random.default_rng(seed=42)
draws = sample_from_empirical_cdf(samples, 10000, rng)
print(np.percentile(draws, [25, 50, 75]))
注意这里的searchsorted查找的是样本值数组，严格来说完整的逆变换应该构造分位数函数。不过两者在离散经验分布的场景下是等价的：每个样本点被抽中的概率都是1/N（或者按重复次数加权）。直接法的优点是零假设、不引入任何平滑，完全忠实于数据；缺点是支持集只有样本点本身，无法生成新值。
二、平滑插值法：让经验CDF变成连续可逆的函数
如果希望抽样结果能取到样本点之间的值，一个常用技巧是把经验CDF的分位数函数当作连续函数来插值。具体做法是：把排序后的样本值对应到(0,1)区间上均匀分布的概率点，然后用保形插值（如PchipInterpolator）拟合这条分位数曲线，之后对均匀随机数求值即可得到连续的抽样结果。
为什么推荐Pchip而不是普通的线性插值？因为普通三次样条在数据分布不均匀时容易产生过冲，插值结果可能超出样本取值范围，甚至出现分位数函数非单调的情况，抽出来的数就会违反分布的基本性质。Pchip（单调三次Hermite插值）能保证曲线单调，非常适合用来拟合CDF或其逆函数。代码如下：
import numpy as np
from scipy.interpolate import PchipInterpolator
def smooth_quantile_sampling(samples, n, rng):
sorted_samples = np.sort(samples)
m = len(sorted_samples)
# 构造分位数函数的自变量：用(i+0.5)/m避免端点0和1
p = (np.arange(m) + 0.5) / m
# 用保形插值拟合分位数函数 Q(p)
quantile_func = PchipInterpolator(p, sorted_samples, extrapolate=False)
u = rng.uniform(low=0.5 / m, high=1 - 0.5 / m, size=n)
return quantile_func(u)
samples = np.array([1.2, 2.5, 3.1, 4.7, 6.3, 8.9])
rng = np.random.default_rng(seed=42)
draws = smooth_quantile_sampling(samples, 5000, rng)
print(draws.min(), draws.max())
这段代码有几个细节值得注意。第一，概率点用的是(i+0.5)/m而不是i/(m-1)，因为分位数函数在p趋近0和1时理论上没有定义（经验分布在两端是开区间），取中点位置可以避开端点处的数值不稳定。第二，均匀随机数的取值范围也被限制在同样的开区间内，配合extrapolate=False防止插值器外推出奇怪的结果。第三，如果希望分布两端有更长的尾巴，可以在首尾人为追加几个扩展点，但这属于建模选择，要看具体业务是否合理。
这种平滑插值法的优点是抽样结果连续、分布形态平滑，且严格保持了样本的分位数结构——插值曲线在每个概率点上精确通过原始样本值，所以中位数、四分位数等统计量都能很好还原。缺点是它假设了样本点之间分位数线性或近线性变化，在数据稀疏的区域（比如尾部）可能产生不够真实的插值。
三、进阶方案：核密度估计与resample方法
除了手动插值，scipy还提供了一个更系统的方案：scipy.stats.gaussian_kde。它用核密度估计（KDE）把离散样本变成一条平滑的概率密度曲线，而且自带resample方法可以直接抽样，不需要自己实现逆变换：
import numpy as np
from scipy.stats import gaussian_kde
samples = np.array([1.2, 2.5, 3.1, 4.7, 6.3, 8.9,
4.7, 4.7, 3.1, 5.5, 7.2])
kde = gaussian_kde(samples)
rng = np.random.default_rng(seed=42)
draws = kde.resample(5000, seed=rng)
print(np.percentile(draws, [25, 50, 75]))
KDE的平滑程度由带宽参数控制，默认的Scott规则在多数情况下表现不错。如果觉得结果过于平滑或过于粗糙，可以手动调整：
# 手动指定带宽因子，越小越贴近原始数据 kde_tight = gaussian_kde(samples, bw_method=0.15) draws_tight = kde_tight.resample(5000, seed=rng) # 也可以基于已有KDE缩放带宽 kde_loose = gaussian_kde(samples, bw_method=kde.factor * 2.0)
三种方法的取舍可以总结成一张表：
| 方法 | 能否取到新值 | 平滑程度 | 主要风险 |
|---|---|---|---|
| choice直接抽样 | 不能 | 无平滑 | 结果只能是已有样本点 |
| Pchip插值分位数函数 | 能 | 保形平滑 | 数据稀疏区域插值可能失真 |
| gaussian_kde | 能 | 核函数平滑 | 带宽选择敏感，尾部可能延伸过度 |
从实践角度给几点建议：如果目的只是做bootstrap重采样，直接用choice就够了，没必要引入平滑；如果数据量较大（几百点以上）且希望生成连续的模拟值，Pchip插值分位数函数是性价比最高的方案，实现简单且单调性有保证；如果需要密度函数本身（比如后续要做条件抽样或混合模型），KDE更合适。另外无论用哪种方法，都建议先用Kolmogorov-Smirnov检验（scipy.stats.kstest配合抽样结果与原始样本对比）验证抽样分布和原始分布的一致性，确认方法没有走偏。随机数生成记得统一使用np.random.default_rng创建的Generator对象，老式的全局随机种子方式在多线程场景下容易出问题。