从自定义经验累积分布函数（CDF）高效抽样：Numpy与Scipy实践指南

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}")

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()