Repository: scipy/scipy
Pull request: ENH: add `fit` method to `multivariate_normal`
#### Reference issue
Closes #15676

#### What does this implement/fix?
Implements a method to fit a multivariate normal distribution to data. Still in very premature stage to get some feedback.

#### Additional information
I am not sure how to generate test cases for the correctness of the fit. Ideally we would compare that the result is actually the MLE. For that we would need to compare against a reference found by numerical optimization (if someone knows a reference implementation in Matlab or R, please let me know) . Optimizing over positive definite matrices is unfortunately difficult, so not sure how to best proceed here.

Problem statement:
ENH: Implement `multivariate_normal.fit`
### Is your feature request related to a problem? Please describe.

`multivariate_normal.fit` is not available, but is trivial to implement.

### Describe the solution you'd like.

Should be available.

### Describe alternatives you've considered.

Doing it yourself. It is inconvenient.

### Additional context (e.g. screenshots, GIFs)

```python
def fitGaussian(points):
    m = np.mean(points, axis=0)
    ms = points - m
    c = ms.T @ ms / (len(points) - 1)
    return scipy.stats.multivariate_normal(mean=m, cov=c)
```

Hints from issue discussion:
This sounds reasonable. Multivariate distributions do not have a lot of shared functionality, so I think it's OK if this is only defined for `multivariate_normal` initially. 

Not that the univariate `fit` interface is perfect, but it is what we have, so I think the interface of a multivariate `fit` method should be similar. 

- One thing is that it should return a tuple of the resulting mean and covariance rather than a frozen distribution.
- Also, if it were to match univariate `fit`, it would accept `fmean` and `fcov` to fix one parameter or the other. In that case, would `c` as calculated above still be the optimal covariance matrix if I were to fix the mean with ` floc`? (and would `m` still the optimal mean if I were to specify `fcov`?)
>One thing is that it should return a tuple of the resulting mean and covariance rather than a frozen distribution.

Mean is available as a property of the distribution object, so is covariance matrix.

>In that case, would c as calculated above still be the optimal covariance matrix if I were to fix the mean with floc? (and would m still the optimal mean if I were to specify fcov?)

You have raised the right issue. I'm not even sure it makes any sense to fix the parameters to any value other than true distribution parameter (garbage in - garbage out). So if we assumme that the parameters provided are true, then they are the same we compute ourselves, then there should be no problem to skip their computation and just use the provided ones.

For mean I'm pretty sure that it would be optimal because covariance only squeezes and rotates the distribution, and because mean is a prereq to MLE covariance.

For covariance I'm not so sure.

https://scikit-learn.org/stable/modules/covariance.html should be related.


> Mean is available as a property of the distribution object, so is covariance matrix.

Certainly. Returning a frozen distribution would be a valid thing to do, but I'm suggesting we follow the precedent set by the `fit` method of univariate distributions, which is to just return the fit parameters.

> I'm not even sure it makes any sense to fix the parameters to any value other than true distribution parameter (garbage in - garbage out).

The true distribution parameters are not necessarily the same as the maximum likelihood estimates. 
For this and other reasons, if we allow the user to fix the mean or covariance, we cannot assume that they are fixed at the maximum likelihood estimates.
The formulas we use for the maximum likelihood estimates need to consider this. Sometimes the MLE of one parameter doesn't depend on the value of another parameter; sometimes it does. Sometimes it does, but you wouldn't know it from looking at the usual formulas because, e.g., the formula for the MLE of the covariance simplifies when the MLE of the mean is assumed. I just think we need to check on this before we give the user the option of fixing one parameter or the other. 

If it's hard to determine (and we don't want to do the math ourselves), maybe initially we just _don't_ allow the user to fix either (don't accept `fmean` or `fcov`).  
Thanks for the clarification.

>maybe initially we just don't allow the user to fix either (don't accept fmean or fcov).

Good idea.


i wrote the method the but the problem boils down to in ability to access cov of a multivariate_normal_frozen during test, do you know how could it be done?

the error given to when i tested the method
AttributeError: 'multivariate_normal_frozen' object has no attribute 'cov'
Poking around shows something called `cov_object`:

```
In [19]: x = np.linspace(0, 5, 10, endpoint=False)
    ...: frozen = multivariate_normal(mean=2.5, cov=0.5)

In [20]: dir(frozen)
Out[20]: 
['__class__',
 '__delattr__',
 '__dict__',
 '__dir__',
 '__doc__',
 '__eq__',
 '__format__',
 '__ge__',
 '__getattribute__',
 '__gt__',
 '__hash__',
 '__init__',
 '__init_subclass__',
 '__le__',
 '__lt__',
 '__module__',
 '__ne__',
 '__new__',
 '__reduce__',
 '__reduce_ex__',
 '__repr__',
 '__setattr__',
 '__sizeof__',
 '__str__',
 '__subclasshook__',
 '__weakref__',
 '_dist',
 'abseps',
 'allow_singular',
 'cdf',
 'cov_object',
 'dim',
 'entropy',
 'logcdf',
 'logpdf',
 'maxpts',
 'mean',
 'pdf',
 'random_state',
 'releps',
 'rvs']

In [21]: frozen.cov_object.covariance
Out[21]: array([[0.5]])
```
I don't see any problem with fixing either parameter for multivariate normal
mean does not depend on cov
cov for given mean just needs to demean with the provided fixed mean. 

If you want to add it also for multivariate t, then I have somewhere a recursive estimator that iterates between mean, cov updates for fixed df, IIRC.

> I don't see any problem with fixing either parameter for multivariate normal mean does not depend on cov cov for given mean just needs to demean with the provided fixed mean.

That could be a follow-up of #18361 indeed.

> 
> If you want to add it also for multivariate t, then I have somewhere a recursive estimator that iterates between mean, cov updates for fixed df, IIRC.

That sounds interesting. Could you provide a reference or even better, that code?

Commit message: TST: properly escape error message

Patch to review:
diff --git a/scipy/stats/_multivariate.py b/scipy/stats/_multivariate.py
index 0fdda070f915..1f0d1f1ad14e 100644
--- a/scipy/stats/_multivariate.py
+++ b/scipy/stats/_multivariate.py
@@ -786,6 +786,39 @@ def entropy(self, mean=None, cov=1):
         dim, mean, cov_object = self._process_parameters(mean, cov)
         return 0.5 * (cov_object.rank * (_LOG_2PI + 1) + cov_object.log_pdet)
 
+    def fit(self, x):
+        """Fit the multivariate normal to data `x`.
+
+        Parameters
+        ----------
+        x : ndarray (m, n)
+            Data the distribution is fitted to. Must have two axes.
+            The first axis of length `m` represents the number of vectors
+            the distribution is fitted to. The second axis of length `n`
+            determines the dimensionality of the fitted distribution.
+
+        Returns
+        -------
+        mean : ndarray (n, )
+            Estimated mean vector
+        cov : ndarray (n, n)
+            Estimated covariance matrix
+
+        """
+        # input validation
+        x = np.asarray(x)
+        if x.ndim != 2:
+            raise ValueError("`x` must be two-dimensional of shape (m, n).")
+
+        n_vectors, dim = x.shape
+
+        # parameter estimation
+        # reference: https://home.ttic.edu/~shubhendu/Slides/Estimation.pdf
+        mean = x.mean(axis=0)
+        centered_data = x - mean
+        cov = centered_data.T @ centered_data / n_vectors
+        return mean, cov
+
 
 multivariate_normal = multivariate_normal_gen()
 
diff --git a/scipy/stats/tests/test_multivariate.py b/scipy/stats/tests/test_multivariate.py
index ec16e94c7775..a2f268301a62 100644
--- a/scipy/stats/tests/test_multivariate.py
+++ b/scipy/stats/tests/test_multivariate.py
@@ -12,6 +12,8 @@
 
 from .test_continuous_basic import check_distribution_rvs
 
+import re
+
 import numpy
 import numpy as np
 
@@ -803,6 +805,18 @@ def test_mean_cov(self):
         ref = multivariate_normal.pdf(x, [1, 1, 1], cov_object)
         assert_equal(multivariate_normal.pdf(x, 1, cov=cov_object), ref)
 
+    def test_fit_error(self):
+        data = [1, 3]
+        error_msg = "`x` must be two-dimensional of shape (m, n)."
+        with pytest.raises(ValueError, match=re.escape(error_msg)):
+            multivariate_normal.fit(data)
+
+    def test_fit_correctness(self):
+        x = np.array([[-1., -1., -1.], [0., 0., 0.], [1., 1., 1.]])
+        mean_est, cov_est = multivariate_normal.fit(x)
+        assert_allclose(mean_est, np.zeros((3, )), atol=1e-15)
+        assert_allclose(cov_est, np.full((3, 3), 2/3), rtol=1e-15)
+
 
 class TestMatrixNormal:
 


Code review comments:
diff_hunk: @@ -803,6 +805,18 @@ def test_mean_cov(self):
         ref = multivariate_normal.pdf(x, [1, 1, 1], cov_object)
         assert_equal(multivariate_normal.pdf(x, 1, cov=cov_object), ref)
 
+    def test_fit_error(self):
+        data = [1, 3]
+        error_msg = "`x` must be two-dimensional of shape (m, n)."
+        with pytest.raises(ValueError, match=re.escape(error_msg)):
path: scipy/stats/tests/test_multivariate.py
text: @user1:
```suggestion
        with pytest.raises(ValueError, match=error_msg):
```

diff_hunk: @@ -786,6 +786,39 @@ def entropy(self, mean=None, cov=1):
         dim, mean, cov_object = self._process_parameters(mean, cov)
         return 0.5 * (cov_object.rank * (_LOG_2PI + 1) + cov_object.log_pdet)
 
+    def fit(self, x):
+        """Fit the multivariate normal to data `x`.
+
+        Parameters
+        ----------
+        x : ndarray (m, n)
+            Data the distribution is fitted to. Must have two axes.
+            The first axis of length `m` represents the number of vectors
+            the distribution is fitted to. The second axis of length `n`
+            determines the dimensionality of the fitted distribution.
+
+        Returns
+        -------
+        mean : ndarray (n, )
+            Estimated mean vector
+        cov : ndarray (n, n)
+            Estimated covariance matrix
+
+        """
+        # input validation
+        x = np.asarray(x)
+        if x.ndim != 2:
+            raise ValueError("`x` must be two-dimensional of shape (m, n).")
path: scipy/stats/_multivariate.py
text: @user1:
`m` and `n` are not defined in the error message (only in the documentation). Rather than defining them, let's just remove them, since they don't add any information.
```suggestion
            raise ValueError("`x` must be two-dimensional.")
```
Incidentally, I think this will also obviate the need to do anything special when matching the message. (e.g. no use of `re`, although there are likely other ways of doing that.)

diff_hunk: @@ -803,6 +805,18 @@ def test_mean_cov(self):
         ref = multivariate_normal.pdf(x, [1, 1, 1], cov_object)
         assert_equal(multivariate_normal.pdf(x, 1, cov=cov_object), ref)
 
+    def test_fit_error(self):
+        data = [1, 3]
+        error_msg = "`x` must be two-dimensional of shape (m, n)."
path: scipy/stats/tests/test_multivariate.py
text: @user1:
```suggestion
        error_msg = "`x` must be two-dimensional."
```

diff_hunk: @@ -786,6 +786,39 @@ def entropy(self, mean=None, cov=1):
         dim, mean, cov_object = self._process_parameters(mean, cov)
         return 0.5 * (cov_object.rank * (_LOG_2PI + 1) + cov_object.log_pdet)
 
+    def fit(self, x):
+        """Fit the multivariate normal to data `x`.
+
+        Parameters
+        ----------
+        x : ndarray (m, n)
+            Data the distribution is fitted to. Must have two axes.
+            The first axis of length `m` represents the number of vectors
+            the distribution is fitted to. The second axis of length `n`
+            determines the dimensionality of the fitted distribution.
+
+        Returns
+        -------
+        mean : ndarray (n, )
+            Estimated mean vector
+        cov : ndarray (n, n)
+            Estimated covariance matrix
path: scipy/stats/_multivariate.py
text: @user1:
We should specify that this is the MLE somewhere.
```suggestion
        mean : ndarray (n, )
            Maximum likelihood estimate of the mean vector
        cov : ndarray (n, n)
            Maximum likelihood estimate of the covariance matrix
```

diff_hunk: @@ -12,6 +12,8 @@
 
 from .test_continuous_basic import check_distribution_rvs
 
+import re
+
path: scipy/stats/tests/test_multivariate.py
text: @user1:
```suggestion
```

diff_hunk: @@ -786,6 +786,39 @@ def entropy(self, mean=None, cov=1):
         dim, mean, cov_object = self._process_parameters(mean, cov)
         return 0.5 * (cov_object.rank * (_LOG_2PI + 1) + cov_object.log_pdet)
 
+    def fit(self, x):
+        """Fit the multivariate normal to data `x`.
path: scipy/stats/_multivariate.py
text: @user1:
```suggestion
        """Fit a multivariate normal distribution to data.
```
The main point is:
<img width="931" alt="image" src="https://github.com/scipy/scipy/assets/6570539/cb8c29b9-bb06-4c8e-a325-55f284249452">