Defining Multivariate Gaussian Function for Quadrature with Scipy
<p>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.</p>

<pre><code>def make_mvn_pdf(mu, sigma):
    def f(x):
        return sp.stats.multivariate_normal.pdf(x, mu, sigma)
    return f
</code></pre>

<p>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:</p>

<pre><code># 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)
</code></pre>

<p>I try to pass in observations to the density in order to get back probabilities:</p>

<pre><code>testarray=np.random.random((2,2))
g_int_func(testarray)
</code></pre>

<p>This returns the following error which I do not understand.</p>

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

/Users/...in f(x)
 17 def make_mvn_pdf(mu, sigma):
 18     def f(x):
---&gt; 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         """
--&gt; 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:
---&gt; 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.
</code></pre>

<p>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.</p>

<p>The value that you give as the mean is <code>np.dot(B, Obs)</code> (taking into account the change you suggested in a comment), where <code>B</code> has shape (2, 2) and <code>Obs</code> has shape (2, 1).  The result of that <code>dot</code> call has shape (2, 1).  The problem is that is a <em>two</em>-dimensional array, and <code>multivariate_normal.pdf</code> expects <code>mu</code> to be a <em>one</em>-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) <em>is</em> a (column) vector.  It would be less ambiguous if the error message said "Array 'mean' must be a one-dimensional array of length 2.")</p>

<p>There are at least two easy ways to fix the problem:</p>

<ul>
<li>You could ensure that <code>Obs</code> has shape (2,) instead of (2, 1), e.g. <code>Obs = np.array([-0.06895746, 0.18778])</code>.  The <code>np.dot(B, Obs)</code> has shape (2,).</li>
<li>"Flatten" the <code>mean</code> argument by using the <code>ravel</code> method: <code>mean=np.dot(B,Obs).ravel()</code>.</li>
</ul>

python scipy statistics gaussian numerical-integration