Add support for the Poisson-Gamma closed-form posterior

[rlouf]

Consider the following observation model:

import aesara
import aesara.tensor as at

srng = at.random.RandomStream(4238)

alpha = at.scalar()
beta = at.scalar()
observed = at.scalar()

Z_rv = srng.gamma(alpha, beta, name="Z_1")
Y_rv = srng.poisson(Z_rv, name="Y")

Z_rv is a conjugate prior of Y_rv and if the observed value of the latter is observed the model is equivalent to the following:

Z_rv = srng.gamma(alpha + observed, beta + 1)