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

[Robert Kern]

On Sun, Aug 4, 2013 at 11:54 AM, Joris Vankerschaver <
joris.vankerschaver at gmail.com> wrote:
>
> 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?

A new scipy.stats.multivariate module would be a fine place for it. This is
common enough that it doesn't need to wait for a full-fledged framework for
multivariate distributions, in my opinion. I'm not sure what such a
framework would really do except collect the various PDFs and log-PDFs.

--
Robert Kern
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20130804/8944be00/attachment.html>