【发布时间】：2018-03-05 22:03:04
【问题描述】：
如何在 Python 中矢量化多元正态 CDF（累积密度函数）？
查看this 帖子时，我发现有一个多变量CDF 的Fortran 实现已“移植”到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
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，但似乎我唯一的选择是使用 scan 功能，这也可能没有太大的改进。
谢谢！！！
【问题讨论】：
-
我知道这是一个老问题，但在生成随机协方差矩阵时要小心。为了有效，您的协方差矩阵需要是半正定的（如果您希望分布是非退化的，甚至是正定的）。 Scikit-learn 提供了一个函数来做到这一点：make_spd_matrix。现在，
mvnun似乎并不关心这一点，无论第四个参数covs[j]作为协方差矩阵是否有意义，它都会返回一个结果，这真的很奇怪。
标签： python statistics vectorization normal-distribution cdf