I have an additional question on this topic:
How would you code the above covariances if the rate of each Poisson depended on say two parameters? i.e. log_rates = constant + coeff * data
Where we can assume constant and coeff are drawn from something like:
constant = pm.Normal('const', mu=0., sigma=1.)
coeff = pm.Normal('b0', mu=0., sigma=1.)