用 Scipy 定义正交的多元高斯函数答案

【问题标题】：Defining Multivariate Gaussian Function for Quadrature with Scipy用 Scipy 定义正交的多元高斯函数
【发布时间】：2017-02-09 12:06:15
【问题描述】：

我在使用 scipy 为正交定义多元高斯 pdf 时遇到了一些麻烦。我写了一个函数，它接受一个均值向量和协方差矩阵作为输入并返回一个高斯函数。

def make_mvn_pdf(mu, sigma):
    def f(x):
        return sp.stats.multivariate_normal.pdf(x, mu, sigma)
    return f

当我使用 make_mvn_pdf 定义高斯并尝试索引到高斯时，我得到一个没有意义的错误。我首先定义一个均值向量和协方差矩阵并将它们传递给 make_mvn_pdf：

# define covariance matrix
Sigma = np.asarray([[1, .15], [.15, 1]])
# define propagator
B = np.diag([2, 2])
# define data
Obs = np.array([[-0.06895746],[ 0.18778   ]])

# define a Gaussian PDF:
g_int_func = make_mvn_pdf(mean = np.dot(B,Obs[t,:]), cov = Sigma)

我尝试将观察结果传递给密度以获取概率：

testarray=np.random.random((2,2))
g_int_func(testarray)

这会返回以下我不理解的错误。

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-50-083a1915f674> in <module>()
  1 g_int_func = make_mvn_pdf(np.dot(B,Obs[t,:]),Gamma)
----> 2 g_int_func(testarray)

/Users/...in f(x)
 17 def make_mvn_pdf(mu, sigma):
 18     def f(x):
---> 19         return sp.stats.multivariate_normal.pdf(x, mu, sigma)
 20     return f
 21 

/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/scipy/stats/_multivariate.pyc in pdf(self, x, mean, cov, allow_singular)
427 
428         """
--> 429         dim, mean, cov = _process_parameters(None, mean, cov)
430         x = _process_quantiles(x, dim)
431         psd = _PSD(cov, allow_singular=allow_singular)

/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/scipy/stats/_multivariate.pyc in _process_parameters(dim, mean, cov)
 54 
 55     if mean.ndim != 1 or mean.shape[0] != dim:
---> 56         raise ValueError("Array 'mean' must be a vector of length %d." % dim)
 57     if cov.ndim == 0:
 58         cov = cov * np.eye(dim)

ValueError: Array 'mean' must be a vector of length 2.

ValueError 指出数组“均值”必须是长度为 2 的向量，但情况确实如此。实际上，均值和协方差矩阵的维数以及传入的数据的长度都是2。

【问题讨论】：

没关系，但 Gamma 应该是 Sigma 以上，Obs[1,:] 应该只是读取 Obs
你能把你的代码重新排列成MCVE吗？目前，t 未定义，Obs 未完全定义。如果您有一个完整的、独立的示例，我们可以轻松地为您提供帮助，我们可以简单地复制和运行而无需编辑。
谢谢（以上更改），现在可以完美运行了。非常令人沮丧...

标签： python scipy statistics gaussian numerical-integration

【解决方案1】：

您给出的平均值为 np.dot(B, Obs)（考虑到您在评论中建议的更改），其中 B 的形状为 (2, 2)，Obs 的形状为 (2, 1) . dot 调用的结果具有形状 (2, 1)。问题是这是一个二维维数组，multivariate_normal.pdf 期望 mu 是一个一维数组，即形状为 (2,) 的数组. （错误信息使用了“向量”这个词，这是一个糟糕的选择，因为对于很多人来说，形状为 (n, 1) 是的数组是一个（列）向量。它会不那么模棱两可如果错误消息说“数组'均值'必须是长度为2的一维数组。”）

至少有两种简单的方法可以解决这个问题：

您可以确保 Obs 的形状为 (2,) 而不是 (2, 1)，例如Obs = np.array([-0.06895746, 0.18778])。 np.dot(B, Obs) 的形状为 (2,)。
使用ravel 方法“展平”mean 参数：mean=np.dot(B,Obs).ravel()。

【讨论】：