[SciPy-Dev] Multivariate normal distribution in scipy.stats

Joris Vankerschaver joris.vankerschaver@gmail....
Sun Aug 4 05:54:45 CDT 2013

Dear all, 

I was wondering if there's any interest to have the multivariate normal distribution integrated into scipy.stats. Its PDF is easily calculated, for instance by something like 

def multinormal_pdf(r, mean, cov):
    """ Probability density function for a multidimensional Gaussian distribution."""
    dim  = r.shape[-1]
    dev  = r - mean
    maha = np.einsum('...k,...kl,...l->...', dev, np.linalg.pinv(cov), dev)
    return (2 * np.pi)**(-0.5 * dim) * np.linalg.det(cov)**(0.5) * np.exp(-0.5 * maha)

Here, r is an array-like of position vectors, with the last axis labeling the components. 

While all of this would be useful and not too hard to code up, it doesn't seem to fit within the scipy.stats.rv_continuous framework, which seems to be explicitly geared to 1D random variables. Is there a natural place where this code could be fitted in, or is this functionality already somewhere else in SciPy?

With best wishes,

More information about the SciPy-Dev mailing list