Fel*_* D. 5 python statistics normal-distribution vectorization cdf
如何在Python中向量化多元标准CDF（累积密度函数）？
在查看这篇文章时，我发现有一个“移植”到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((dim^2 - dim)/2+dim)
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语言编写（原代码在这里），并使用f2pye可以看出“移植”到Python 在这里。
谁能帮我这个忙吗？我已经开始研究theano了，但是看来我唯一的选择就是使用扫描功能，这可能也没有太大的改进。
谢谢！！！
这只是部分答案，但是如果多元正态分布的维数较小（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,)
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
在维度 3 中，增益更为显着：建议的代码速度快了 15 倍。
不幸的是，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，并且这些计算中的每一个都比在较小维度中慢。