已知X、Y时如何模拟三元高斯分布Z值？求Python简便方案

不需要手动实现Gibbs采样或复杂的条件参数计算，我们可以直接基于多元正态分布的条件分布特性，结合Scipy的工具快速完成采样，甚至可以封装成通用函数，比自定义类更灵活。

核心逻辑 #

对于联合正态分布 $(\mathbf{X}, \mathbf{Y}) \sim \mathcal{N}(\boldsymbol{\mu}, \boldsymbol{\Sigma})$，给定已知变量 $\mathbf{X}=\mathbf{x}$ 时，未知变量 $\mathbf{Y}$ 的条件分布满足：
$$\mathbf{Y} \mid \mathbf{X} = \mathbf{x} \sim \mathcal{N}\left( \boldsymbol{\mu}Y + \boldsymbol{\Sigma}{YX}\boldsymbol{\Sigma}{XX}^{-1}(\mathbf{x} - \boldsymbol{\mu}X), \boldsymbol{\Sigma}{YY} - \boldsymbol{\Sigma}{YX}\boldsymbol{\Sigma}{XX}^{-1}\boldsymbol{\Sigma}{XY} \right)$$
其中 $\boldsymbol{\mu}$ 和 $\boldsymbol{\Sigma}$ 是联合分布的均值向量与协方差矩阵，按已知/未知变量拆分成分块矩阵即可计算条件分布的参数。

简便实现代码 #

import numpy as np
from scipy.stats import multivariate_normal

def sample_conditional_normal(given_indices, given_values, mean, cov, size=1):
    """
    从条件多元正态分布中采样
    参数:
        given_indices: 已知变量的索引列表（如X、Y对应索引[0,1]）
        given_values: 已知变量的取值数组（如[10,20]）
        mean: 联合分布的均值向量
        cov: 联合分布的协方差矩阵
        size: 采样数量
    返回:
        采样结果数组，形状为(size, 未知变量数量)
    """
    # 划分已知/未知变量的索引
    all_indices = np.arange(len(mean))
    unknown_indices = np.setdiff1d(all_indices, given_indices)
    
    # 拆分均值向量与协方差矩阵
    mu_x = mean[given_indices]
    mu_y = mean[unknown_indices]
    sigma_xx = cov[np.ix_(given_indices, given_indices)]
    sigma_yx = cov[np.ix_(unknown_indices, given_indices)]
    sigma_yy = cov[np.ix_(unknown_indices, unknown_indices)]
    
    # 计算条件分布的均值与协方差
    sigma_xx_inv = np.linalg.inv(sigma_xx)
    cond_mean = mu_y + sigma_yx @ sigma_xx_inv @ (given_values - mu_x)
    cond_cov = sigma_yy - sigma_yx @ sigma_xx_inv @ cov[np.ix_(given_indices, unknown_indices)]
    
    # 生成采样结果
    return multivariate_normal(mean=cond_mean, cov=cond_cov).rvs(size=size)

# 示例使用
mean4 = np.array([2, 3, 4, 5])
cov_matrix4 = np.array([[1, 0.5, 0.3, 0.2],
                       [0.5, 1, 0.4, 0.1],
                       [0.3, 0.4, 1, 0.15],
                       [0.2, 0.1, 0.15, 1]])

# 给定X(索引0)=10、Y(索引1)=20，模拟Z(索引2)
sim_z = sample_conditional_normal(given_indices=[0,1], given_values=[10,20], mean=mean4, cov=cov_matrix4, size=1000)

# 给定X(0)=10、Y(1)=20、Z(2)=5，模拟W(3)
sim_w = sample_conditional_normal(given_indices=[0,1,2], given_values=[10,20,5], mean=mean4, cov=cov_matrix4, size=1000)

更省心的现成包 #

如果连分块矩阵都不想手动处理，可以用概率编程库自动完成条件分布的参数计算：

• Pyro：Facebook开源的概率编程库，支持直接定义条件分布

  import pyro
  import pyro.distributions as dist
  
  def pyro_conditional_sample(given_indices, given_values, mean, cov, size=1):
      joint_dist = dist.MultivariateNormal(loc=mean, covariance_matrix=cov)
      # 固定已知变量的取值，生成条件分布
      cond_dist = joint_dist.mask(False).condition(dict(zip(given_indices, given_values)))
      return cond_dist.sample((size,)).numpy()
  

• PyMC3：另一个主流概率编程库，同样支持一键生成条件分布并采样

方案优势 #

• 通用性强：支持任意数量的已知/未知变量，无需修改核心逻辑
• 效率更高：依赖Numpy矩阵运算，比Python循环实现快数倍
• 可维护性好：逻辑清晰，避免手动计算的出错风险

内容的提问来源于stack exchange，提问作者Víťa Horák