user2687863 - 1 year ago 135
Python Question

# Defining Multivariate Gaussian Function for Quadrature with Scipy

I'm having some trouble defining a multivariate gaussian pdf for quadrature using scipy. I write a function that takes a mean vector and covariance matrix as input and returns a gaussian function.

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

When I use make_mvn_pdf to define a Gaussian and try to index into the Gaussian I get an error that does not make sense. I begin by defining a mean vector and covariance matrix and passing them into 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)
``````

I try to pass in observations to the density in order to get back probabilities:

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

This returns the following error which I do not understand.

``````---------------------------------------------------------------------------
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.
``````

The ValueError states that the array 'mean' must be a vector of length 2 but this is the case. In fact, the dimension of the mean and covariance matrix and data passed in are all of length 2.

The value that you give as the mean is `np.dot(B, Obs)` (taking into account the change you suggested in a comment), where `B` has shape (2, 2) and `Obs` has shape (2, 1). The result of that `dot` call has shape (2, 1). The problem is that is a two-dimensional array, and `multivariate_normal.pdf` expects `mu` to be a one-dimensional array, i.e. an array with shape (2,). (The error message uses the word "vector", which is a poor choice, because for many people, an array with shape (n, 1) is a (column) vector. It would be less ambiguous if the error message said "Array 'mean' must be a one-dimensional array of length 2.")
• You could ensure that `Obs` has shape (2,) instead of (2, 1), e.g. `Obs = np.array([-0.06895746, 0.18778])`. The `np.dot(B, Obs)` has shape (2,).
• "Flatten" the `mean` argument by using the `ravel` method: `mean=np.dot(B,Obs).ravel()`.