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.