本文详细阐述了如何从自定义的经验累积分布函数(cdf)中生成随机样本。我们将探讨两种主要方法：一是利用numpy的`interp`函数进行基于线性插值的直接抽样，该方法高效且易于实现；二是借助scipy的`interp1d`函数，通过选择不同的插值类型（如线性、三次样条等）实现更平滑的抽样。文章将通过具体的代码示例，指导读者如何在python环境中实现这些抽样技术，并讨论不同方法的适用场景与注意事项。
1. 逆变换抽样原理概述
从任意累积分布函数(CDF)中抽样的核心方法是逆变换抽样（Inverse Transform Sampling）。其基本原理是：如果随机变量$X$的CDF为$F(x)$，且$U$是一个服从$[0, 1]$区间均匀分布的随机变量，那么$Y = F^{-1}(U)$（其中$F^{-1}$是$F$的逆函数）将服从与$X$相同的分布。
对于离散或经验CDF，我们通常没有一个解析的逆函数。在这种情况下，我们可以通过插值来近似$F^{-1}(U)$。具体来说，给定一组$(x_i, F(x_i))$数据点，我们可以将$U$值视为新的$F(x)$值，然后通过插值找到对应的$x$值。
2. 数据准备：定义经验CDF
首先，我们需要定义一个经验累积分布函数。通常，这以一组离散的$(x, \text{cdf})$对的形式给出，其中$x$是某个值，$\text{cdf}$是小于或等于$x$的观测值的比例。
import pandas as pd
import numpy as np
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt
# 定义一个经验CDF
# 'x' 列表示随机变量的值
# 'cdf' 列表示累积概率
cdf_data = pd.DataFrame.from_dict(
{'x':[10e6, 20e6, 50e6, 100e6, 250e6],
'cdf':[0.4, 0.6, 0.7, 0.8, 1]}
)
print("定义的经验CDF数据:")
print(cdf_data)3. 方法一：基于线性插值的CDF抽样 (使用 numpy.interp)
numpy.interp 函数可以用于一维线性插值。在逆变换抽样中，我们将均匀分布的随机数作为新的CDF值（xp），然后根据已知的CDF值（fp）和对应的$x$值（x），来查找插值后的$x$值。
立即学习“Python免费学习笔记（深入）”；
原理： np.interp(x_new, xp, fp) 实际上是在已知点 (xp[i], fp[i]) 的基础上，对 x_new 进行线性插值，以得到对应的 fp 值。在我们的场景中，xp是CDF的概率值，fp是对应的$x$值。我们将均匀分布的随机数作为x_new（即我们希望的CDF概率），然后通过插值找到对应的$x$值。
# 生成10,000个服从[0, 1)均匀分布的随机数
num_samples = 10000
uniform_samples = np.random.uniform(0, 1, num_samples)
# 使用numpy.interp进行线性插值抽样
# uniform_samples 是我们希望的CDF值 (y轴)
# cdf_data['cdf'] 是已知的CDF值 (x轴)
# cdf_data['x'] 是已知的x值 (y轴)
# np.interp 会找到与 uniform_samples 对应的 x 值
samples_method1 = np.interp(uniform_samples, cdf_data['cdf'], cdf_data['x'])
print("\n方法一（线性插值）抽样结果示例 (前5个):")
print(samples_method1[:5])
# 可视化抽样结果的直方图
plt.figure(figsize=(10, 6))
plt.hist(samples_method1, bins=50, density=True, alpha=0.7, color='skyblue', label='线性插值抽样')
plt.title('基于numpy.interp的CDF抽样结果直方图')
plt.xlabel('X 值')
plt.ylabel('概率密度')
plt.grid(True, linestyle='--', alpha=0.6)
plt.legend()
plt.show()4. 方法二：基于更高级插值的CDF抽样 (使用 scipy.interpolate.interp1d)
scipy.interpolate.interp1d 提供了更灵活的插值选项，包括线性、最近邻、零阶、Slerp以及三次样条插值等。这允许我们在抽样时引入不同程度的平滑。
原理： interp1d 首先创建一个插值函数，然后我们可以用这个函数来计算新的点。在逆变换抽样中，我们创建的插值函数将CDF概率映射到$x$值。
# 创建一个插值函数，将CDF概率映射到X值
# x轴是CDF值，y轴是对应的X值
# kind 参数指定插值类型，例如 'linear' (线性), 'cubic' (三次样条)
# fill_value='extrapolate' 允许在定义域之外进行外推，但通常在CDF抽样中应避免
# 或者设置 fill_value=(cdf_data['x'].iloc[0], cdf_data['x'].iloc[-1]) 进行边界钳制
# 线性插值 (与np.interp效果类似)
f_linear = interp1d(cdf_data['cdf'], cdf_data['x'], kind='linear',
bounds_error=False, fill_value=(cdf_data['x'].iloc[0], cdf_data['x'].iloc[-1]))
samples_linear_interp1d = f_linear(uniform_samples)
print("\n方法二（scipy.interpolate.interp1d, 线性插值）抽样结果示例 (前5个):")
print(samples_linear_interp1d[:5])
# 三次样条插值 (更平滑的曲线)
# 注意：三次样条插值需要至少4个数据点
if len(cdf_data) >= 4:
f_cubic = interp1d(cdf_data['cdf'], cdf_data['x'], kind='cubic',
bounds_error=False, fill_value=(cdf_data['x'].iloc[0], cdf_data['x'].iloc[-1]))
samples_cubic_interp1d = f_cubic(uniform_samples)
print("\n方法二（scipy.interpolate.interp1d, 三次样条插值）抽样结果示例 (前5个):")
print(samples_cubic_interp1d[:5])
# 可视化两种插值方法的CDF曲线和抽样结果的直方图
plt.figure(figsize=(14, 6))
# 绘制原始CDF和插值CDF
plt.subplot(1, 2, 1)
plt.plot(cdf_data['x'], cdf_data['cdf'], 'o', label='原始CDF数据点')
x_interp = np.linspace(cdf_data['x'].min(), cdf_data['x'].max(), 500)
cdf_interp_linear = f_linear(np.interp(x_interp, cdf_data['x'], cdf_data['cdf'])) # 逆向插值得到CDF曲线
cdf_interp_cubic = f_cubic(np.interp(x_interp, cdf_data['x'], cdf_data['cdf'])) # 逆向插值得到CDF曲线
# 为了正确绘制逆CDF，我们需要插值CDF值到X值
# 更直接的方式是绘制插值函数本身
cdf_points = np.linspace(0, 1, 100)
plt.plot(f_linear(cdf_points), cdf_points, '-', label='线性插值CDF (逆函数)')
plt.plot(f_cubic(cdf_points), cdf_points, '--', label='三次样条插值CDF (逆函数)')
plt.title('原始与插值CDF的逆函数')
plt.xlabel('X 值')
plt.ylabel('累积概率')
plt.grid(True, linestyle='--', alpha=0.6)
plt.legend()
# 绘制抽样结果直方图
plt.subplot(1, 2, 2)
plt.hist(samples_linear_interp1d, bins=50, density=True, alpha=0.5, color='skyblue', label='线性插值抽样')
plt.hist(samples_cubic_interp1d, bins=50, density=True, alpha=0.5, color='lightcoral', label='三次样条插值抽样')
plt.title('基于interp1d的不同插值抽样结果直方图')
plt.xlabel('X 值')
plt.ylabel('概率密度')
plt.grid(True, linestyle='--', alpha=0.6)
plt.legend()
plt.tight_layout()
plt.show()
else:
print("\n数据点不足，无法进行三次样条插值 (至少需要4个点)。")5. 注意事项与选择建议
-
插值类型选择：
- 线性插值 (numpy.interp 或 interp1d(kind='linear'))： 简单、计算效率高，适用于数据点之间变化相对平稳，或对精度要求不极致的场景。它会产生分段线性的逆CDF，导致抽样结果的概率密度函数（PDF）是分段常数。
- 三次样条插值 (interp1d(kind='cubic'))： 产生更平滑的逆CDF，进而生成更平滑的PDF。适用于需要模拟连续且平滑分布的场景，但计算成本略高，且需要更多的数据点（至少4个）。
- 其他插值类型（如quadratic、nearest等）可根据具体需求选择。
-
边界处理 (fill_value 和 bounds_error)：
- 在interp1d中，bounds_error=False并结合fill_value=(low, high)可以确保当抽样的均匀随机数落在CDF定义域之外时，结果会被钳制在CDF的最小和最大$x$值。这对于避免不合理的抽样结果至关重要。numpy.interp默认会进行边界钳制。
- 对于CDF抽样，均匀随机数理论上总在$[0, 1]$之间，而我们的CDF数据通常覆盖整个$[0, 1]$区间（从0到1）。因此，外推（fill_value='extrapolate'）通常是不必要的，甚至可能导致不合理的样本值。
-
数据点数量与质量：
- 经验CDF的数据点越多，插值结果越接近真实的潜在分布。
- 数据点的分布和精度会直接影响抽样结果的准确性。
-
计算效率：
- numpy.interp通常比scipy.interpolate.interp1d在处理大量样本时更高效，尤其是在只需要线性插值的情况下。
- 如果需要高级插值，interp1d是更合适的选择。
总结
本文介绍了在Python中从自定义经验累积分布函数(CDF)进行抽样的两种主要方法。numpy.interp提供了一种高效的线性插值方法，适用于快速生成基于分段线性逆CDF的样本。而scipy.interpolate.interp1d则提供了更灵活的插值选项，特别是三次样条插值，能够生成更平滑的样本分布，更适合需要模拟连续分布特性的场景。选择哪种方法取决于对抽样结果平滑度和计算效率的具体需求。理解逆变换抽样原理和不同插值方法的特性是有效应用这些技术的关键。
