Repository: sympy/sympy
Pull request: [GSoC] Added `is_random` to check random expressions of stats
<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->
Added `is_random` function to check if expression contains random variables.
#### 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. -->
Fixes #12283
References from #19299 

#### Brief description of what is fixed or changed
* Added a function `is_random` to check if the expression contains random variables in it.
For example:
```python
>>> from sympy.stats.rv import is_random
>>> from sympy.stats import Normal, BernoulliProcess
>>> from sympy.abc import x
>>> X = Normal('X',1,2)
>>> is_random(X)
True
>>> is_random(X>3)
True
>>> is_random(x>3)
False
>>> is_random(x)
False
>>> Y = BernoulliProcess('Y',0.3)
>>> is_random(Y[2])
True
```
#### Other comments


#### Release Notes

<!-- Write the release notes for this release below. 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 `is_random` in `sympy.stats.rv` to check if expression contains random variables
<!-- END RELEASE NOTES -->
ping @czgdp1807 @Upabjojr @jmig5776 

Problem statement:
SymPy stats expansion
Add

- `is_random` to every SymPy object (returns true if expression contains random symbols).
- test `is_random` on undefined functions and indexed bases.
- stats functions such as `P`, `E`, `var`, and `cov` should return their respective objects in case the `RandomSymbol` has no specified probability distribution.

Hints from issue discussion:
May be, `is_random` property can be added to `Basic` class(or whatever is the top most parent class) which can compute it's value by using the `OR` over all the `is_random` values of the objects args. Default value of `Symbol.is_random` should be `False` and `RandomSymbol.is_random` should be `True` and since every expression contains symbols, so for an expression `is_random` can be computed in O(number of args) complexity.
@czgdp1807 since Python 3.5 methods can be added externally to the class through singledispatch.

We have dropped support for Python 2.7 and 3.4 this year, so we should start using this way in order not to overload too much `sympy/core`.
I see, that's also possible. However, any thoughts on the approach? 
ping @Smit-create 
One way is to add a method in class Expr
```python
    @property
    def is_random(self):
        from sympy.stats.rv import random_symbols
        if random_symbols(self):
            return True
        return False
```
I am trying to get an exact approach and will come up with that soon
Both Basic and Expr have far too many methods, attributes and properties. Why can't this just be a normal function?
>  Why can't this just be a normal function?

We presently have a function `random_symbols` to get the random symbols from an expression.
So why is a method needed?
https://github.com/sympy/sympy/pull/19299/files#diff-91b28bab28724762f84e77df3facdb54R15
I will update it in my PR #19290
I would recommend to add `is_random` in a separate PR, it's easier to review.

Commit message: return objects when no pspace

Patch to review:
diff --git a/sympy/stats/crv.py b/sympy/stats/crv.py
index 61a8fb7a1935..c38f6966ef3f 100644
--- a/sympy/stats/crv.py
+++ b/sympy/stats/crv.py
@@ -21,7 +21,7 @@
 from sympy.solvers.inequalities import reduce_rational_inequalities
 from sympy.core.sympify import _sympify
 from sympy.external import import_module
-from sympy.stats.rv import (RandomDomain, SingleDomain, ConditionalDomain,
+from sympy.stats.rv import (RandomDomain, SingleDomain, ConditionalDomain, is_random,
         ProductDomain, PSpace, SinglePSpace, random_symbols, NamedArgsMixin)
 
 scipy = import_module('scipy')
@@ -555,7 +555,7 @@ def probability(self, condition, **kwargs):
         except NotImplementedError:
             from sympy.stats.rv import density
             expr = condition.lhs - condition.rhs
-            if not random_symbols(expr):
+            if not is_random(expr):
                 dens = self.density
                 comp = condition.rhs
             else:
diff --git a/sympy/stats/crv_types.py b/sympy/stats/crv_types.py
index 576ad128b952..7120fc5bcf95 100644
--- a/sympy/stats/crv_types.py
+++ b/sympy/stats/crv_types.py
@@ -65,7 +65,7 @@
 from sympy.stats.crv import SingleContinuousPSpace, SingleContinuousDistribution
 from sympy.stats.joint_rv import JointPSpace, CompoundDistribution
 from sympy.stats.joint_rv_types import multivariate_rv
-from sympy.stats.rv import _value_check, RandomSymbol
+from sympy.stats.rv import _value_check, is_random
 
 oo = S.Infinity
 
@@ -127,7 +127,7 @@ def rv(symbol, cls, args):
     dist = cls(*args)
     dist.check(*args)
     pspace = SingleContinuousPSpace(symbol, dist)
-    if any(isinstance(arg, RandomSymbol) for arg in args):
+    if any(is_random(arg) for arg in args):
         pspace = JointPSpace(symbol, CompoundDistribution(dist))
     return pspace.value
 
diff --git a/sympy/stats/drv_types.py b/sympy/stats/drv_types.py
index 14914df119b1..399336b015a3 100644
--- a/sympy/stats/drv_types.py
+++ b/sympy/stats/drv_types.py
@@ -21,7 +21,7 @@
 from sympy.stats import density
 from sympy.stats.drv import SingleDiscreteDistribution, SingleDiscretePSpace
 from sympy.stats.joint_rv import JointPSpace, CompoundDistribution
-from sympy.stats.rv import _value_check, RandomSymbol
+from sympy.stats.rv import _value_check, is_random
 from sympy.external import import_module
 
 numpy = import_module('numpy')
@@ -45,7 +45,7 @@ def rv(symbol, cls, *args):
     dist = cls(*args)
     dist.check(*args)
     pspace = SingleDiscretePSpace(symbol, dist)
-    if any(isinstance(arg, RandomSymbol) for arg in args):
+    if any(is_random(arg) for arg in args):
         pspace = JointPSpace(symbol, CompoundDistribution(dist))
     return pspace.value
 
diff --git a/sympy/stats/error_prop.py b/sympy/stats/error_prop.py
index 44c6862e52e9..ecc24dec1adf 100644
--- a/sympy/stats/error_prop.py
+++ b/sympy/stats/error_prop.py
@@ -5,6 +5,7 @@
 
 from sympy import S, Symbol, Add, Mul, simplify, Pow, exp
 from sympy.stats.symbolic_probability import RandomSymbol, Variance, Covariance
+from sympy.stats.rv import is_random
 
 _arg0_or_var = lambda var: var.args[0] if len(var.args) > 0 else var
 
@@ -59,7 +60,7 @@ def variance_prop(expr, consts=(), include_covar=False):
     if len(args) == 0:
         if expr in consts:
             return S.Zero
-        elif isinstance(expr, RandomSymbol):
+        elif is_random(expr):
             return Variance(expr).doit()
         elif isinstance(expr, Symbol):
             return Variance(RandomSymbol(expr)).doit()
diff --git a/sympy/stats/frv_types.py b/sympy/stats/frv_types.py
index 5f4a3185f920..10c8db3c1d25 100644
--- a/sympy/stats/frv_types.py
+++ b/sympy/stats/frv_types.py
@@ -26,7 +26,7 @@
 from sympy.tensor.array import ArrayComprehensionMap
 from sympy.stats.frv import (SingleFiniteDistribution,
                              SingleFinitePSpace)
-from sympy.stats.rv import _value_check, Density, RandomSymbol
+from sympy.stats.rv import _value_check, Density, is_random
 
 numpy = import_module('numpy')
 scipy = import_module('scipy')
@@ -214,7 +214,7 @@ def set(self):
 
     def pmf(self, x):
         x = sympify(x)
-        if not (x.is_number or x.is_Symbol or isinstance(x, RandomSymbol)):
+        if not (x.is_number or x.is_Symbol or is_random(x)):
             raise ValueError("'x' expected as an argument of type 'number' or 'Symbol' or , "
                         "'RandomSymbol' not %s" % (type(x)))
         cond = Ge(x, 1) & Le(x, self.sides) & Contains(x, S.Integers)
@@ -396,7 +396,7 @@ def set(self):
     def pmf(self, x):
         n, p = self.n, self.p
         x = sympify(x)
-        if not (x.is_number or x.is_Symbol or isinstance(x, RandomSymbol)):
+        if not (x.is_number or x.is_Symbol or is_random(x)):
             raise ValueError("'x' expected as an argument of type 'number' or 'Symbol' or , "
                         "'RandomSymbol' not %s" % (type(x)))
         cond = Ge(x, 0) & Le(x, n) & Contains(x, S.Integers)
diff --git a/sympy/stats/random_matrix_models.py b/sympy/stats/random_matrix_models.py
index 0735a7bcbc2d..871a92bdf978 100644
--- a/sympy/stats/random_matrix_models.py
+++ b/sympy/stats/random_matrix_models.py
@@ -3,12 +3,12 @@
 from sympy import (Basic, exp, pi, Lambda, Trace, S, MatrixSymbol, Integral,
                    gamma, Product, Dummy, Sum, Abs, IndexedBase, I)
 from sympy.core.sympify import _sympify
-from sympy.stats.rv import (_symbol_converter, Density, RandomMatrixSymbol,
-                            RandomSymbol)
+from sympy.stats.rv import _symbol_converter, Density, RandomMatrixSymbol, is_random
 from sympy.stats.joint_rv_types import JointDistributionHandmade
 from sympy.stats.random_matrix import RandomMatrixPSpace
 from sympy.tensor.array import ArrayComprehension
 
+
 __all__ = [
     'CircularEnsemble',
     'CircularUnitaryEnsemble',
@@ -351,7 +351,7 @@ def JointEigenDistribution(mat):
 
     """
     eigenvals = mat.eigenvals(multiple=True)
-    if any(not eigenval.has(RandomSymbol) for eigenval in set(eigenvals)):
+    if any(not is_random(eigenval) for eigenval in set(eigenvals)):
         raise ValueError("Eigen values don't have any random expression, "
                          "joint distribution cannot be generated.")
     return JointDistributionHandmade(*eigenvals)
diff --git a/sympy/stats/rv.py b/sympy/stats/rv.py
index 2a427915cd47..aa7a3478c07a 100644
--- a/sympy/stats/rv.py
+++ b/sympy/stats/rv.py
@@ -15,11 +15,12 @@
 
 from __future__ import print_function, division
 
+from functools import singledispatch
 from typing import Tuple as tTuple
 
 from sympy import (Basic, S, Expr, Symbol, Tuple, And, Add, Eq, lambdify,
                    Equality, Lambda, sympify, Dummy, Ne, KroneckerDelta,
-                   DiracDelta, Mul, Indexed, MatrixSymbol, Function)
+                   DiracDelta, Mul, Indexed, MatrixSymbol, Function, MatrixBase)
 from sympy.core.relational import Relational
 from sympy.core.sympify import _sympify
 from sympy.logic.boolalg import Boolean
@@ -243,8 +244,10 @@ def __new__(cls, symbol, pspace=None):
         if pspace is None:
             # Allow single arg, representing pspace == PSpace()
             pspace = PSpace()
+        if isinstance(symbol, str):
+            symbol = Symbol(symbol)
         if not isinstance(symbol, Symbol):
-            raise TypeError("symbol should be of type Symbol")
+            raise TypeError("symbol should be of type Symbol or string")
         if not isinstance(pspace, PSpace):
             raise TypeError("pspace variable should be of type PSpace")
         if cls == JointRandomSymbol and isinstance(pspace, SinglePSpace):
@@ -640,7 +643,7 @@ def given(expr, condition=None, **kwargs):
          2*\/ pi
     """
 
-    if not random_symbols(condition) or pspace_independent(expr, condition):
+    if not is_random(condition) or pspace_independent(expr, condition):
         return expr
 
     if isinstance(condition, RandomSymbol):
@@ -718,7 +721,7 @@ def expectation(expr, condition=None, numsamples=None, evaluate=True, **kwargs):
     5
     """
 
-    if not random_symbols(expr):  # expr isn't random?
+    if not is_random(expr):  # expr isn't random?
         return expr
     if numsamples:  # Computing by monte carlo sampling?
         evalf = kwargs.get('evalf', True)
@@ -738,6 +741,10 @@ def expectation(expr, condition=None, numsamples=None, evaluate=True, **kwargs):
                      for arg in expr.args])
 
     # Otherwise case is simple, pass work off to the ProbabilitySpace
+    if pspace(expr) == PSpace():
+        from sympy.stats.symbolic_probability import Expectation
+        return Expectation(expr)
+
     result = pspace(expr).compute_expectation(expr, evaluate=evaluate, **kwargs)
     if evaluate and hasattr(result, 'doit'):
         return result.doit(**kwargs)
@@ -816,6 +823,10 @@ def probability(condition, given_condition=None, numsamples=None,
         return probability(given(condition, given_condition, **kwargs), **kwargs)
 
     # Otherwise pass work off to the ProbabilitySpace
+    if pspace(condition) == PSpace():
+        from sympy.stats.symbolic_probability import Probability
+        return Probability(condition, given_condition)
+
     result = pspace(condition).probability(condition, **kwargs)
     if evaluate and hasattr(result, 'doit'):
         return result.doit()
@@ -1514,3 +1525,40 @@ def _symbol_converter(sym):
     if not isinstance(sym, Symbol):
         raise TypeError("%s is neither a Symbol nor a string"%(sym))
     return sym
+
+
+@singledispatch
+def is_random(x):
+    return False
+
+@is_random.register(RandomSymbol)
+def _(x):
+    return True
+
+@is_random.register(Expr)
+def _(x):
+    atoms = x.free_symbols
+    if len(atoms) == 1 and next(iter(atoms)) == x:
+        return False
+    return any([is_random(i) for i in atoms])
+
+@is_random.register(Boolean)
+def _(x):
+    atoms = x.free_symbols
+    return any([is_random(i) for i in atoms])
+
+@is_random.register(Indexed)
+def _(x):
+    return is_random(x.base)
+
+@is_random.register(MatrixBase)
+def _(x):
+    return any([is_random(i) for i in x])
+
+@is_random.register(RandomMatrixSymbol)
+def _(x):
+    return True
+
+@is_random.register(RandomIndexedSymbol)
+def _(x):
+    return True
diff --git a/sympy/stats/rv_interface.py b/sympy/stats/rv_interface.py
index fa63a6ca429a..f305c61da9d1 100644
--- a/sympy/stats/rv_interface.py
+++ b/sympy/stats/rv_interface.py
@@ -1,7 +1,7 @@
 from __future__ import print_function, division
 from sympy.sets import FiniteSet
 from sympy import sqrt, log, exp, FallingFactorial, Rational, Eq, Dummy, piecewise_fold, solveset
-from .rv import (probability, expectation, density, where, given, pspace, cdf,
+from .rv import (probability, expectation, density, where, given, pspace, cdf, PSpace,
                  characteristic_function, sample, sample_iter, random_symbols, independent, dependent,
                  sampling_density, moment_generating_function, quantile)
 
@@ -57,6 +57,10 @@ def variance(X, condition=None, **kwargs):
     >>> simplify(variance(B))
     p*(1 - p)
     """
+    if pspace(X) == PSpace():
+        from sympy.stats.symbolic_probability import Variance
+        return Variance(X, condition)
+
     return cmoment(X, 2, condition, **kwargs)
 
 
@@ -148,6 +152,10 @@ def covariance(X, Y, condition=None, **kwargs):
     >>> covariance(X, Y + rate*X)
     1/lambda
     """
+    if pspace(X) == PSpace() or pspace(Y) == PSpace():
+        from sympy.stats.symbolic_probability import Covariance
+        return Covariance(X, Y, condition)
+
     return expectation(
         (X - expectation(X, condition, **kwargs)) *
         (Y - expectation(Y, condition, **kwargs)),
diff --git a/sympy/stats/symbolic_probability.py b/sympy/stats/symbolic_probability.py
index 0f114bc1c39c..4bd53d122724 100644
--- a/sympy/stats/symbolic_probability.py
+++ b/sympy/stats/symbolic_probability.py
@@ -5,7 +5,8 @@
 from sympy.core.parameters import global_parameters
 from sympy.core.sympify import _sympify
 from sympy.stats import variance, covariance
-from sympy.stats.rv import RandomSymbol, probability, expectation
+from sympy.stats.rv import RandomSymbol, probability, expectation, is_random
+
 
 __all__ = ['Probability', 'Expectation', 'Variance', 'Covariance']
 
@@ -102,7 +103,7 @@ class Expectation(Expr):
     def __new__(cls, expr, condition=None, **kwargs):
         expr = _sympify(expr)
         if condition is None:
-            if not expr.has(RandomSymbol):
+            if not is_random(expr):
                 return expr
             obj = Expr.__new__(cls, expr)
         else:
@@ -115,7 +116,7 @@ def doit(self, **hints):
         expr = self.args[0]
         condition = self._condition
 
-        if not expr.has(RandomSymbol):
+        if not is_random(expr):
             return expr
 
         if isinstance(expr, Add):
@@ -124,7 +125,7 @@ def doit(self, **hints):
             rv = []
             nonrv = []
             for a in expr.args:
-                if isinstance(a, RandomSymbol) or a.has(RandomSymbol):
+                if is_random(a):
                     rv.append(a)
                 else:
                     nonrv.append(a)
@@ -229,7 +230,7 @@ def doit(self, **hints):
         arg = self.args[0]
         condition = self._condition
 
-        if not arg.has(RandomSymbol):
+        if not is_random(arg):
             return S.Zero
 
         if isinstance(arg, RandomSymbol):
@@ -237,7 +238,7 @@ def doit(self, **hints):
         elif isinstance(arg, Add):
             rv = []
             for a in arg.args:
-                if a.has(RandomSymbol):
+                if is_random(a):
                     rv.append(a)
             variances = Add(*map(lambda xv: Variance(xv, condition).doit(), rv))
             map_to_covar = lambda x: 2*Covariance(*x, condition=condition).doit()
@@ -247,7 +248,7 @@ def doit(self, **hints):
             nonrv = []
             rv = []
             for a in arg.args:
-                if a.has(RandomSymbol):
+                if is_random(a):
                     rv.append(a)
                 else:
                     nonrv.append(a**2)
@@ -343,9 +344,9 @@ def doit(self, **hints):
         if arg1 == arg2:
             return Variance(arg1, condition).doit()
 
-        if not arg1.has(RandomSymbol):
+        if not is_random(arg1):
             return S.Zero
-        if not arg2.has(RandomSymbol):
+        if not is_random(arg2):
             return S.Zero
 
         arg1, arg2 = sorted([arg1, arg2], key=default_sort_key)
@@ -370,13 +371,13 @@ def _expand_single_argument(cls, expr):
             for a in expr.args:
                 if isinstance(a, Mul):
                     outval.append(cls._get_mul_nonrv_rv_tuple(a))
-                elif isinstance(a, RandomSymbol):
+                elif is_random(a):
                     outval.append((S.One, a))
 
             return outval
         elif isinstance(expr, Mul):
             return [cls._get_mul_nonrv_rv_tuple(expr)]
-        elif expr.has(RandomSymbol):
+        elif is_random(expr):
             return [(S.One, expr)]
 
     @classmethod
@@ -384,7 +385,7 @@ def _get_mul_nonrv_rv_tuple(cls, m):
         rv = []
         nonrv = []
         for a in m.args:
-            if a.has(RandomSymbol):
+            if is_random(a):
                 rv.append(a)
             else:
                 nonrv.append(a)
diff --git a/sympy/stats/tests/test_rv.py b/sympy/stats/tests/test_rv.py
index 6d91e3097d4a..035bc01e31be 100644
--- a/sympy/stats/tests/test_rv.py
+++ b/sympy/stats/tests/test_rv.py
@@ -1,13 +1,14 @@
 from __future__ import unicode_literals
 from sympy import (S, Symbol, Interval, exp,
         symbols, Eq, cos, And, Tuple, integrate, oo, sin, Sum, Basic,
-        DiracDelta, Lambda, log, pi, FallingFactorial, Rational)
+        DiracDelta, Lambda, log, pi, FallingFactorial, Rational, Matrix)
 from sympy.stats import (Die, Normal, Exponential, FiniteRV, P, E, H, variance,
-        density, given, independent, dependent, where, pspace,
+        density, given, independent, dependent, where, pspace, GaussianUnitaryEnsemble,
         random_symbols, sample, Geometric, factorial_moment, Binomial, Hypergeometric,
-        DiscreteUniform, Poisson, characteristic_function, moment_generating_function)
+        DiscreteUniform, Poisson, characteristic_function, moment_generating_function,
+        BernoulliProcess, Variance, Expectation, Probability, Covariance, covariance)
 from sympy.stats.rv import (IndependentProductPSpace, rs_swap, Density, NamedArgsMixin,
-        RandomSymbol, sample_iter, PSpace)
+        RandomSymbol, sample_iter, PSpace, is_random)
 from sympy.testing.pytest import raises, skip, XFAIL, ignore_warnings
 from sympy.external import import_module
 from sympy.core.numbers import comp
@@ -338,3 +339,33 @@ def test_issue_12237():
     assert W == P(X + Y > 0, X)
     assert U == BernoulliDistribution(S.Half, S.Zero, S.One)
     assert V == S.Half
+
+def test_is_random():
+    X = Normal('X', 0, 1)
+    Y = Normal('Y', 0, 1)
+    a, b = symbols('a, b')
+    G = GaussianUnitaryEnsemble('U', 2)
+    B = BernoulliProcess('B', 0.9)
+    assert not is_random(a)
+    assert not is_random(a + b)
+    assert not is_random(a * b)
+    assert not is_random(Matrix([a**2, b**2]))
+    assert is_random(X)
+    assert is_random(X**2 + Y)
+    assert is_random(Y + b**2)
+    assert is_random(Y > 5)
+    assert is_random(B[3] < 1)
+    assert is_random(G)
+    assert is_random(X * Y * B[1])
+    assert is_random(Matrix([[X, B[2]], [G, Y]]))
+    assert is_random(Eq(X, 4))
+
+def test_issue_12283():
+    x = symbols('x')
+    X = RandomSymbol(x)
+    Y = RandomSymbol('Y')
+    assert pspace(X) == PSpace()
+    assert E(X) == Expectation(X)
+    assert P(X > 3) == Probability(X > 3)
+    assert variance(X) == Variance(X)
+    assert covariance(X, Y) == Covariance(X, Y)
diff --git a/sympy/stats/tests/test_symbolic_probability.py b/sympy/stats/tests/test_symbolic_probability.py
index 4dca7d1e3d77..01ae5b899d6c 100644
--- a/sympy/stats/tests/test_symbolic_probability.py
+++ b/sympy/stats/tests/test_symbolic_probability.py
@@ -66,7 +66,6 @@ def test_literal_probability():
     assert Covariance(X, X**2).doit() == Covariance(X, X**2)
     assert Covariance(X, sin(X)).doit() == Covariance(sin(X), X)
     assert Covariance(X**2, sin(X)*Y).doit() == Covariance(sin(X)*Y, X**2)
-    assert Covariance(w, X).evaluate_integral() == 0
 
 
 def test_probability_rewrite():
@@ -102,9 +101,5 @@ def test_probability_rewrite():
     assert Variance(X).rewrite(Sum) == Variance(X).rewrite(Integral)
     assert Expectation(X).rewrite(Sum) == Expectation(X).rewrite(Integral)
 
-    assert Covariance(w, X).rewrite(Sum) == 0
-
-    assert Covariance(w, X).rewrite(Integral) == 0
-
     assert Variance(X, condition=Y).rewrite(Probability) == Integral(x**2*Probability(Eq(X, x), Y), (x, -oo, oo)) - \
                                                             Integral(x*Probability(Eq(X, x), Y), (x, -oo, oo))**2


Code review comments:
diff_hunk: @@ -66,7 +66,6 @@ def test_literal_probability():
     assert Covariance(X, X**2).doit() == Covariance(X, X**2)
     assert Covariance(X, sin(X)).doit() == Covariance(sin(X), X)
     assert Covariance(X**2, sin(X)*Y).doit() == Covariance(sin(X)*Y, X**2)
-    assert Covariance(w, X).evaluate_integral() == 0
path: sympy/stats/tests/test_symbolic_probability.py
text: @user1:
Were the deleted results wrong?

@author:
Oops, will add them. Actually while debugging an error, I deleted it mistakenly