给定X和Y时如何从高斯三变量分布进行模拟？

Question

有没有简单的Python包解决以下问题？

我有 3 个变量，假设它们共享三变量正态分布 - 所以我们知道均值和协方差矩阵。当给定 X、Y 时，有没有简单的方法来模拟 Z 值？

我已经通过使用吉布斯采样器编码了双变量模拟，它可以扩展到从多变量分布中采样，可惜它有点令人困惑和烦人。 真的，没有更简单的方法从条件分布中采样吗？

谢谢你。

这是我尝试过的，我相信它有效。但我宁愿直接使用一些包，而不需要自己计算：

import numpy as np
import scipy

class GibbsMultivariateSampler():
    def __init__(
            self, 
            data: np.ndarray,
            means: np.ndarray, 
            covariance: np.ndarray):

            self.means=means
            self.covariance=covariance
            self.data = data
                        
    def conditioned_mean(self,x: np.ndarray):
        i=0
        conditioned_means=[]
        for element in x:           
            i=i+1
            _h=len(x)
            _l = _h - i
            conditioned_mean=(self.covariance[_h, _l] / self.covariance[_l, _l]) * (element - self.means[_l])
            conditioned_means.append(conditioned_mean)
        return self.means[_h]+sum(conditioned_means)

    def conditioned_covariance(self,x: np.ndarray):
        i=0
        conditioned_covariances=[]
        for element in x:
            i=i+1
            _h=len(x)
            _l = _h - i
            conditioned_covariance=-1*((self.covariance[_h, _l] ** 2  / self.covariance[_l, _l]))
            conditioned_covariances.append(conditioned_covariance)
        return self.covariance[_h, _h]+ sum(conditioned_covariances)


    def simulate(self,size):
        self.conditioned_covariance(self.data)
        conditioned_distribution = scipy.stats.multivariate_normal(mean=self.conditioned_mean(self.data), cov=self.conditioned_covariance(self.data))
        return conditioned_distribution.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]])                


#simulate Z given X,Y
sampler=GibbsMultivariateSampler(data=np.array([10,20]), means=mean4, covariance=cov_matrix4)
simulation=sampler.simulate(1000)

#simulate W given X,Y,Z
sampler=GibbsMultivariateSampler(data=np.array([10,20,5]), means=mean4, covariance=cov_matrix4)
simulation=sampler.simulate(1000)

Answer 1

我相信您想知道如何生成具有正确均值和方差/协方差结构的高斯随机变量。以下是如何使用 Cholesky 分解执行此操作的 python/numpy 实现。

import numpy as np
import random

M = np.array([2, 3, 4, 5])
V = 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]])

print("Show the covariance matrix\n")
print(V)               

print("\nShow Cholesky factorization of covariance matrix\n")                
L = np.linalg.cholesky(V)
print(L)

print("\nConfirm that L L.T = V\n")
print(np.dot(L, L.T))

print("\nShow input vector of standard normals\n")
Z = np.random.default_rng().normal(size = 4)
print(Z)

print("\nSample of correlated Gaussian results\n")
print(np.dot(L, Z) + M)

---

如果您确实想知道如何在给定 X₁,...,X_k-1 的情况下生成 X_k，可以反转此方法以确定相应的标准法线 Z₁,.. .,Z_k-1，生成Z_k，并将扩展的Z向量乘以L的第k行。