多元伪随机生成:如何从多元分布生成伪随机数据?
Great question! 单变量用逆CDF生成伪随机数据确实是经典操作,但到了多元场景,变量之间的相关性是绕不开的核心问题——毕竟多元分布的关键就是变量间的依赖结构。下面我给你梳理几个最常用的实用方法,覆盖不同的应用场景:
方法1:Copula函数法(分离边际与相关性)
这是处理多元分布最灵活的思路之一,核心是把边际分布和变量间的相关性结构拆分开来处理:
- 步骤1:为每个变量生成独立的均匀分布样本 ( u_i \sim U(0,1) )
- 步骤2:用Copula函数把这些独立的均匀样本转换为具有目标相关性结构的联合均匀样本 ( (u_1, u_2, ..., u_d) )(比如高斯Copula对应正态相关性,t-Copula对应厚尾相关性)
- 步骤3:对每个维度的均匀样本,应用目标边际分布的逆CDF,得到最终的多元样本 ( x_i = F_i^{-1}(u_i) )
举个例子:如果要生成二元“正态边际+高斯Copula”的样本,先生成二元标准正态样本(带指定相关系数),再用单变量正态CDF把每个维度转成均匀分布,最后用目标边际(比如对数正态)的逆CDF得到结果。
方法2:直接联合逆变换(仅限简单可分离场景)
如果多元分布是独立变量的联合,那直接把单变量逆CDF的思路扩展就行:每个维度单独生成符合边际分布的样本,组合起来就是多元样本——毕竟独立变量的联合CDF是各边际CDF的乘积。
但如果变量之间有依赖,这种方法只适用于少数有显式联合逆CDF的分布,比如二维均匀分布在单位正方形,直接生成两个独立均匀变量即可;但大多数复杂依赖的多元分布(比如多元Gamma)没有显式的联合逆,所以这个方法局限性很大。
方法3:接受-拒绝采样(低维场景适用)
和单变量的接受-拒绝思路一致,适合目标分布的PDF容易计算,但不好直接采样的低维多元分布:
- 步骤1:找一个容易采样的提议分布Q(比如多元均匀、多元正态),确保Q的覆盖范围包含目标分布的支撑集,且存在常数M使得 ( f(x) \leq M \cdot Q(x) )(f是目标分布的PDF)
- 步骤2:从Q中采样一个候选样本x,再生成均匀变量 ( u \sim U(0,1) )
- 步骤3:如果 ( u \leq \frac{f(x)}{M \cdot Q(x)} ),就接受x作为目标样本;否则重复步骤2
⚠️ 注意:高维场景下这个方法效率极低(维度诅咒),因为接受概率会随着维度增加急剧下降,所以只适合2-3维的简单分布。
方法4:马尔可夫链蒙特卡洛(MCMC)——高维复杂分布的首选
当面对高维、无显式采样方法的多元分布(比如贝叶斯建模中的后验分布),MCMC是目前最实用的解决方案。核心是构造一个马尔可夫链,让它的平稳分布正好是目标分布,常见的变体有:
- Gibbs采样:如果能方便地从每个变量的条件分布 ( p(x_i | x_{-i}) ) 采样,就轮流固定其他变量,从当前变量的条件分布中采样,重复迭代后得到的样本就近似服从联合分布。比如二元正态分布,条件分布还是正态,非常好采样。
- Metropolis-Hastings(MH)算法:不需要条件分布,只需要计算目标分布的相对概率(比如后验概率的比例),通过提议分布生成候选样本,再根据接受概率决定是否保留。
MCMC的关键是要验证链的收敛性(比如用Gelman-Rubin诊断),确保采样的样本确实来自目标分布。
方法5:特定分布的专用变换法
针对一些常见的多元分布,有专门的高效采样方法,比如:
- 多元正态分布:用Cholesky分解把协方差矩阵Σ拆成 ( L L^T )(L是下三角矩阵),先生成独立的标准正态样本 ( z \sim N(0, I) ),然后通过线性变换 ( x = \mu + L z ) 得到符合 ( N(\mu, \Sigma) ) 的样本。伪代码示例:
import numpy as np def sample_multivariate_normal(mu, sigma, n_samples): d = len(mu) # Cholesky分解协方差矩阵 L = np.linalg.cholesky(sigma) # 生成独立标准正态样本 z = np.random.normal(size=(n_samples, d)) # 线性变换得到目标样本 return mu + z @ L
- Dirichlet分布:可以通过采样d个独立的Gamma变量,再归一化得到;
- 多元均匀分布在凸多面体:可以用拒绝采样或者变换法生成。
总结一下不同方法的适用场景:
- 灵活控制边际和相关性 → Copula法
- 独立变量或简单可分离分布 → 单变量逆CDF组合
- 低维、PDF易计算的分布 → 接受-拒绝采样
- 高维、复杂无显式采样方法的分布 → MCMC
- 常见标准多元分布 → 专用变换法
内容的提问来源于stack exchange,提问作者user25004

