Repository: sympy/sympy
Pull request: [GSoC]Change return type of `P` and `E` with `evaluate=False`
<!-- 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. -->
Fixes #19798

#### Brief description of what is fixed or changed
Change in the return types of P and E with `evaluate=False`. Examples:
```py
>>> from sympy.stats import *
>>> X = Normal('X',1,2)
>>> P(X > 1, evaluate=False)
Probability(X > 1)
>>> E(X, evaluate=False)
Expectation(X)
```

#### 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
   * Change in return type of `P` and `E` with `evaluate=False`. With `evaluate=False`, `P` and `E` are made to return `Probability` and `Expectation` object respectively.
<!-- END RELEASE NOTES -->
ping @czgdp1807 @Upabjojr @jmig5776 

Problem statement:
stats: Should `P(X > 1, evaluate=False)` return `Probability(X > 1)` or unevaluated Integral?
On current master, `P(X > 1, evaluate=False)` produces unevaluated Integral. However, according to the discussions on Gitter, it was suggested by @czgdp1807 that we should make `P(X > 1, evaluate=False)` to return a Symbolic probability object as `Probability(X > 1)`. Symbolic `Probability` class has been added with `doit` method recently and so, `P(X > 1, evaluate=True)` should return `Probability(X > 1).doit()`.
More clearly the current scenario is:
```py
In [1]: from sympy.stats import *

In [2]: X = Normal('X', 1, 2)

In [3]: P(X > 1, evaluate=False)
Out[3]:

∞
⌠
⎮             2
⎮     -(z - 1)
⎮     ──────────
⎮         8
⎮ √2⋅ℯ
⎮ ────────────── dz
⎮      4⋅√π
⌡
1

```
Proposal to change the behavior to:

```py
In [5]: P(X > 1, evaluate=False)
Out[5]: Probability(X > 1)
```
ping @czgdp1807 @Upabjojr

Hints from issue discussion:
What about:

1. `P(X > 1, evaluate=False)` ==> return `Probability` object,
2. `P(X > 1).doit(evaluate=False)` ==> return unevaluated integral.
Agreed with @Upabjojr , seems like a good idea and will also help in making APIs unambiguous. 
> `P(X > 1).doit(evaluate=False)` ==> return unevaluated integral.

Do you mean `Probability(X > 1).doit(evaluate=False)` ==> return unevaluated integral ?
as `P(X > 1)` will already be evaluated(as `evaluate=True` is by default) in `P`
I think he meant to say what you are suggesting. What `doit` does is a different issue, what we need right now is `evaluate=False` to return `Probability(X > 1)`. Output of `doit` depends on how much the input can be processed before giving results and return type is variable so we cannot do much about it.
Maybe we don't need to deal with `.doit()` at all. We already have `Probability(X > 1).rewrite(Integral)`... which is supposed to return the `Integral` object.
Maybe `P(X > 1, evaluate=False)` should raise a warning about the change.

Commit message: remove import Integral

Patch to review:
diff --git a/sympy/stats/rv.py b/sympy/stats/rv.py
index 2d2c6f67ada1..b84c6a895b5c 100644
--- a/sympy/stats/rv.py
+++ b/sympy/stats/rv.py
@@ -20,7 +20,7 @@
 
 from sympy import (Basic, S, Expr, Symbol, Tuple, And, Add, Eq, lambdify, Or,
                    Equality, Lambda, sympify, Dummy, Ne, KroneckerDelta,
-                   DiracDelta, Mul, Indexed, MatrixSymbol, Function, Integral)
+                   DiracDelta, Mul, Indexed, MatrixSymbol, Function)
 from sympy.core.relational import Relational
 from sympy.core.sympify import _sympify
 from sympy.sets.sets import FiniteSet, ProductSet, Intersection
@@ -757,7 +757,11 @@ def expectation(expr, condition=None, numsamples=None, evaluate=True, **kwargs):
     from sympy.stats.symbolic_probability import Expectation
     if evaluate:
         return Expectation(expr, condition).doit(**kwargs)
-    return Expectation(expr, condition).rewrite(Integral) # will return Sum in case of discrete RV
+    message = ("Using `evaluate=False` now returns `Expectation` object "
+              "since version 1.7. If you want unevaluated Integral/Sum use "
+              "`E(expr, condition, evaluate=False).rewrite(Integral)`")
+    warnings.warn(filldedent(message))
+    return Expectation(expr, condition)
 
 
 def probability(condition, given_condition=None, numsamples=None,
@@ -796,7 +800,11 @@ def probability(condition, given_condition=None, numsamples=None,
     from sympy.stats.symbolic_probability import Probability
     if evaluate:
         return Probability(condition, given_condition).doit(**kwargs)
-    return Probability(condition, given_condition).rewrite(Integral) # will return Sum in case of discrete RV
+    message = ("Using `evaluate=False` returns `Probability` object since version 1.7. "
+              "If you want unevaluated Integral/Sum use "
+              "`P(condition, given_condition, evaluate=False).rewrite(Integral)`")
+    warnings.warn(filldedent(message))
+    return Probability(condition, given_condition)
 
 
 class Density(Basic):
diff --git a/sympy/stats/tests/test_compound_rv.py b/sympy/stats/tests/test_compound_rv.py
index 5b32abeb1e45..3b88bb1a1344 100644
--- a/sympy/stats/tests/test_compound_rv.py
+++ b/sympy/stats/tests/test_compound_rv.py
@@ -6,7 +6,7 @@
 from sympy.stats.crv_types import NormalDistribution
 from sympy.stats.drv_types import PoissonDistribution
 from sympy.stats.frv_types import BernoulliDistribution
-from sympy.testing.pytest import raises
+from sympy.testing.pytest import raises, ignore_warnings
 from sympy.stats.joint_rv_types import MultivariateNormalDistribution
 
 
@@ -66,7 +66,8 @@ def test_unevaluated_CompoundDist():
     expre = Integral(Piecewise((_k*exp(S(3)/4 - _k/4)/8, 2*Abs(arg(_k - 3)) <= pi/2),
     (sqrt(2)*_k*Integral(exp(-(_k**4 + 16*(_k - 3)**2)/(32*_k**2)),
     (_k, 0, oo))/(32*sqrt(pi)), True)), (_k, -oo, oo))
-    assert (E(X, evaluate=False).simplify()).dummy_eq(expre.simplify())
+    with ignore_warnings(UserWarning):
+        assert (E(X, evaluate=False).rewrite(Integral).simplify()).dummy_eq(expre.simplify())
 
 
 def test_Compound_Distribution():
diff --git a/sympy/stats/tests/test_continuous_rv.py b/sympy/stats/tests/test_continuous_rv.py
index 1e249c68eaa8..a978c17e2f95 100644
--- a/sympy/stats/tests/test_continuous_rv.py
+++ b/sympy/stats/tests/test_continuous_rv.py
@@ -25,6 +25,7 @@
 from sympy.stats.joint_rv_types import MultivariateLaplaceDistribution, MultivariateNormalDistribution
 from sympy.stats.crv import SingleContinuousPSpace, SingleContinuousDomain
 from sympy.stats.compound_rv import CompoundPSpace
+from sympy.stats.symbolic_probability import Probability
 from sympy.testing.pytest import raises, XFAIL, slow, skip, ignore_warnings
 from sympy.testing.randtest import verify_numerically as tn
 
@@ -643,12 +644,12 @@ def test_exponential():
     _z = Dummy('_z')
     b = SingleContinuousPSpace(x, ExponentialDistribution(2))
 
-    expected1 = Integral(2*exp(-2*_z), (_z, 3, oo))
-    assert b.probability(x > 3, evaluate=False).dummy_eq(expected1) is True
-
-    expected2 = Integral(2*exp(-2*_z), (_z, 0, 4))
-    assert b.probability(x < 4, evaluate=False).dummy_eq(expected2) is True
+    with ignore_warnings(UserWarning):
+        expected1 = Integral(2*exp(-2*_z), (_z, 3, oo))
+        assert b.probability(x > 3, evaluate=False).rewrite(Integral).dummy_eq(expected1)
 
+        expected2 = Integral(2*exp(-2*_z), (_z, 0, 4))
+        assert b.probability(x < 4, evaluate=False).rewrite(Integral).dummy_eq(expected2)
     Y = Exponential('y', 2*rate)
     assert coskewness(X, X, X) == skewness(X)
     assert coskewness(X, Y + rate*X, Y + 2*rate*X) == \
@@ -1342,21 +1343,23 @@ def test_input_value_assertions():
 
 def test_unevaluated():
     X = Normal('x', 0, 1)
-    assert str(E(X, evaluate=False)) == ("Integral(sqrt(2)*x*exp(-x**2/2)/"
-    "(2*sqrt(pi)), (x, -oo, oo))")
+    with ignore_warnings(UserWarning):
+        assert str(E(X, evaluate=False).rewrite(Integral)) == ("Integral(sqrt(2)*x*exp(-x**2/2)/"
+        "(2*sqrt(pi)), (x, -oo, oo))")
 
-    assert str(E(X + 1, evaluate=False)) == ("Integral(sqrt(2)*x*exp(-x**2/2)/"
-    "(2*sqrt(pi)), (x, -oo, oo)) + 1")
+        assert str(E(X + 1, evaluate=False).rewrite(Integral)) == ("Integral(sqrt(2)*x*exp(-x**2/2)/"
+        "(2*sqrt(pi)), (x, -oo, oo)) + 1")
 
-    assert str(P(X > 0, evaluate=False)) == ("Integral(sqrt(2)*exp(-_z**2/2)/"
-    "(2*sqrt(pi)), (_z, 0, oo))")
+        assert str(P(X > 0, evaluate=False).rewrite(Integral)) == ("Integral(sqrt(2)*exp(-_z**2/2)/"
+        "(2*sqrt(pi)), (_z, 0, oo))")
 
-    assert P(X > 0, X**2 < 1, evaluate=False) == S.Half
+    assert P(X > 0, X**2 < 1) == S.Half
 
 
 def test_probability_unevaluated():
     T = Normal('T', 30, 3)
-    assert type(P(T > 33, evaluate=False)) == Integral
+    with ignore_warnings(UserWarning):
+        assert type(P(T > 33, evaluate=False)) == Probability
 
 
 def test_density_unevaluated():
diff --git a/sympy/stats/tests/test_discrete_rv.py b/sympy/stats/tests/test_discrete_rv.py
index b03b67c4686b..1c4e94070b29 100644
--- a/sympy/stats/tests/test_discrete_rv.py
+++ b/sympy/stats/tests/test_discrete_rv.py
@@ -15,6 +15,7 @@
 from sympy.stats.rv import sample
 from sympy.testing.pytest import slow, nocache_fail, raises, skip, ignore_warnings
 from sympy.external import import_module
+from sympy.stats.symbolic_probability import Expectation
 
 x = Symbol('x')
 
@@ -33,8 +34,9 @@ def test_Poisson():
     assert E(x) == l
     assert variance(x) == l
     assert density(x) == PoissonDistribution(l)
-    assert isinstance(E(x, evaluate=False), Sum)
-    assert isinstance(E(2*x, evaluate=False), Sum)
+    with ignore_warnings(UserWarning):
+        assert isinstance(E(x, evaluate=False), Expectation)
+        assert isinstance(E(2*x, evaluate=False), Expectation)
     # issue 8248
     assert x.pspace.compute_expectation(1) == 1
 
@@ -82,7 +84,8 @@ def test_Logarithmic():
     assert E(x) == -p / ((1 - p) * log(1 - p))
     assert variance(x) == -1/log(2)**2 + 2/log(2)
     assert E(2*x**2 + 3*x + 4) == 4 + 7 / log(2)
-    assert isinstance(E(x, evaluate=False), Sum)
+    with ignore_warnings(UserWarning):
+        assert isinstance(E(x, evaluate=False), Expectation)
 
 
 @nocache_fail
@@ -94,7 +97,8 @@ def test_negative_binomial():
     # This hangs when run with the cache disabled:
     assert variance(x) == p*r / (1-p)**2
     assert E(x**5 + 2*x + 3) == Rational(9207, 4)
-    assert isinstance(E(x, evaluate=False), Sum)
+    with ignore_warnings(UserWarning):
+        assert isinstance(E(x, evaluate=False), Expectation)
 
 
 def test_skellam():
@@ -121,7 +125,8 @@ def test_yule_simon():
     x = YuleSimon('x', rho)
     assert simplify(E(x)) == rho / (rho - 1)
     assert simplify(variance(x)) == rho**2 / ((rho - 1)**2 * (rho - 2))
-    assert isinstance(E(x, evaluate=False), Sum)
+    with ignore_warnings(UserWarning):
+        assert isinstance(E(x, evaluate=False), Expectation)
     # To test the cdf function
     assert cdf(x)(x) == Piecewise((-beta(floor(x), 4)*floor(x) + 1, x >= 1), (0, True))
 
@@ -272,7 +277,8 @@ def test_product_spaces():
     #assert str(P(X1 + X2 < 3, evaluate=False)) == """Sum(Piecewise((2**(X2 - n - 2)*(2/3)**(X2 - 1)/6, """\
     #    + """(-X2 + n + 3 >= 1) & (-X2 + n + 3 < oo)), (0, True)), (X2, 1, oo), (n, -oo, -1))"""
     n = Dummy('n')
-    assert P(X1 + X2 < 3, evaluate=False).dummy_eq(Sum(Piecewise((2**(-n)/4,
+    with ignore_warnings(UserWarning):
+        assert P(X1 + X2 < 3, evaluate=False).rewrite(Sum).dummy_eq(Sum(Piecewise((2**(-n)/4,
          n + 2 >= 1), (0, True)), (n, -oo, -1))/3)
     #assert str(P(X1 + X2 > 3)) == """Sum(Piecewise((2**(X2 - n - 2)*(2/3)**(X2 - 1)/6, """ +\
     #    """(-X2 + n + 3 >= 1) & (-X2 + n + 3 < oo)), (0, True)), (X2, 1, oo), (n, 1, oo))"""
diff --git a/sympy/stats/tests/test_mix.py b/sympy/stats/tests/test_mix.py
index 9cbd336646b5..b2686f26093e 100644
--- a/sympy/stats/tests/test_mix.py
+++ b/sympy/stats/tests/test_mix.py
@@ -1,6 +1,6 @@
 from sympy import (Symbol, Eq, Ne, simplify, sqrt, exp, pi, symbols,
                 Piecewise, factorial, gamma, IndexedBase, Add, Pow, Mul,
-                Indexed, Integer)
+                Indexed, Integer, Integral)
 from sympy.functions.elementary.piecewise import ExprCondPair
 from sympy.stats import (Poisson, Beta, Exponential, P,
                         Multinomial, MultivariateBeta)
@@ -9,6 +9,7 @@
 from sympy.stats.compound_rv import CompoundPSpace, CompoundDistribution
 from sympy.stats.joint_rv import MarginalDistribution
 from sympy.stats.rv import pspace, density
+from sympy.testing.pytest import ignore_warnings
 
 def test_density():
     x = Symbol('x')
@@ -58,7 +59,8 @@ def test_mix_expression():
     Y, E = Poisson('Y', 1), Exponential('E', 1)
     assert P(Eq(Y + E, 1)) == 0
     assert P(Ne(Y + E, 2)) == 1
-    assert str(P(E + Y < 2, evaluate=False)) == """Integral(Sum(exp(-1)*Integral"""\
-+"""(exp(-E)*DiracDelta(-_z + E + Y - 2), (E, 0, oo))/factorial(Y), (Y, 0, oo)), (_z, -oo, 0))"""
-    assert str(P(E + Y > 2, evaluate=False)) == """Integral(Sum(exp(-1)*Integral"""\
-+"""(exp(-E)*DiracDelta(-_z + E + Y - 2), (E, 0, oo))/factorial(Y), (Y, 0, oo)), (_z, 0, oo))"""
+    with ignore_warnings(UserWarning):
+        assert str(P(E + Y < 2, evaluate=False).rewrite(Integral)) == """Integral(Sum(exp(-1)*Integral"""\
+    +"""(exp(-E)*DiracDelta(-_z + E + Y - 2), (E, 0, oo))/factorial(Y), (Y, 0, oo)), (_z, -oo, 0))"""
+        assert str(P(E + Y > 2, evaluate=False).rewrite(Integral)) == """Integral(Sum(exp(-1)*Integral"""\
+    +"""(exp(-E)*DiracDelta(-_z + E + Y - 2), (E, 0, oo))/factorial(Y), (Y, 0, oo)), (_z, 0, oo))"""
diff --git a/sympy/stats/tests/test_stochastic_process.py b/sympy/stats/tests/test_stochastic_process.py
index f4ff19a0e14d..a3c281c92c65 100644
--- a/sympy/stats/tests/test_stochastic_process.py
+++ b/sympy/stats/tests/test_stochastic_process.py
@@ -152,7 +152,7 @@ def test_ContinuousMarkovChain():
     assert C2.transition_probabilities(A)(t) == Matrix([[S.Half + exp(-2*t)/2, S.Half - exp(-2*t)/2, 0],
                                                        [S.Half - exp(-2*t)/2, S.Half + exp(-2*t)/2, 0],
                                                        [S.Half - exp(-t) + exp(-2*t)/2, S.Half - exp(-2*t)/2, exp(-t)]])
-    assert P(Eq(C2(1), 1), Eq(C2(0), 1), evaluate=False) == Probability(Eq(C2(1), 1))
+    assert P(Eq(C2(1), 1), Eq(C2(0), 1), evaluate=False) == Probability(Eq(C2(1), 1), Eq(C2(0), 1))
     assert P(Eq(C2(1), 1), Eq(C2(0), 1)) == exp(-2)/2 + S.Half
     assert P(Eq(C2(1), 0) & Eq(C2(2), 1) & Eq(C2(3), 1),
                 Eq(P(Eq(C2(1), 0)), S.Half)) == (Rational(1, 4) - exp(-2)/4)*(exp(-2)/2 + S.Half)


Code review comments:
diff_hunk: @@ -757,7 +757,11 @@ def expectation(expr, condition=None, numsamples=None, evaluate=True, **kwargs):
     from sympy.stats.symbolic_probability import Expectation
     if evaluate:
         return Expectation(expr, condition).doit(**kwargs)
-    return Expectation(expr, condition).rewrite(Integral) # will return Sum in case of discrete RV
+    message = ("Using `evaluate=False` now returns `Expectation` object "
+              "since version 1.7. If you want unevaluated Integral/Sum use "
+              "`E(expr, condition, evaluate=False).rewrite(Integral)`")
path: sympy/stats/rv.py
text: @user1:
```suggestion
    message = ("Since version 1.7, using `evaluate=False` returns `Expectation` object. "
              "If you want unevaluated Integral/Sum use "
              "`E(expr, condition, evaluate=False).rewrite(Integral)`")
```

diff_hunk: @@ -796,7 +800,11 @@ def probability(condition, given_condition=None, numsamples=None,
     from sympy.stats.symbolic_probability import Probability
     if evaluate:
         return Probability(condition, given_condition).doit(**kwargs)
-    return Probability(condition, given_condition).rewrite(Integral) # will return Sum in case of discrete RV
+    message = ("Using `evaluate=False` returns `Probability` object since version 1.7. "
+              "If you want unevaluated Integral/Sum use "
+              "`P(condition, given_condition, evaluate=False).rewrite(Integral)`")
path: sympy/stats/rv.py
text: @user1:
```suggestion
    message = ("Since version version 1.7, using `evaluate=False` returns `Probability` object. "
              "If you want unevaluated Integral/Sum use "
              "`P(condition, given_condition, evaluate=False).rewrite(Integral)`")
```