CODEWITHSUNDEEP
pythonscipystatisticsgaussiannumerical-integration
user2687863
user2687863

Reputation: 13

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.

Upvotes: 1

Views: 4056

Answers (1)

Warren Weckesser
Warren Weckesser

Reputation: 114781

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.")

There are at least two easy ways to fix the problem:

  • 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().

Upvotes: 5

Related Questions