I spent a major part of this week working on product probability spaces, in #14730. This would allow SymPy to compute the probability/expectation of expressions with multiple random variables. For example:

```
>>> from sympy.stats import Poisson, P, Geometric
>>> from sympy import S
>>> X1, X2 = Geometric('X1', S(1)/2), Geometric('X2', S(1)/3)
>>> P(X1 + X2 < 5) #This was not possible for discrete random variables earlier.
```

While the PR gives the correct result, one thing that requires improvement is that the result is often a complicated expression. In most cases, it comprises of multiple `summations`

that cannot be simplified further by SymPy as of now. For instance, the result of the above is as follows.

```
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)))
```

Francesco pointed out that this calls for a change in API, where the `Probability`

object can be returned as per the convinience of the user, but this needs to undergo the deprecation process, which cannot be done unless a clear alternative is available.
Meanwhile, I started working on Joint probability spaces, which is what I plan on spending the next week on as well.