Repository: sympy/sympy

Issue:
Sum.doit() gives incorrect results with Random Indexed Variables
While working on #19428, I came across the following error
```python
>>> from sympy.stats import *
>>> from sympy import symbols, S, Eq, Sum
>>> B = BernoulliProcess('B', S(1)/2, 1, -1)
>>> t = symbols('t', positive=True, integer=True)
>>> expr = Sum(B[t], (t,0,5))
>>> expr
Sum(B[t], (t, 0, 5))
>>> expr.doit()  # this should not consider B[t] equal in all the cases as it is random variable
6*B[t]          # incorrect result
```
Expected behaviour, it should produce: `B[0] + B[1] + B[2] + B[3] + B[4] + B[5]`

Discussion:

Patch:
diff --git a/sympy/stats/rv.py b/sympy/stats/rv.py
index aa81d18908..4369832cae 100644
--- a/sympy/stats/rv.py
+++ b/sympy/stats/rv.py
@@ -310,6 +310,14 @@ def key(self):
         elif isinstance(self.symbol, Function):
             return self.symbol.args[0]
 
+    @property
+    def free_symbols(self):
+        if self.key.free_symbols:
+            free_syms = self.key.free_symbols
+            free_syms.add(self)
+            return free_syms
+        return {self}
+
 class RandomMatrixSymbol(MatrixSymbol):
     def __new__(cls, symbol, n, m, pspace=None):
         n, m = _sympify(n), _sympify(m)


Test patch:
diff --git a/sympy/stats/tests/test_stochastic_process.py b/sympy/stats/tests/test_stochastic_process.py
index 4157a4f50b..ae8800bc0a 100644
--- a/sympy/stats/tests/test_stochastic_process.py
+++ b/sympy/stats/tests/test_stochastic_process.py
@@ -1,5 +1,5 @@
 from sympy import (S, symbols, FiniteSet, Eq, Matrix, MatrixSymbol, Float, And,
-                   ImmutableMatrix, Ne, Lt, Gt, exp, Not, Rational, Lambda,
+                   ImmutableMatrix, Ne, Lt, Gt, exp, Not, Rational, Lambda, Sum,
                    Piecewise)
 from sympy.stats import (DiscreteMarkovChain, P, TransitionMatrixOf, E,
                          StochasticStateSpaceOf, variance, ContinuousMarkovChain,
@@ -151,7 +151,7 @@ def test_BernoulliProcess():
     H, T = symbols("H,T")
     assert E(X[1]+X[2]*X[3]) == H**2/9 + 4*H*T/9 + H/3 + 4*T**2/9 + 2*T/3
 
-    t = symbols('t', positive=True, integer=True)
+    t, x = symbols('t, x', positive=True, integer=True)
     assert isinstance(B[t], RandomIndexedSymbol)
 
     raises(ValueError, lambda: BernoulliProcess("X", p=1.1, success=1, failure=0))
@@ -194,3 +194,14 @@ def test_BernoulliProcess():
     assert P(B[5] > 0, B[5]) == BernoulliDistribution(0.6, 0, 1)
     raises(ValueError, lambda: P(3))
     raises(ValueError, lambda: P(B[3] > 0, 3))
+
+    # test issue 19456
+    expr = Sum(B[t], (t, 0, 4))
+    expr2 = Sum(B[t], (t, 1, 3))
+    expr3 = Sum(B[t]**2, (t, 1, 3))
+    assert expr.doit() == B[0] + B[1] + B[2] + B[3] + B[4]
+    assert expr2.doit() == Y
+    assert expr3.doit() == B[1]**2 + B[2]**2 + B[3]**2
+    assert B[2*t].free_symbols == {B[2*t], t}
+    assert B[4].free_symbols == {B[4]}
+    assert B[x*t].free_symbols == {B[x*t], x, t}