在Python中向量化多元正态分布的CDF（累积密度函数）

Question

我如何在Python中矢量化多元正态分布的CDF（累积密度函数）？

当查看this帖子时，我发现有一个Fortran实现的多元CDF已经“移植”到Python中。这意味着我可以轻松地评估一种特定情况下的CDF。

然而，我遇到了很多问题，无法高效地将此函数应用于多个条目。

具体而言，我需要“矢量化”的函数需要4个参数：

• 积分下限（向量）
• 积分上限（向量）
• 正态随机变量的均值（向量）
• 正态随机变量的协方差矩阵（矩阵）

但是，我正在尝试在1000多个元素的列表上高效地评估此函数许多次。

以下是一些代码，以说明我的问题。在下面的示例中，我只是使用随机数据来说明我的观点。

import time
import numpy as np
from scipy.stats.mvn import mvnun # library that calculates MVN CDF
    
np.random.seed(666)
    
iters = 1000 # number of times the whole dataset will be evaluated
obs = 1500 # number of elements in the dataset
dim = 2 # dimension of multivariate normal distribution
     
# Creates a random correlation matrix
def gen_random_corr_mtx(d,k):
    # Source: https://stats.stackexchange.com/a/125020/215330
    #d = number of dimensions
    #k = number of factors. more factors -> smaller correlations

    W = np.random.randn(d*k).reshape((d,k))

    S_temp = np.dot(W,W.T) + np.diag(np.random.rand(d))
    S_norm = np.diag(1/np.sqrt(np.diag(S_temp)))
    S = np.dot(np.dot(S_norm,S_temp),S_norm)
   
    return(S)

# Creates a covariance matrix from a correlation matrix and 
# an array of std devs
def from_cor_to_cov(cor_mtx,st_devs):
    cor_mtx = np.array(cor_mtx)
    st_devs = np.array(st_devs)
    st_devs_diag = np.diag(st_devs)
    cov_mtx = st_devs_diag.dot(cor_mtx).dot(st_devs_diag)
    return(cov_mtx)

# Creating an array with the lower bounds of the integration
lower = np.random.rand(obs,dim)

# Creating an array with the upper bounds of the integration
upper = lower + np.random.rand(obs,dim)

# Creating an array with the means of the distributions
means = np.random.rand(obs,dim)

# Generating the random covariance matrices
covs = []
for i in range(obs):
    # Making sure the covariance matrix is positive semi-definite
    while True:
        cor_mtx = gen_random_corr_mtx(dim,2)
        st_devs = np.abs(np.random.randn(dim))
        cov_mtx = from_cor_to_cov(cor_mtx,st_devs)
        if np.all(np.linalg.eigvals(cov_mtx) > 0):
            break
    covs.append(cov_mtx)
covs = np.array(covs)
    
# Here is where the trouble starts.
time_start = time.time()
for i in range(iters):
    results = []
    for j in range(obs):
        this_p, this_i = mvnun(lower[j],upper[j],means[j],covs[j])
        results.append(this_p)
time_end = time.time()
    
print(time_end-time_start)
# > 4.090040922164917

我有一个包含1500个观测值的数据集，需要进行1000次评估。在我的机器上，这需要4.090040922164917秒来计算。

请注意，我并不试图摆脱外部for循环（i循环）。我只是为了模拟我的真实问题而创建了它。我真正想要消除的是内部循环（j循环）。

如果我找到一种有效地评估整个数据集的CDF的方法，执行时间可以大大缩短。

我知道mvnun函数最初是用Fortran编写的（原始代码here），并使用f2pye“移植”到Python中，如here所示。

有人能帮我吗？我已经开始研究theano，但似乎我唯一的选择是使用scan函数，这可能也不是很好的改进。

谢谢！！！

更新（2023年4月14日）

scipy库最近进行了更新，现在用户可以确定multivariate_normal.cdf()方法中积分下限，而不总是默认的n维数组-inf值。

以下是使用新方法更新后的代码部分：

import scipy 

time_start = time.time()
for i in range(iters):
    results = []
    for j in range(obs):
        this_p = scipy.stats.multivariate_normal.cdf(lower_limit=lower[j],
                                                     x=upper[j],
                                                     mean=means[j],
                                                     cov=covs[j])
        results.append(this_p)
time_end = time.time()

print(time_end-time_start)
# > 268.0528531074524

上述方法比我最初提出的方法慢得多，需要超过4分钟才能计算完成。但由于大多数scipy的类和方法已经设置为矢量化操作，我猜测/希望我在这里发布的代码只需要轻微调整即可真正实现矢量化。

- Felipe D.

---

我知道这是一个老问题，但在生成随机协方差矩阵时要小心。为了有效，您的协方差矩阵需要是半正定的（如果您想使分布非退化，则需要是正定的）。Scikit-learn提供了一个函数来实现这一点：make_spd_matrix。现在，mvnun似乎并不关心第四个参数covs[j]是否作为协方差矩阵有意义，这真的很奇怪。 - josephmure

谢谢提醒！我添加了一些更改以确保矩阵是半正定的 =) - Felipe D.

Answer 1

这只是部分答案，但是如果多元正态分布的维数较小（2或3）且协方差矩阵保持不变，则有一种增加速度的方法。

import numpy as np
import openturns as ot

def computeRectangularDomainProbability(lower, upper, means, cov_matrix):
    """
    Compute the probability of a rectangular solid
    under a multinormal distribution.
    
    """
    # Center the bounds of the rectangular solid on the mean
    lower -= means
    upper -= means
    
    # The same covariance matrix for all rectangular solids.
    cov_matrix = ot.CovarianceMatrix(cov_matrix)
    
    # This way, we only need to define one multivariate normal distribution.
    # That is the trick that allows vectorization.
    dimension = lower.shape[1]   
    multinormal = ot.Normal([0.0] * dimension, cov_matrix)
    
    # The probability of the rectangular solid is a weighted sum
    # of the CDF of the vertices (with weights equal to 1 or -1).
    # The following block computes the CDFs and applies the correct weight.
    full_reverse_binary = np.array(list(bin(2**dimension)[:1:-1]), dtype=int)
    prob = 0.0
    for i in range(2**dimension):
        reverse_binary = np.array(list(bin(i)[:1:-1]), dtype=int)
        reverse_binary = np.append(reverse_binary, 
                                   np.zeros(len(full_reverse_binary) - 
                                            len(reverse_binary) -
                                            1)).astype(int)
        point = np.zeros(lower.shape)
        for num, digit in enumerate(reverse_binary):
            if digit:
                point[:, num] = upper[:, num]
            else:
                point[:, num] = lower[:, num]
        cdf = np.array(multinormal.computeCDF(point))
        if (reverse_binary.sum() % 2) == (dimension % 2):
            prob += cdf
        else:
            prob -= cdf
    
    return prob.reshape(-1,)

测试脚本：维度2

iters = 1000 # loop size
obs = 1500 # number of rectangular solids
dim = 2 # dimension of multivariate normal distribution

import time
import numpy as np
from scipy.stats.mvn import mvnun # library that calculates MVN CDF
from sklearn.datasets import make_spd_matrix
import openturns as ot

time_mvnun = 0.0
time_openturns = 0.0
discrepancy = 0.0
np.random.seed(0)

for iteration in range(iters):

    lower = np.random.rand(obs,dim)
    upper = lower + np.random.rand(obs,dim)
    means = np.random.rand(obs,dim)
    
    # Generating the random covariance matrices with sklearn
    # to make sure they are positive semi-definite        
    cov_mtx = make_spd_matrix(dim)
        
    time_start = time.time()
    results = []
    for j in range(obs):
        this_p, this_i = mvnun(lower[j],upper[j],means[j],cov_mtx)
        results.append(this_p)
    results = np.array(results)
    time_end = time.time()
    time_mvnun += time_end - time_start
    
    
    time_start = time.time()       
    otparallel = computeRectangularDomainProbability(lower, upper, means, cov_mtx)
    time_end = time.time()
    time_openturns += time_end - time_start
    
    mvnorm_vs_otparallel = np.abs(results - otparallel).sum()
    discrepancy += mvnorm_vs_otparallel

print('Dimension {}'.format(dim))

# Print computation time
print('mvnun     time: {0:e}'.format(time_mvnun))
print('openturns time: {0:e}'.format(time_openturns))
print('ratio mvnun/ot: {0:f}'.format(time_mvnun / time_openturns))

# Check that the results are the same for mvnum and openturns
print('mvnun-openturns result discrepancy: {0:e}'.format(discrepancy))

我的机器上的输出结果：

Dimension 2
mvnun     time: 4.040635e+00
openturns time: 3.588211e+00
ratio mvnun/ot: 1.126086
mvnun-openturns result discrepancy: 8.057912e-11

有轻微的加速：略高于10%。

三维

让我们更改控制脚本的全局变量。

iters = 100 # loop size
obs = 1500 # number of rectangular solids
dim = 3 # dimension of multivariate normal distribution

我的电脑上的输出：

Dimension 3
mvnun     time: 2.378337e+01
openturns time: 1.596872e+00
ratio mvnun/ot: 14.893725
mvnun-openturns result discrepancy: 4.537064e-03

在三维空间中，提出的代码速度提高了15倍，效果显著。

四维空间

不幸的是，当维度为4时，openturns变得非常缓慢。它包含了对于1、2和3维空间的CDF的智能实现，但对于大于3维的空间，则退回到一个更慢、更通用的实现。

iters = 1 # loop size
obs = 15 # number of rectangular solids
dim = 4 # dimension of multivariate normal distribution

Dimension 4
mvnun     time: 7.289171e-03
openturns time: 3.689714e+01
ratio mvnun/ot: 0.000198
mvnun-openturns result discrepancy: 6.297527e-07

在四维空间中，该代码的速度慢了约四个数量级！这可能是因为在四维空间中，每个矩形体需要计算 16=2^4 个累积分布函数（CDF），而每次计算的速度比较慢。

- josephmure