Python如何根据自定义概率密度函数生成符合对应分布的随机数

首先你需要先把已得到的CubicSpline拟合结果修正为合法的概率密度函数（PDF），再通过两种常用方法生成随机数，具体实现如下：

前置操作：归一化PDF #

你原始拟合用的输入是频率值，因此得到的样条函数积分大概率不等于1，需要先做归一化，同时处理区间外的取值：

import numpy as np
from scipy.interpolate import CubicSpline
from scipy.integrate import quad, cumulative_trapezoid

# 假设你的原始数据为随机变量序列X、对应频率序列freq
# 你已有的拟合代码示例
pdf_raw = CubicSpline(X, freq)
x_min, x_max = X.min(), X.max()

# 计算原始样条在有效区间的总积分
total_area, _ = quad(pdf_raw, x_min, x_max)

# 得到归一化后的合法PDF，超出有效区间的概率密度为0
def pdf(x):
    val = pdf_raw(x) / total_area
    # 避免拟合出负的概率密度，做截断处理
    return np.where((x >= x_min) & (x <= x_max) & (val > 0), val, 0)

方法1：逆变换采样（效率更高，推荐） #

利用累积分布函数（CDF）的逆变换实现采样，适合绝大多数场景：

• 第一步：生成CDF查找表
  先在有效区间内生成密集的x网格，计算每个网格点对应的CDF值（即PDF从区间左端点到当前点的积分）：

  # 生成10000个点的网格，密度可根据精度需求调整
  x_grid = np.linspace(x_min, x_max, 10000)
  # 用累积梯形积分快速计算CDF序列
  cdf_grid = cumulative_trapezoid(pdf(x_grid), x_grid, initial=0)
  

• 第二步：拟合逆CDF函数
  因为CDF是严格单调递增函数，可以直接通过插值得到逆CDF（输入为0~1的均匀分布值，输出为对应随机变量取值）：

  from scipy.interpolate import interp1d
  inv_cdf = interp1d(cdf_grid, x_grid, bounds_error=False, fill_value=(x_min, x_max))
  

• 第三步：生成随机数

  # 生成1000个符合目标分布的随机数，修改size参数可调整生成数量
  u = np.random.uniform(0, 1, size=1000)
  random_samples = inv_cdf(u)
  

方法2：拒绝采样（实现更简单） #

不需要计算CDF，实现逻辑简单但采样效率低于逆变换采样，适合小批量采样场景：

• 第一步：先找到PDF在有效区间内的最大值

  pdf_max = pdf(x_grid).max()
  

• 第二步：实现采样逻辑

  def reject_sample(n):
      samples = []
      while len(samples) < n:
          # 随机生成x候选值
          x_cand = np.random.uniform(x_min, x_max)
          # 生成接受阈值
          u = np.random.uniform(0, pdf_max)
          if u <= pdf(x_cand):
              samples.append(x_cand)
      return np.array(samples)
  
  # 生成1000个样本
  random_samples = reject_sample(1000)
  

注意事项 #

• 若你拟合的三次样条在有效区间内出现大量负值，说明原始频率数据的拟合效果不佳，可尝试调整CubicSpline的边界条件，或改用核密度估计拟合PDF。
• 逆变换采样的x网格密度越高，生成样本的精度越高，常规场景下1e4个点足够使用。

内容的提问来源于stack exchange，提问作者ranky123