Repository: sympy/sympy

Issue:
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.])
```

Discussion:
sympy-bot: :white_check_mark:

Hi, I am the [SymPy bot](https://github.com/sympy/sympy-bot) (v161). I'm here to help you write a release notes entry. Please read the [guide on how to write release notes](https://github.com/sympy/sympy/wiki/Writing-Release-Notes).



* No release notes entry will be added for this pull request.

<details><summary>Click here to see the pull request description that was parsed.</summary>

    …i definite matrix

    <!-- 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 -->
    NO ENTRY
    <!-- END RELEASE NOTES -->


</details><p>

Maelstrom6: Sorry about the late response, studying is getting hectic.

Some other tests that we can add to make sure that it's usable:
- Check that E(X) is a vector of zeros
- Check that Var(X) is a matrix of zeros
- It might be a good idea to test that the sample produces a vector of zeros too
- This might be a big problem but preferably `E(X*Y.T)` should be able to work too. But I think that might be beyond the scope of the stats module at the moment

Are there any other things that people might use a random variable with a semidefinite matrix for?
praneethratna: @Maelstrom6 i'm kinda struck up here , E(X) returns an error stating sympy.matrices.common.NonInvertibleMatrixError: Matrix det == 0; not invertible. and similar error for variance(X) , could you help me figure this out?
Maelstrom6: So I found out today that SymPy computes expectations by first principles every time for every kind of distribution. Nothing is precomputed. This module is long overdue for an overhaul and, when it was first designed, it didn't have multivariate distributions in mind. We might just have to leave out expectation calculations for now.

The sampling seems to work fine so it might be helpful to add a test for that. You can look at similar tests done by other distributions in the testing module. 

I'm worried that with general use-cases other than sampling, the user will get lots of `NonInvertibleMatrixError`s.
praneethratna: > So I found out today that SymPy computes expectations by first principles every time for every kind of distribution. Nothing is precomputed. This module is long overdue for an overhaul and, when it was first designed, it didn't have multivariate distributions in mind. We might just have to leave out expectation calculations for now.
> 
> The sampling seems to work fine so it might be helpful to add a test for that. You can look at similar tests done by other distributions in the testing module.
> 
> I'm worried that with general use-cases other than sampling, the user will get lots of `NonInvertibleMatrixError`s.

But does sympy support  a pdf of dirac delta function ? Then should'nt it be feasible to create a Normal distribution with 0 Standard deviation and expectation?
github-actions[bot]: Benchmark results from GitHub Actions

Lower numbers are good, higher numbers are bad. A ratio less than 1
means a speed up and greater than 1 means a slowdown. Green lines
beginning with `+` are slowdowns (the PR is slower then master or
master is slower than the previous release). Red lines beginning
with `-` are speedups.

Significantly changed benchmark results (PR vs master)
```diff

```
Significantly changed benchmark results (master vs previous release)
```diff
       before           after         ratio
     [ed9a550f]       [f57d98d8]
     <sympy-1.8^0>                 
-      1.20±0.01s          162±3ms     0.13  dsolve.TimeDsolve01.time_dsolve
-      10.3±0.03s       4.97±0.02s     0.48  integrate.TimeIntegrationRisch02.time_doit(100)
-      10.4±0.02s       5.02±0.07s     0.48  integrate.TimeIntegrationRisch02.time_doit_risch(100)
-      18.2±0.6ms       12.0±0.2ms     0.66  logic.LogicSuite.time_load_file
-        89.0±3μs       34.1±0.9μs     0.38  matrices.TimeDiagonalEigenvals.time_eigenvals
-         101±3μs         64.6±2μs     0.64  matrices.TimeMatrixGetItem.time_ImmutableDenseMatrix_getitem
-       104±0.8μs         68.1±2μs     0.66  matrices.TimeMatrixGetItem.time_ImmutableSparseMatrix_getitem
-        109±10μs         68.5±2μs     0.63  matrices.TimeMatrixGetItem.time_MutableSparseMatrix_getitem
-        95.0±4μs         61.4±2μs     0.65  solve.TimeMatrixArithmetic.time_dense_add(10, 0)
+      13.9±0.7μs       21.1±0.8μs     1.52  solve.TimeMatrixArithmetic.time_dense_add(3, 0)
+        16.6±1μs         36.4±1μs     2.19  solve.TimeMatrixArithmetic.time_dense_add(3, 5)
+      26.1±0.3μs         47.9±2μs     1.84  solve.TimeMatrixArithmetic.time_dense_add(4, 5)
+        44.9±2μs         73.3±4μs     1.63  solve.TimeMatrixArithmetic.time_dense_add(6, 5)
-     1.38±0.03ms          252±3μs     0.18  solve.TimeMatrixArithmetic.time_dense_multiply(10, 0)
-        58.3±1μs       28.3±0.9μs     0.49  solve.TimeMatrixArithmetic.time_dense_multiply(3, 0)
+        90.9±5μs          149±5μs     1.64  solve.TimeMatrixArithmetic.time_dense_multiply(3, 5)
-         119±4μs         35.3±2μs     0.30  solve.TimeMatrixArithmetic.time_dense_multiply(4, 0)
-         340±7μs         72.7±2μs     0.21  solve.TimeMatrixArithmetic.time_dense_multiply(6, 0)
-         108±4μs         48.4±1μs     0.45  solve.TimeMatrixOperations.time_det(3, 0)
-         210±6μs          122±2μs     0.58  solve.TimeMatrixOperations.time_det(3, 2)
-         205±5μs          111±2μs     0.54  solve.TimeMatrixOperations.time_det(3, 5)
-         111±2μs         50.3±3μs     0.45  solve.TimeMatrixOperations.time_det_bareiss(3, 0)
-         196±8μs          120±5μs     0.61  solve.TimeMatrixOperations.time_det_bareiss(3, 2)
-         196±7μs          110±3μs     0.56  solve.TimeMatrixOperations.time_det_bareiss(3, 5)
-         116±2μs         49.8±2μs     0.43  solve.TimeMatrixOperations.time_det_berkowitz(3, 0)
-         203±4μs          122±2μs     0.60  solve.TimeMatrixOperations.time_det_berkowitz(3, 2)
-         197±4μs          109±5μs     0.55  solve.TimeMatrixOperations.time_det_berkowitz(3, 5)
+         865±9μs      1.39±0.04ms     1.61  solve.TimeMatrixOperations.time_det_berkowitz(4, 0)
+     1.05±0.06ms      2.05±0.06ms     1.96  solve.TimeMatrixOperations.time_det_berkowitz(4, 2)
+     1.06±0.05ms      2.09±0.05ms     1.97  solve.TimeMatrixOperations.time_det_berkowitz(4, 5)
+         327±5μs         547±10μs     1.68  solve.TimeMatrixOperations.time_rank(3, 0)
+        523±20μs         812±20μs     1.55  solve.TimeMatrixOperations.time_rank(4, 0)
+         142±3μs         231±10μs     1.63  solve.TimeMatrixOperations.time_rref(3, 0)
+         286±5μs         444±20μs     1.55  solve.TimeMatrixOperations.time_rref(4, 0)
-      5.87±0.2ms       3.22±0.2ms     0.55  solve.TimeRationalSystem.time_linsolve(10)
-     1.28±0.07ms         773±20μs     0.60  solve.TimeRationalSystem.time_linsolve(5)
+      7.37±0.2ms       11.3±0.7ms     1.53  solve.TimeSparseSystem.time_linear_eq_to_matrix(30)
-     1.47±0.06ms         909±30μs     0.62  solve.TimeSparseSystem.time_linsolve_eqs(10)
-      2.78±0.2ms      1.63±0.06ms     0.59  solve.TimeSparseSystem.time_linsolve_eqs(20)
-      3.99±0.2ms      2.42±0.04ms     0.61  solve.TimeSparseSystem.time_linsolve_eqs(30)

```
Full benchmark results can be found as artifacts in GitHub Actions
(click on checks at the top of the PR).
Maelstrom6: SymPy does have a Dirac delta function and we could return it as a PDF but then again, Wikipedia says that if the matrix is noninvertible, then the PDF does not exist. I'm not sure how correct they are there.

I'd guess a last thing that you can change for now is to change the docstring of `MultivariateNormal`. Sigma can now be positive semidefinite but only supports sampling.
praneethratna: > SymPy does have a Dirac delta function and we could return it as a PDF but then again, Wikipedia says that if the matrix is noninvertible, then the PDF does not exist. I'm not sure how correct they are there.
> 
> Edit:
> With `Normal('X', [0, 0], [0, 0, 0, 1])`, I don't think it is possible to express the PDF even with Dirac Deltas. It might be useful to add a test for some normal distribution of that nature for sampling as well just to make sure it works.
> 
> I'd guess a last thing that you can change for now is to change the docstring of `MultivariateNormal`. Sigma can now be positive semidefinite but only supports sampling.

I've modified the docstring , is it fine or should i add that statement in parameters part?
Maelstrom6: Maybe combine both statements so that we get something like:

```
sigma : Positive semidefinite square matrix
        Represents covariance Matrix
        If `sigma` is noninvertible then only sampling is supported currently
```

Once that's done, we could get @oscarbenjamin to see if they have any other suggestions.
oscarbenjamin: Maybe if the covariance matrix is a zero matrix then a constant expression can be returned instead of a random variable.

Patch:
diff --git a/sympy/stats/joint_rv_types.py b/sympy/stats/joint_rv_types.py
index 51ccb7089b..26e1d84a5d 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
@@ -180,8 +180,9 @@ def MultivariateNormal(name, mu, sigma):
     ==========
 
     mu : List representing the mean or the mean vector
-    sigma : Positive definite square matrix
+    sigma : Positive semidefinite square matrix
         Represents covariance Matrix
+        If `sigma` is noninvertible then only sampling is supported currently
 
     Returns
     =======


Test patch:
diff --git a/sympy/stats/tests/test_joint_rv.py b/sympy/stats/tests/test_joint_rv.py
index c05161c5ba..72c7086a43 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,23 @@ 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]])
+    p = Normal("x", [0, 0], [[0, 0], [0, 1]])
+    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)
+                s3 = sample(p, size=8)
+                assert tuple(s1.flatten()) == tuple(s2.flatten())
+                for s in s3:
+                    assert tuple(s.flatten()) in p.pspace.distribution.set
+        except NotImplementedError:
+            continue