Repository: scipy/scipy
Pull request: BUG: multivariate_normal returns a pdf for values outside its support when covariance is singular
Closes #5286.

Before returning a value from `pdf` or `logpdf`, these methods now test that `x` is in the support of the distribution, by removing the projection and seeing if anything is left. The tolerance for this test is the same as that for deciding if the eigenvalues are zero.

The corresponding test, `test_degenerate_distributions` has been modified so that the test points are always in the support. Points guaranteed not to be in the support are also generated and we test that these get a pdf of `0.0`.


Problem statement:
BUG: multivariate_normal returns a pdf for values outside its support when covariance is singular
When a multivariate normal distribution has a singular covariance matrix, its support (i.e. the set of possible values the random variable can take) is restricted to a manifold. In this case, the density with respect to the usual Lebesgue measure is not defined. There IS a density with respect to the Hausdorff measure on the support, and this is exactly what stats.multivariate_normal returns. The problem is that it returns a pdf for any input. Specifically, it returns the density of the projection of the input onto the support. I don't think this behaviour is correct. To demonstrate:

``` python
>>> import numpy as np
>>> from scipy.stats import multivariate_normal
>>> mean = np.zeros(2)
>>> cov = np.ones((2,2))
>>> rv = multivariate_normal(mean=mean, cov=cov, allow_singular=True)
>>> x = np.array([[0.0, 0.0], [1.0, -1.0], [2.0, -2.0]])
>>> x
array([[ 0.,  0.],
       [ 1., -1.],
       [ 2., -2.]])
>>> rv.pdf(x)
array([ 0.28209479,  0.28209479,  0.28209479])
```

For this distribution, the random variable must lie on the line `x[0]-x[1]=0`. `[0,0]` is within the support, so that's all well and good. The other two points are not. Technically, the density is simply undefined at these points. For consistency with other scipy.stats distributions, the output should probably be `0.0`.

Hints from issue discussion:
I've never contributed to scipy before. I've patched this. Do you want it? Do I just submit a pull request? Should I add tests?

https://github.com/scipy/scipy/compare/master...drpeteb:bug-mvn-singular-covariance

This is kind of a tricky issue and I'm not sure how it should have been handled.  I think I've understood what you've said about densities and manifolds and projections, and honestly it's kind of reassuring to hear that the current code at least behaves in an internally consistent (even if poorly chosen) way.

> Technically, the density is simply undefined at these points. 

Just to make sure I understand your proposed patch, the idea is that if a point is farther than some distance `eps` from the subspace then with your patch the pdf will return 0 (instead of the current behavior of using the projection onto the distant subspace), and if the point is closer than `eps` to the subspace then the patch will not change the current behavior?

Yes. That's exactly it. `eps` is the same value as is used to determine which eigenvalues are `0`.

A simpler way to describe the issue is that it currently assigns a probability density to points that could never be generated by the `rvs()` method. A 1D equivalent would be like a gamma distribution giving a density to negative values.

A reasonable behavior when the covariance matrix is singular is to simply raise a `ValueError` (or some other appropriate error).  For example, consider the univariate case, where a singular covariance matrix corresonds to a variance of 0.  This case can be thought of as the limit of the distributions with the variance approaching 0.  As the variance goes to 0, the PDF at 0 goes to infinity.  The limiting case of the PDF would be a delta "function".  I don't see how it could be considered correct to return 1 at 0 when the variance is 0.

@WarrenWeckesser Yes, `allow_singular` is `False` by default, so it raises a LinAlgError on singular matrices. If someone wants to do exotic things involving Hausdorff measures they have to give the keyword argument explicitly.

@drpeteb I think I like the idea.  Technically it would break backwards compatibility; if someone was counting on the function to do the projection implicitly then this patch breaks their code.  Deprecations may be required.

> Do I just submit a pull request? Should I add tests?

Yes, that would be great!

@argriffing Whoops--sorry, I read this issue without going back to look at the code (and without noticing the `allow_singular` argument in the example).  The behavior being discussed for the case `allow_singular=True` seems a bit dicey, numerically speaking, given normal floating point imprecision.   A tiny change in a near-singular covariance matrix that makes the matrix closer to singular can result in PDF values changing from huge to tiny.  I don't really see the use-cases for this.  However, as you point out, it isn't the default behavior, so it won't affect most users.

pinging `multivariate_normal` contributor @jvkersch ...

@WarrenWeckesser I agree - dicey. I don't like the fact that it gets to decide what the rank is for exactly that reason.

@drpeteb  What is the usecase for the degenerate pdf?

I thought the initial proposal was quite reasonable, but I start to see the problems in the interpretation that Warren is pointing out.
For practical purposes I guess that the initial solution is still more likely the expected and useful solution.

I recently saw a similar case where the behavior at the threshold that Warren points out causes weird behavior
http://stackoverflow.com/questions/32599159/robustness-issue-of-statsmodel-linear-regression-ols-python  in that case it is between inv and pinv at the rcond threshold when the reduced rank treatment kicks in.

Also, we have several cases with degenerate multivariate normal distribution (reduced rank cov_params) in statsmodels. I started to worry about the behavior similar to this issue where we would have to look at the appropriate subspace. But because the behavior that we actually use (like the non-degenerate distribution of linear combinations), I didn't need to make up my mind about the desired behavior.

I think the solution is theoretically clean if we define it as the distribution on the subspace, which is the main application in statsmodels of similar degenerate models, e.g. models with linear restrictions, or zero restrictions which have zero variance for some dimensions or subspaces.

The only problem is the numerical issue in the threshold. But for this I think "practicality beats purity".  An almost degenerate distribution has very bad numerical properties and users have to make up their mind whether they want near-singular (eps=0) or singular behavior (eps>0). Same problem is in pinv, linalg.matrix_rank and so on. But the defaults work in 99.99 % (?) of the cases.

Commit message: Apply suggestions from code review

Co-authored-by: Tirth Patel <tirthasheshpatel@gmail.com>

Patch to review:
diff --git a/scipy/stats/_multivariate.py b/scipy/stats/_multivariate.py
index a5777aab2e44..9ffbed44af97 100644
--- a/scipy/stats/_multivariate.py
+++ b/scipy/stats/_multivariate.py
@@ -168,14 +168,26 @@ def __init__(self, M, cond=None, rcond=None, lower=True,
         s_pinv = _pinv_1d(s, eps)
         U = np.multiply(u, np.sqrt(s_pinv))
 
+        # Save the eigenvector basis, and tolerance for testing support
+        self.eps = 1e3*eps
+        self.V = u[:, s <= eps]
+
         # Initialize the eagerly precomputed attributes.
         self.rank = len(d)
         self.U = U
         self.log_pdet = np.sum(np.log(d))
 
-        # Initialize an attribute to be lazily computed.
+        # Initialize attributes to be lazily computed.
         self._pinv = None
 
+    def _support_mask(self, x):
+        """
+        Check whether x lies in the support of the distribution.
+        """
+        residual = np.linalg.norm(x @ self.V, axis=-1)
+        in_support = residual < self.eps
+        return in_support
+
     @property
     def pinv(self):
         if self._pinv is None:
@@ -312,10 +324,16 @@ class multivariate_normal_gen(multi_rv_generic):
                \exp\left( -\frac{1}{2} (x - \mu)^T \Sigma^{-1} (x - \mu) \right),
 
     where :math:`\mu` is the mean, :math:`\Sigma` the covariance matrix,
-    and :math:`k` is the dimension of the space where :math:`x` takes values.
+    :math:`k` the rank of :math:`\Sigma`. In case of singular :math:`\Sigma`,
+    SciPy extends this definition according to [1]_.
 
     .. versionadded:: 0.14.0
 
+    References
+    ----------
+    .. [1] Multivariate Normal Distribution - Degenerate Case, Wikipedia,
+           https://en.wikipedia.org/wiki/Multivariate_normal_distribution#Degenerate_case
+
     Examples
     --------
     >>> import matplotlib.pyplot as plt
@@ -494,6 +512,9 @@ def logpdf(self, x, mean=None, cov=1, allow_singular=False):
         x = self._process_quantiles(x, dim)
         psd = _PSD(cov, allow_singular=allow_singular)
         out = self._logpdf(x, mean, psd.U, psd.log_pdet, psd.rank)
+        if allow_singular and (psd.rank < dim):
+            out_of_bounds = ~psd._support_mask(x-mean)
+            out[out_of_bounds] = -np.inf
         return _squeeze_output(out)
 
     def pdf(self, x, mean=None, cov=1, allow_singular=False):
@@ -519,6 +540,9 @@ def pdf(self, x, mean=None, cov=1, allow_singular=False):
         x = self._process_quantiles(x, dim)
         psd = _PSD(cov, allow_singular=allow_singular)
         out = np.exp(self._logpdf(x, mean, psd.U, psd.log_pdet, psd.rank))
+        if allow_singular and (psd.rank < dim):
+            out_of_bounds = ~psd._support_mask(x-mean)
+            out[out_of_bounds] = 0.0
         return _squeeze_output(out)
 
     def _cdf(self, x, mean, cov, maxpts, abseps, releps):
@@ -728,6 +752,7 @@ def __init__(self, mean=None, cov=1, allow_singular=False, seed=None,
         array([[1.]])
 
         """
+        self.allow_singular = allow_singular
         self._dist = multivariate_normal_gen(seed)
         self.dim, self.mean, self.cov = self._dist._process_parameters(
                                                             None, mean, cov)
@@ -742,6 +767,9 @@ def logpdf(self, x):
         x = self._dist._process_quantiles(x, self.dim)
         out = self._dist._logpdf(x, self.mean, self.cov_info.U,
                                  self.cov_info.log_pdet, self.cov_info.rank)
+        if self.allow_singular and (self.cov_info.rank < self.dim):
+            out_of_bounds = ~self.cov_info._support_mask(x-self.mean)
+            out[out_of_bounds] = -np.inf
         return _squeeze_output(out)
 
     def pdf(self, x):
diff --git a/scipy/stats/tests/test_multivariate.py b/scipy/stats/tests/test_multivariate.py
index b46719e41603..f5f930187e49 100644
--- a/scipy/stats/tests/test_multivariate.py
+++ b/scipy/stats/tests/test_multivariate.py
@@ -35,6 +35,12 @@
 from unittest.mock import patch
 
 
+def _sample_orthonormal_matrix(n):
+    M = np.random.randn(n, n)
+    u, s, v = scipy.linalg.svd(M)
+    return u
+
+
 class TestMultivariateNormal:
     def test_input_shape(self):
         mu = np.arange(3)
@@ -130,22 +136,21 @@ def test_rank(self):
 
     def test_degenerate_distributions(self):
 
-        def _sample_orthonormal_matrix(n):
-            M = np.random.randn(n, n)
-            u, s, v = scipy.linalg.svd(M)
-            return u
-
         for n in range(1, 5):
-            x = np.random.randn(n)
-            for k in range(1, n + 1):
+            z = np.random.randn(n)
+            for k in range(1, n):
                 # Sample a small covariance matrix.
                 s = np.random.randn(k, k)
                 cov_kk = np.dot(s, s.T)
 
-                # Embed the small covariance matrix into a larger low rank matrix.
+                # Embed the small covariance matrix into a larger singular one.
                 cov_nn = np.zeros((n, n))
                 cov_nn[:k, :k] = cov_kk
 
+                # Embed part of the vector in the same way
+                x = np.zeros(n)
+                x[:k] = z[:k]
+
                 # Define a rotation of the larger low rank matrix.
                 u = _sample_orthonormal_matrix(n)
                 cov_rr = np.dot(u, np.dot(cov_nn, u.T))
@@ -172,6 +177,37 @@ def _sample_orthonormal_matrix(n):
                 assert_allclose(logpdf_kk, logpdf_nn)
                 assert_allclose(logpdf_kk, logpdf_rr)
 
+                # Add an orthogonal component and find the density
+                y_orth = y + u[:, -1]
+                pdf_rr_orth = distn_rr.pdf(y_orth)
+                logpdf_rr_orth = distn_rr.logpdf(y_orth)
+
+                # Ensure that this has zero probability
+                assert_equal(pdf_rr_orth, 0.0)
+                assert_equal(logpdf_rr_orth, -np.inf)
+
+    def test_degenerate_array(self):
+        # Test that we can generate arrays of random variate from a degenerate
+        # multivariate normal, and that the pdf for these samples is non-zero
+        # (i.e. samples from the distribution lie on the subspace)
+        k = 10
+        for n in range(2, 6):
+            for r in range(1, n):
+                mn = np.zeros(n)
+                u = _sample_orthonormal_matrix(n)[:, :r]
+                vr = np.dot(u, u.T)
+                X = multivariate_normal.rvs(mean=mn, cov=vr, size=k)
+
+                pdf = multivariate_normal.pdf(X, mean=mn, cov=vr,
+                                              allow_singular=True)
+                assert_equal(pdf.size, k)
+                assert pdf > 0.0
+
+                logpdf = multivariate_normal.logpdf(X, mean=mn, cov=vr,
+                                                    allow_singular=True)
+                assert_equal(logpdf.size, k)
+                assert logpdf > -np.inf
+
     def test_large_pseudo_determinant(self):
         # Check that large pseudo-determinants are handled appropriately.
 


Code review comments:
diff_hunk: @@ -172,6 +177,37 @@ def _sample_orthonormal_matrix(n):
                 assert_allclose(logpdf_kk, logpdf_nn)
                 assert_allclose(logpdf_kk, logpdf_rr)
 
+                # Add an orthogonal component and find the density
+                y_orth = y + u[:, -1]
+                pdf_rr_orth = distn_rr.pdf(y_orth)
+                logpdf_rr_orth = distn_rr.logpdf(y_orth)
+
+                # Ensure that this has zero probability
+                assert_equal(pdf_rr_orth, 0.0)
+                assert_equal(logpdf_rr_orth, -np.inf)
+
+    def test_degenerate_array(self):
+        # Test that we can generate arrays of random variate from a degenerate
+        # multivariate normal, and that the pdf for these samples is non-zero
+        # (i.e. samples from the distribution lie on the subspace)
+        k = 10
+        for n in range(2, 6):
+            for r in range(1, n):
+                mn = np.zeros(n)
+                u = _sample_orthonormal_matrix(n)[:, :r]
+                vr = np.dot(u, u.T)
+                X = multivariate_normal.rvs(mean=mn, cov=vr, size=k)
+
+                pdf = multivariate_normal.pdf(X, mean=mn, cov=vr,
+                                              allow_singular=True)
+                assert_equal(pdf.size, k)
+                assert pdf > 0.0
path: scipy/stats/tests/test_multivariate.py
text: @user1:
```suggestion
                assert np.all(pdf > 0.0)
```

Oops, X is an array of quantiles.

diff_hunk: @@ -172,6 +177,37 @@ def _sample_orthonormal_matrix(n):
                 assert_allclose(logpdf_kk, logpdf_nn)
                 assert_allclose(logpdf_kk, logpdf_rr)
 
+                # Add an orthogonal component and find the density
+                y_orth = y + u[:, -1]
+                pdf_rr_orth = distn_rr.pdf(y_orth)
+                logpdf_rr_orth = distn_rr.logpdf(y_orth)
+
+                # Ensure that this has zero probability
+                assert_equal(pdf_rr_orth, 0.0)
+                assert_equal(logpdf_rr_orth, -np.inf)
+
+    def test_degenerate_array(self):
+        # Test that we can generate arrays of random variate from a degenerate
+        # multivariate normal, and that the pdf for these samples is non-zero
+        # (i.e. samples from the distribution lie on the subspace)
+        k = 10
+        for n in range(2, 6):
+            for r in range(1, n):
+                mn = np.zeros(n)
+                u = _sample_orthonormal_matrix(n)[:, :r]
+                vr = np.dot(u, u.T)
+                X = multivariate_normal.rvs(mean=mn, cov=vr, size=k)
+
+                pdf = multivariate_normal.pdf(X, mean=mn, cov=vr,
+                                              allow_singular=True)
+                assert_equal(pdf.size, k)
+                assert pdf > 0.0
+
+                logpdf = multivariate_normal.logpdf(X, mean=mn, cov=vr,
+                                                    allow_singular=True)
+                assert_equal(logpdf.size, k)
+                assert logpdf > -np.inf
path: scipy/stats/tests/test_multivariate.py
text: @user1:
```suggestion
                assert np.all(logpdf > -np.inf)
```