Implement RV convolutions

[ricardoV94]

We could try to wrap some method from scipy.integrate and use it to convolve continuous variables in a joint graph of z = x + y where neither variable {x, y} is being conditioned on.

I don't know how feasible this would be, since scipy requires a callable. Can we use something like OpFromGraph?

Obviously we can (and should) implement the optimized graphs for variables that have simple analytical form such as the convolution of two gaussians (which I think pymc-symbolic could already do)?

[brandonwillard]

Are you talking about creating a numerical integration Op? We could definitely do that, but that approach falls well outside of the symbolic stuff we're doing here. Plus, that approach has severe limitations (e.g. the kind that often lead one to MCMC), so it's difficult to imagine such a thing being applied automatically.

We could reproduce some simple symbolic convolutions and build on those over time. There's a lot that can be done with only Laplace and Mellin transforms in the area of probability theory, so we wouldn't really need to implement symbolic integration in its more general forms.

[rlouf]

We should start with reasoning at the level of random variables. We can already implement relations for convolutions of known distributions that result in a known distribution in AeMCMC, see aesara-devs/aemcmc#57.

Efforts in this direction should be pursued following the same approach.