Python中从自定义经验累积分布函数(CDF)抽样：直接与平滑插值方法

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)

# 生成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()

# 创建一个插值函数，将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个点)。")