python - 在 Python 中矢量化多元正态 CDF(累积密度函数)

Question

如何在 Python 中矢量化多元正态 CDF(累积密度函数)？
当看 this发布后，我发现有一个“移植”到 Python 的多变量 CDF 的 Fortran 实现。这意味着我可以针对一个特定案例轻松评估 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

lower = np.random.rand(obs,dim)
upper = lower + np.random.rand(obs,dim)
means = np.random.rand(obs,dim)

# Creates a symmetric matrix - used for the random covariance matrices
def make_sym_matrix(dim,vals):
    m = np.zeros([dim,dim])
    xs,ys = np.triu_indices(dim,k=1)
    m[xs,ys] = vals[:-dim]
    m[ys,xs] = vals[:-dim]
    m[ np.diag_indices(dim) ] = vals[-dim:]
    return m

# Generating the random covariance matrices
covs = []
for i in range(obs):
    cov_vals = np.random.rand(int((dim**2 + dim)/2))
    cov_mtx = make_sym_matrix(dim,cov_vals)
    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)

在这里，我有一个包含 1500 个观察值的数据集，我对其进行了 1000 次评估。在我的机器上，这需要 6.74399995804 秒来计算。
请注意，我并不是要摆脱外部 for 循环(在 i 上)。我只是创建它来模仿我的真正问题。我的 for 循环 真的试图消除的是内部的(超过 j)。
执行时间可能是 大大如果我找到了一种有效评估整个数据集的 CDF 的方法，则减少。
我知道 mvnun 函数最初是用 Fortran 编写的(原始代码 here )并使用 f2pye“移植”到 Python，如 here 所示.
谁能帮我解决这个问题？我已经开始研究 theano，但似乎我唯一的选择是使用扫描功能，这也可能没有太大的改进。
谢谢!!!

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%。
维度 3
让我们更改控制脚本的全局变量。

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

在维度 3 中的增益更为显着:建议的代码快了 15 倍。
维度 4
不幸的是，openturns 在维度 4 上慢了很多。它包含了维度 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

在第 4 维中，建议的代码慢了大约 4 个数量级!这可能是因为在维度 4 中，它需要为每个长方体计算 16=2^4 个 CDF，并且这些计算中的每一个都比在较小维度中慢。