CDF高效抽样教程：Numpy与Scipy实战解析

小伙伴们有没有觉得学习文章很有意思？有意思就对了！今天就给大家带来《自定义CDF高效抽样：Numpy与Scipy实战教程》，以下内容将会涉及到，若是在学习中对其中部分知识点有疑问，或许看了本文就能帮到你！

本文详细介绍了如何从自定义的经验累积分布函数（CDF）中进行数据抽样。我们将利用逆变换抽样原理，结合Python的Numpy和Scipy库，实现两种抽样方法：一是直接基于CDF离散点进行阶梯式插值抽样，二是采用样条等平滑技术对CDF进行插值后抽样，以生成更平滑、更符合实际分布的样本。

理解经验累积分布函数 (Empirical CDF)

经验累积分布函数（Empirical CDF, ECDF）是根据观测数据而非理论模型构建的CDF。它通过统计样本中小于或等于某个特定值的观测比例来估计随机变量的CDF。与理论CDF通常是连续函数不同，经验CDF通常由一系列离散点定义，这些点代表了观测到的值及其对应的累积概率。

例如，一个由DataFrame定义的经验CDF可能包含两列：x 值（随机变量的取值）和 cdf 值（对应 x 的累积概率）。

逆变换抽样原理 (Inverse Transform Sampling)

从任何给定的累积分布函数 F(x) 中进行抽样的核心方法是逆变换抽样（Inverse Transform Sampling）。其基本原理如下：

1. 生成一个在 (0, 1) 区间上均匀分布的随机数 U。
2. 计算 X = F⁻¹(U)，其中 F⁻¹ 是 F 的逆函数。
3. 得到的 X 值将服从原始CDF F(x) 所描述的分布。

在处理经验CDF时，由于 F(x) 通常不是一个解析函数，我们需要通过插值来近似其逆函数 F⁻¹(U)。

方法一：直接从经验CDF抽样（阶梯式插值）

当不需要对CDF进行平滑处理时，可以直接利用经验CDF的离散点进行抽样。这种方法通常采用分段线性插值，将均匀随机数映射到CDF的 x 值。Numpy库的 np.interp 函数非常适合这种场景。

np.interp(x, xp, fp) 函数的作用是：给定一组已知的数据点 (xp, fp)，它会根据 xp 和 fp 之间的关系，对 x 中的每个值进行线性插值，并返回对应的 fp 值。在这里，我们将均匀随机数作为 x，CDF的概率值作为 xp，CDF的 x 值作为 fp，从而实现逆变换抽样。

示例代码：

import pandas as pd
import numpy as np

# 定义一个自定义的经验CDF
cdf_data = pd.DataFrame.from_dict(
    {'x':[10e6, 20e6, 50e6, 100e6, 250e6],
     'cdf':[0.4, 0.6, 0.7, 0.8, 1]}
)

# 1. 生成10,000个在(0, 1)区间上的均匀随机数
num_samples = 10000
uniform_samples = np.random.uniform(0, 1, num_samples)

# 2. 使用numpy.interp进行抽样
# uniform_samples 作为待查找的概率值 (x)
# cdf_data['cdf'] 作为已知概率点 (xp)
# cdf_data['x'] 作为已知x值 (fp)
# np.interp 会根据 uniform_samples 在 cdf_data['cdf'] 中的位置，
# 线性插值出 cdf_data['x'] 中对应的值。
samples_direct = np.interp(uniform_samples, cdf_data['cdf'], cdf_data['x'])

print("直接抽样结果示例 (前10个):")
print(samples_direct[:10])
print(f"抽样结果的最小值: {samples_direct.min():.2f}, 最大值: {samples_direct.max():.2f}")

代码解释：np.interp 函数通过在 cdf_data['cdf'] 中查找 uniform_samples 对应的位置，然后返回 cdf_data['x'] 中相应的线性插值结果。这种方法简单高效，生成的样本值将严格落在原始CDF的 x 范围内。

方法二：平滑经验CDF后抽样（样条插值）

在某些情况下，如果原始CDF数据点较少，或者希望生成的样本分布更平滑、更连续，可以对经验CDF进行平滑插值。Scipy库的 scipy.interpolate.interp1d 函数提供了多种插值方法，包括线性、二次、三次样条等，可以用来构建一个更连续的CDF逆函数。

示例代码：

from scipy.interpolate import interp1d
import matplotlib.pyplot as plt # 可选：用于可视化

# 沿用之前的 cdf_data 和 uniform_samples

# 1. 使用scipy.interpolate.interp1d创建插值函数
# kind参数指定插值类型：'linear'（线性）、'quadratic'（二次）、'cubic'（三次样条）等
# bounds_error=False 允许插值点超出原始数据范围，此时会使用 fill_value 进行外推或填充
# fill_value=(cdf_data['x'].iloc[0], cdf_data['x'].iloc[-1])
#   - 如果插值点小于最小xp值，使用cdf_data['x'].iloc[0]
#   - 如果插值点大于最大xp值，使用cdf_data['x'].iloc[-1]
#   这确保了均匀随机数0和1也能得到有效的x值，且不会超出原始x的物理范围。
cdf_interpolator_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]))
cdf_interpolator_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]))

# 2. 使用插值函数进行抽样
samples_smoothed_linear = cdf_interpolator_linear(uniform_samples)
samples_smoothed_cubic = cdf_interpolator_cubic(uniform_samples)

print("\n平滑抽样结果示例 (线性插值, 前10个):")
print(samples_smoothed_linear[:10])
print(f"抽样结果的最小值: {samples_smoothed_linear.min():.2f}, 最大值: {samples_smoothed_linear.max():.2f}")

print("\n平滑抽样结果示例 (三次样条插值, 前10个):")
print(samples_smoothed_cubic[:10])
print(f"抽样结果的最小值: {samples_smoothed_cubic.min():.2f}, 最大值: {samples_smoothed_cubic.max():.2f}")

# 可选：通过直方图可视化比较不同方法的样本分布
# plt.figure(figsize=(12, 6))
# plt.hist(samples_direct, bins=50, density=True, alpha=0.5, label='直接抽样 (np.interp)')
# plt.hist(samples_smoothed_linear, bins=50, density=True, alpha=0.5, label='平滑抽样 (interp1d, linear)')
# plt.hist(samples_smoothed_cubic, bins=50, density=True, alpha=0.5, label='平滑抽样 (interp1d, cubic)')
# plt.title('不同抽样方法下的样本分布')
# plt.xlabel('X值')
# plt.ylabel('密度')
# plt.legend()
# plt.grid(True)
# plt.show()

代码解释：interp1d 创建了一个可调用的插值函数。通过设置 kind 参数，可以选择不同的插值算法。'linear' 行为与 np.interp 类似，而 'cubic'（三次样条）会生成一条更平滑的曲线，从而产生更连续的样本分布。bounds_error=False 和 fill_value 参数是关键，它们确保了即使均匀随机数落在CDF定义域的边缘（例如，非常接近0或1），也能得到一个有效的 x 值，并将其限制在原始 x 数据的物理范围内。

注意事项与最佳实践

1. 数据准备：确保输入的CDF数据是单调递增的，且 cdf 值范围从0到1。任何违反这些条件的数据都可能导致插值错误或不合理的抽样结果。
2. 插值方法的选择：
   • numpy.interp (线性插值)：简单、高效，适用于对精度要求不高或数据点足够密集的情况。它不会生成超出原始 x 范围的值。
   • scipy.interpolate.interp1d：提供更灵活的插值 kind 选项。
     • kind='linear'：与 np.interp 行为类似，但提供更多控制。
     • kind='cubic' (三次样条)：生成更平滑的曲线，样本分布也更连续，尤其适用于数据点较少但需要平滑过渡的场景。然而，样条插值在某些情况下可能会引入轻微的“过冲”或“欠冲”，导致生成的样本值略微超出原始 x 的最小/最大范围（尽管通过 fill_value 可以有效控制）。
3. fill_value 参数：在使用 interp1d 时，务必设置 bounds_error=False 并合理配置 fill_value。将其设置为CDF的最小和最大 x 值是一个稳健的选择，可以确保所有均匀随机数（包括0和1）都能得到有效的映射，并且样本不会超出合理的物理范围。
4. 样本量：生成足够多的样本（例如10,000或更多）才能更好地反映底层分布的特征。样本量过小可能无法捕捉到分布的细节。
5. 计算效率：对于大规模抽样，numpy.interp 通常比 scipy.interpolate.interp1d 稍快，因为它是一个C实现的简单线性插值。如果对性能有极高要求且仅需线性插值，np.interp 是首选。

总结

本文详细阐述了从自定义经验CDF进行抽样的两种主要方法：直接阶梯式插值和通过平滑函数（如样条）插值。numpy.interp 适用于快速、直接的线性插值抽样，生成的样本严格在原始 x 范围内。scipy.interpolate.interp1d 则提供了更丰富的插值选项，特别是样条插值，能够生成更平滑、更连续的样本分布，适用于对分布平滑性有更高要求的场景。在实际应用中，应根据数据特性、对样本平滑度的要求以及计算效率等因素，选择最合适的抽样策略。

今天关于《CDF高效抽样教程：Numpy与Scipy实战解析》的内容就介绍到这里了，是不是学起来一目了然！想要了解更多关于的内容请关注golang学习网公众号！