Repository: sympy/sympy

Issue:
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 

Discussion:

Patch:
diff --git a/sympy/stats/rv.py b/sympy/stats/rv.py
index 2d2c6f67ad..c1bc926766 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,12 @@ 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
+    ### TODO: Remove the user warnings in the future releases
+    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)`")
+    warnings.warn(filldedent(message))
+    return Expectation(expr, condition)
 
 
 def probability(condition, given_condition=None, numsamples=None,
@@ -796,7 +801,12 @@ 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
+    ### TODO: Remove the user warnings in the future releases
+    message = ("Since version 1.7, using `evaluate=False` returns `Probability` "
+              "object. 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):
@@ -1048,6 +1058,7 @@ def sample(expr, condition=None, size=(), library='scipy', numsamples=1,
         iterator object containing the sample/samples of given expr
 
     """
+    ### TODO: Remove the user warnings in the future releases
     message = ("The return type of sample has been changed to return an "
                   "iterator object since version 1.7. For more information see "
                   "https://github.com/sympy/sympy/issues/19061")


Test patch:
diff --git a/sympy/stats/tests/test_compound_rv.py b/sympy/stats/tests/test_compound_rv.py
index 5b32abeb1e..40b6616594 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): ### TODO: Restore tests once warnings are removed
+        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 1e249c68ea..5b32b68878 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
 
@@ -327,7 +328,7 @@ def test_sample_continuous():
     scipy = import_module('scipy')
     if not scipy:
         skip('Scipy is not installed. Abort tests')
-    with ignore_warnings(UserWarning):
+    with ignore_warnings(UserWarning): ### TODO: Restore tests once warnings are removed
         assert next(sample(Z)) in Z.pspace.domain.set
     sym, val = list(Z.pspace.sample().items())[0]
     assert sym == Z and val in Interval(0, oo)
@@ -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): ### TODO: Restore tests once warnings are removed
+        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) == \
@@ -758,7 +759,7 @@ def test_sampling_gamma_inverse():
     if not scipy:
         skip('Scipy not installed. Abort tests for sampling of gamma inverse.')
     X = GammaInverse("x", 1, 1)
-    with ignore_warnings(UserWarning):
+    with ignore_warnings(UserWarning): ### TODO: Restore tests once warnings are removed
         assert next(sample(X)) in X.pspace.domain.set
 
 def test_gompertz():
@@ -885,13 +886,13 @@ def test_lognormal():
     scipy = import_module('scipy')
     if not scipy:
         skip('Scipy is not installed. Abort tests')
-    with ignore_warnings(UserWarning):
+    with ignore_warnings(UserWarning): ### TODO: Restore tests once warnings are removed
         for i in range(3):
             X = LogNormal('x', i, 1)
             assert next(sample(X)) in X.pspace.domain.set
 
     size = 5
-    with ignore_warnings(UserWarning):
+    with ignore_warnings(UserWarning): ### TODO: Restore tests once warnings are removed
         samps = next(sample(X, size=size))
         for samp in samps:
             assert samp in X.pspace.domain.set
@@ -1014,7 +1015,7 @@ def test_sampling_gaussian_inverse():
     if not scipy:
         skip('Scipy not installed. Abort tests for sampling of Gaussian inverse.')
     X = GaussianInverse("x", 1, 1)
-    with ignore_warnings(UserWarning):
+    with ignore_warnings(UserWarning): ### TODO: Restore tests once warnings are removed
         assert next(sample(X, library='scipy')) in X.pspace.domain.set
 
 def test_pareto():
@@ -1316,7 +1317,7 @@ def test_prefab_sampling():
     variables = [N, L, E, P, W, U, B, G]
     niter = 10
     size = 5
-    with ignore_warnings(UserWarning):
+    with ignore_warnings(UserWarning): ### TODO: Restore tests once warnings are removed
         for var in variables:
             for i in range(niter):
                 assert next(sample(var)) in var.pspace.domain.set
@@ -1342,21 +1343,21 @@ 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))")
-
-    assert str(E(X + 1, evaluate=False)) == ("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))")
+    k = Dummy('k')
+    expr1 = Integral(sqrt(2)*k*exp(-k**2/2)/(2*sqrt(pi)), (k, -oo, oo))
+    expr2 = Integral(sqrt(2)*exp(-k**2/2)/(2*sqrt(pi)), (k, 0, oo))
+    with ignore_warnings(UserWarning): ### TODO: Restore tests once warnings are removed
+        assert E(X, evaluate=False).rewrite(Integral).dummy_eq(expr1)
+        assert E(X + 1, evaluate=False).rewrite(Integral).dummy_eq(expr1 + 1)
+        assert P(X > 0, evaluate=False).rewrite(Integral).dummy_eq(expr2)
 
-    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): ### TODO: Restore tests once warnings are removed
+        assert type(P(T > 33, evaluate=False)) == Probability
 
 
 def test_density_unevaluated():
@@ -1583,7 +1584,7 @@ def test_sample_numpy():
     if not numpy:
         skip('Numpy is not installed. Abort tests for _sample_numpy.')
     else:
-        with ignore_warnings(UserWarning):
+        with ignore_warnings(UserWarning): ### TODO: Restore tests once warnings are removed
             for X in distribs_numpy:
                 samps = next(sample(X, size=size, library='numpy'))
                 for sam in samps:
@@ -1616,7 +1617,7 @@ def test_sample_scipy():
     if not scipy:
         skip('Scipy is not installed. Abort tests for _sample_scipy.')
     else:
-        with ignore_warnings(UserWarning):
+        with ignore_warnings(UserWarning): ### TODO: Restore tests once warnings are removed
             g_sample = list(sample(Gamma("G", 2, 7), size=size, numsamples=numsamples))
             assert len(g_sample) == numsamples
             for X in distribs_scipy:
@@ -1646,7 +1647,7 @@ def test_sample_pymc3():
     if not pymc3:
         skip('PyMC3 is not installed. Abort tests for _sample_pymc3.')
     else:
-        with ignore_warnings(UserWarning):
+        with ignore_warnings(UserWarning): ### TODO: Restore tests once warnings are removed
             for X in distribs_pymc3:
                 samps = next(sample(X, size=size, library='pymc3'))
                 for sam in samps:
diff --git a/sympy/stats/tests/test_discrete_rv.py b/sympy/stats/tests/test_discrete_rv.py
index b03b67c468..30388483e6 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): ### TODO: Restore tests once warnings are removed
+        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): ### TODO: Restore tests once warnings are removed
+        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): ### TODO: Restore tests once warnings are removed
+        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): ### TODO: Restore tests once warnings are removed
+        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))
 
@@ -140,7 +145,7 @@ def test_sample_discrete():
     scipy = import_module('scipy')
     if not scipy:
         skip('Scipy not installed. Abort tests')
-    with ignore_warnings(UserWarning):
+    with ignore_warnings(UserWarning): ### TODO: Restore tests once warnings are removed
         assert next(sample(X)) in X.pspace.domain.set
         samps = next(sample(X, size=2)) # This takes long time if ran without scipy
         for samp in samps:
@@ -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): ### TODO: Restore tests once warnings are removed
+        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))"""
@@ -295,7 +301,7 @@ def test_sample_numpy():
     if not numpy:
         skip('Numpy is not installed. Abort tests for _sample_numpy.')
     else:
-        with ignore_warnings(UserWarning):
+        with ignore_warnings(UserWarning): ### TODO: Restore tests once warnings are removed
             for X in distribs_numpy:
                 samps = next(sample(X, size=size, library='numpy'))
                 for sam in samps:
@@ -325,7 +331,7 @@ def test_sample_scipy():
     if not scipy:
         skip('Scipy is not installed. Abort tests for _sample_scipy.')
     else:
-        with ignore_warnings(UserWarning):
+        with ignore_warnings(UserWarning): ### TODO: Restore tests once warnings are removed
             z_sample = list(sample(Zeta("G", 7), size=size, numsamples=numsamples))
             assert len(z_sample) == numsamples
             for X in distribs_scipy:
@@ -348,7 +354,7 @@ def test_sample_pymc3():
     if not pymc3:
         skip('PyMC3 is not installed. Abort tests for _sample_pymc3.')
     else:
-        with ignore_warnings(UserWarning):
+        with ignore_warnings(UserWarning): ### TODO: Restore tests once warnings are removed
             for X in distribs_pymc3:
                 samps = next(sample(X, size=size, library='pymc3'))
                 for sam in samps:
diff --git a/sympy/stats/tests/test_finite_rv.py b/sympy/stats/tests/test_finite_rv.py
index 92f7749132..5598df07b4 100644
--- a/sympy/stats/tests/test_finite_rv.py
+++ b/sympy/stats/tests/test_finite_rv.py
@@ -149,7 +149,7 @@ def test_given():
     scipy = import_module('scipy')
     if not scipy:
         skip('Scipy is not installed. Abort tests')
-    with ignore_warnings(UserWarning):
+    with ignore_warnings(UserWarning): ### TODO: Restore tests once warnings are removed
         assert next(sample(X, X > 5)) == 6
 
 
@@ -452,7 +452,7 @@ def test_sample_numpy():
     if not numpy:
         skip('Numpy is not installed. Abort tests for _sample_numpy.')
     else:
-        with ignore_warnings(UserWarning):
+        with ignore_warnings(UserWarning): ### TODO: Restore tests once warnings are removed
             for X in distribs_numpy:
                 samps = next(sample(X, size=size, library='numpy'))
                 for sam in samps:
@@ -480,7 +480,7 @@ def test_sample_scipy():
     if not scipy:
         skip('Scipy not installed. Abort tests for _sample_scipy.')
     else:
-        with ignore_warnings(UserWarning):
+        with ignore_warnings(UserWarning): ### TODO: Restore tests once warnings are removed
             h_sample = list(sample(Hypergeometric("H", 1, 1, 1), size=size, numsamples=numsamples))
             assert len(h_sample) == numsamples
             for X in distribs_scipy:
@@ -503,7 +503,7 @@ def test_sample_pymc3():
     if not pymc3:
         skip('PyMC3 is not installed. Abort tests for _sample_pymc3.')
     else:
-        with ignore_warnings(UserWarning):
+        with ignore_warnings(UserWarning): ### TODO: Restore tests once warnings are removed
             for X in distribs_pymc3:
                 samps = next(sample(X, size=size, library='pymc3'))
                 for sam in samps:
diff --git a/sympy/stats/tests/test_mix.py b/sympy/stats/tests/test_mix.py
index 9cbd336646..a47f7beae8 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, DiracDelta, Dummy, Sum, oo)
 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')
@@ -56,9 +57,13 @@ def test_compound_distribution():
 
 def test_mix_expression():
     Y, E = Poisson('Y', 1), Exponential('E', 1)
+    k = Dummy('k')
+    expr1 = Integral(Sum(exp(-1)*Integral(exp(-k)*DiracDelta(k - 2), (k, 0, oo)
+    )/factorial(k), (k, 0, oo)), (k, -oo, 0))
+    expr2 = Integral(Sum(exp(-1)*Integral(exp(-k)*DiracDelta(k - 2), (k, 0, oo)
+    )/factorial(k), (k, 0, oo)), (k, 0, oo))
     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): ### TODO: Restore tests once warnings are removed
+        assert P(E + Y < 2, evaluate=False).rewrite(Integral).dummy_eq(expr1)
+        assert P(E + Y > 2, evaluate=False).rewrite(Integral).dummy_eq(expr2)
diff --git a/sympy/stats/tests/test_rv.py b/sympy/stats/tests/test_rv.py
index 2c50fdedc9..26071dbc6f 100644
--- a/sympy/stats/tests/test_rv.py
+++ b/sympy/stats/tests/test_rv.py
@@ -191,7 +191,7 @@ def test_Sample():
     scipy = import_module('scipy')
     if not scipy:
         skip('Scipy is not installed. Abort tests')
-    with ignore_warnings(UserWarning):
+    with ignore_warnings(UserWarning): ### TODO: Restore tests once warnings are removed
         assert next(sample(X)) in [1, 2, 3, 4, 5, 6]
         assert isinstance(next(sample(X + Y)), float)
 
diff --git a/sympy/stats/tests/test_stochastic_process.py b/sympy/stats/tests/test_stochastic_process.py
index f4ff19a0e1..a3c281c92c 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)