Repository: sympy/sympy
Pull request: stats - fixes issue#21057 by change value check to positive semi-definite
<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->


#### Brief description of what is fixed or changed
Fixes #21057 and now joint normal distribution with a positive semi-definite covariance matrix

#### Other comments


#### Release Notes

<!-- Write the release notes for this release below between the BEGIN and END
statements. The basic format is a bulleted list with the name of the subpackage
and the release note for this PR. For example:

* solvers
  * Added a new solver for logarithmic equations.

* functions
  * Fixed a bug with log of integers.

or if no release note(s) should be included use:

NO ENTRY

See https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more
information on how to write release notes. The bot will check your release
notes automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
* stats
    * Added support for positive semi-definite covariance matrix in joint Multivariate normal distribution.
<!-- END RELEASE NOTES -->


Problem statement:
Cannot create joint normal distribution with positive semidefinite covariance matrix
I want to create a joint normal distributed random variable with a covariance matrix that is positive semidefinite but not positive definite. This seems to be impossible:

```python
import sympy.stats

sympy.stats.Normal("x", [0, 0], [[0,0], [0,0]])
```
gives

```
  File "<ipython-input-28-3bd23a2d609a>", line 1, in <module>
    sympy.stats.Normal("x", [0, 0], [[0,0], [0,0]])

  File "/home/.local/lib/python3.8/site-packages/sympy/stats/crv_types.py", line 2966, in Normal
    return MultivariateNormal(name, mean, std)

  File "/home/.local/lib/python3.8/site-packages/sympy/stats/joint_rv_types.py", line 209, in MultivariateNormal
    return multivariate_rv(MultivariateNormalDistribution, name, mu, sigma)

  File "/home/.local/lib/python3.8/site-packages/sympy/stats/joint_rv_types.py", line 28, in multivariate_rv
    dist.check(*args)

  File "/home/.local/lib/python3.8/site-packages/sympy/stats/joint_rv_types.py", line 143, in check
    _value_check(sigma.is_positive_definite,

  File "/home/.local/lib/python3.8/site-packages/sympy/stats/rv.py", line 1487, in _value_check
    raise ValueError(message)

ValueError: The covariance matrix must be positive definite. 
```

I don't know if there is some implementation reason for this, but mathematically any positive semidefinite matrix should be accepted. I checked numpy and there it works to do:

```python
import numpy as np

np.random.multivariate_normal([0, 0], [[0, 0], [0, 0]])
```

with the expected output:

```
Out[30]: array([0., 0.])
```

Hints from issue discussion:
It seems like the fix is simply changing `sigma.is_positive_definite` to `sigma.is_positive_semidefinite`  in https://github.com/sympy/sympy/blob/2346054bb4888ef7eec2f6dad6c3dd52bf1fe927/sympy/stats/joint_rv_types.py#L138-L144

Reminds me of [this](https://github.com/sympy/sympy/pull/19935) PR. I haven't checked the code for the multivariate distributions but I suspect `pdf`, `cdf` and a whole lot of other functions would have to be changed as well.

It might be simpler just to work with constants (since I think their `pdf` and `cdf` are defined somewhere).
@Maelstrom6  Hey , I am bit new to open source contribution and would like to work on this issue, Could you give me an headstart on how to start working on this issue?
I don't have any experience with the random vectors part of SymPy but there would probably be at least 4 things that need to change:
- The assertion to check if the covariance Matrix is positive should be changed to positive semidefinite.
- The pdf should also change. According to Wikipedia:
  > In the degenerate case where the covariance matrix is singular, the corresponding distribution has no density; see the section below for details. This case arises frequently in statistics; for example, in the distribution of the vector of residuals in the ordinary least squares regression.
  
  This means that we should probably have a piecewise or if statement in the pdf function. Considering that a piecewise expression is not what people are expecting and it slows things down a bit, an if statement to return some kind of SymPy object representing a Nan would be good.
- cdf should similarly be changed if the code defines the cdf.
- The sampling which would be the main reason why people use a positive semidefinite matrix should be changed. I'm not sure how the code looks so that could either require no change or it could be more complicated. 
- Testing should be done on every applicable function for random vectors. 

There are probably other things to figure out along the way. 
A normal distribution with zero variance is also not permitted:
```python
In [1]: from sympy.stats import *

In [2]: X = Normal('X', 0, 0)
---------------------------------------------------------------------------
ValueError: Standard deviation must be positive 
```
> A normal distribution with zero variance is also not permitted

The scalar normal distribution did have a PR for zero variance but it was closed because we couldn't reach an agreement on what the PDF and CDF should be and because the Wikipedia page says that the variance must be strictly positive. 

The user could also replace `X = Normal('X', mu, sigma)` with `X = mu` to get the same functionality (except for sampling capabilities). 
The case of zero variance for a single variate normal distribution is analogous to the case of positive semi-definite but not positive definite for the joint normal distribution. In my opinion sympy should allow it. Some arguments as for why it should:

1. Numpy allows it

```
np.random.normal(1, 0)
Out[1]: 1.0
```

2. It arises naturally from numerical situations. For example if you have two datasets like for example

```
a = [1, 1, 1, 1, 1]
b = [1, 1, 1, 1, 2]
```

The sample variance would be 0 for ```a``` and 0.16 for ```b```. If we want to model both as normally distributed (this assumption might be wrong but it might also be correct) in sympy we only do it directly for b and not for a, which would lead to an unnecessary special case.

3. Mathematically the 0 variance case or perhaps the limit when the variance approaching zero is usually considered to be the dirac distribution.

Perhaps sympy could support it by letting `Normal('X', 1 , 0)` return a dirac distribution.
Given that SymPy represents random variables as expressions wouldn't a dirac distribution just be `mu`?
Not really. Let's say we want to calculate the probability of the random variable falling in some range: `P(a < X < b)` This would then be the integral from `a` to `b` of the pdf of the distribution. If the range includes `mu` the integral should become 1 else 0. A constant expression really cannot capture this. I think the dirac distribution needs to implemented as a distribution.
It seems to work as far as I can tell:
```python
In [19]: mu = Symbol('mu')

In [20]: density(mu)
Out[20]: x ↦ δ(-μ + x)

In [21]: E(mu)
Out[21]: μ

In [22]: integrate(density(mu)(z), (z, 0, 1))
Out[22]: -θ(-μ) + θ(1 - μ)

In [23]: integrate(density(mu)(z), (z, mu-1, mu+1))
Out[23]: 1

In [24]: integrate(density(mu)(z), (z, mu-1, mu-0.5))
Out[24]: 0
```
Apart from these:
```python
In [25]: variance(mu)   # <--- incorrect
Out[25]: μ

In [26]: E(mu**2) - E(mu)**2
Out[26]: 0

In [27]: P(mu < 1)
---------------------------------------------------------------------------
ValueError: Expression containing Random Variable expected, not mu < 1

In [28]: P(mu < mu+1)
---------------------------------------------------------------------------
ValueError: Expression containing Random Variable expected, not mu < mu + 1
```

> I think the dirac distribution needs to implemented as a distribution.

I'm not sure what this means. SymPy's stats module does not really provide for working with distributions rather than expressions involving random variables e.g. https://github.com/sympy/sympy/issues/20643#issuecomment-749473694 https://github.com/sympy/sympy/pull/20560#issuecomment-740792093

Commit message: changed docstrings

Patch to review:
diff --git a/sympy/stats/joint_rv_types.py b/sympy/stats/joint_rv_types.py
index 51ccb7089ba5..57b2ab20bd84 100644
--- a/sympy/stats/joint_rv_types.py
+++ b/sympy/stats/joint_rv_types.py
@@ -139,10 +139,10 @@ def set(self):
     def check(mu, sigma):
         _value_check(mu.shape[0] == sigma.shape[0],
             "Size of the mean vector and covariance matrix are incorrect.")
-        #check if covariance matrix is positive definite or not.
+        #check if covariance matrix is positive semi definite or not.
         if not isinstance(sigma, MatrixSymbol):
-            _value_check(sigma.is_positive_definite,
-            "The covariance matrix must be positive definite. ")
+            _value_check(sigma.is_positive_semidefinite,
+            "The covariance matrix must be positive semi definite. ")
 
     def pdf(self, *args):
         mu, sigma = self.mu, self.sigma
@@ -182,6 +182,7 @@ def MultivariateNormal(name, mu, sigma):
     mu : List representing the mean or the mean vector
     sigma : Positive definite square matrix
         Represents covariance Matrix
+        If `sigma` is noninvertible then only sampling is supported currently
 
     Returns
     =======
diff --git a/sympy/stats/tests/test_joint_rv.py b/sympy/stats/tests/test_joint_rv.py
index c05161c5bad6..2f9339c08a59 100644
--- a/sympy/stats/tests/test_joint_rv.py
+++ b/sympy/stats/tests/test_joint_rv.py
@@ -25,8 +25,6 @@ def test_Normal():
     assert density(m)(1, 2) == 1/(2*pi)
     assert m.pspace.distribution.set == ProductSet(S.Reals, S.Reals)
     raises (ValueError, lambda:m[2])
-    raises (ValueError,\
-        lambda: Normal('M',[1, 2], [[0, 0], [0, 1]]))
     n = Normal('B', [1, 2, 3], [[1, 0, 0], [0, 1, 0], [0, 0, 1]])
     p = Normal('C',  Matrix([1, 2]), Matrix([[1, 0], [0, 1]]))
     assert density(m)(x, y) == density(p)(x, y)
@@ -354,3 +352,19 @@ def test_sample_seed():
                 assert all(s1 != s2)
         except NotImplementedError:
             continue
+
+
+def test_issue_21057():
+    m = Normal("x", [0, 0], [[0, 0], [0, 0]])
+    n = MultivariateNormal("x", [0, 0], [[0, 0], [0, 0]])
+    assert m == n
+    libraries = ['scipy', 'numpy', 'pymc3']
+    for library in libraries:
+        try:
+            imported_lib = import_module(library)
+            if imported_lib:
+                s1 = sample(m, size=8)
+                s2 = sample(n, size=8)
+                assert tuple(s1.flatten()) == tuple(s2.flatten())
+        except NotImplementedError:
+            continue


Code review comments:
diff_hunk: @@ -182,6 +182,7 @@ def MultivariateNormal(name, mu, sigma):
     mu : List representing the mean or the mean vector
     sigma : Positive definite square matrix
path: sympy/stats/joint_rv_types.py
text: @user1:
Should it say `sigma : Positive semidefinite square matrix` instead? 

@author:
Well , yes i guess since all Positive definite square matrix are subset of Positive semidefinite square matrix's , I'll update the change.