Bivariate Poisson

[s-morten]

Hello, I would like to incorporate a bivariate Poisson distribution into my model.
I couldn’t find the possibility for bivariate Poisson in Pymc, so I need to model it myself.
According to Wikipedia:

A simple way to generate a bivariate Poisson distribution X_{1},X_{2} is to take three independent Poisson distributions Y_1 , Y_2 , Y_3 with means λ_1 , λ_2 , λ_3 and then set X_1 = Y_1 + Y_3 , X_2 = Y_2 + Y_3.

So far, so good. Sounds simple enough to implement. The route I wanted to take:

with pymc.Model() as bivariate_poisson:
    y1 = pymc.Poisson("Y1", θ_1)
    y2 = pymc.Poisson("Y2", θ_2)
    y3 = pymc.Poisson("Y3", θ_3)

    x1 = pymc.Deterministic("X1", y1 + y3)
    x2 = pymc.Deterministic("X5", y2 + y3)

Edit: By trying to implement, I discovered that Deterministic can’t be an observed variable. Therefore, it probably also needs to be Poisson, but I am even more confused.

Is this correct? Firstly, is Deterministic correct to use? The Poisson distribution would result from the combination of the other Poissons, correct?

And the second question:
I logically only have two Poisson distributions, and therefore mu. Because of the bivariate Poisson I need three mu values.

with pymc.Model() as bivariate_poisson:
    y1 = pymc.Poisson("Y1", θ_1)
    y2 = pymc.Poisson("Y2", θ_2)
    y3 = pymc.Poisson("Y3", ???)

How do I choose the mu for the third independent Poisson? Do I just use a Normal distribution with low prior information?

Thank you for your help and your answers

[ckrapu]

Could you talk a little more about why you’d want to have correlated Poisson variates? Is there a specific reason you wouldn’t want to model (\log \theta_1, \log \theta_2) as draws from a bivariate normal? This construction is used quite a bit more frequently.

[s-morten]

Hello, thank you for your answer.
You are right, this is probably exactly what I want.

Could you maybe guide me with implementing the bivariate normal, as I don’t find it to be easy :D.
As I mentioned previously, right now, I have two independent Poisson.
\theta_1 = exp(N_1 * data_1 + N_2 * data_2)
\theta_2 = exp(N_3 * data_3 + N_4 * data_4)

Is it possible to use the results of the \theta 's as the mu of the bivariate normal(with or without the exp), because I don’t seem to be able to do so. And how to choose the parameter chol?

Thank you very much for your help

[ckrapu]

Here’s an example notebook with a few ways to implement the correlated Poisson model. There’s more discussion on this thread which is essentially the same topic:

Dependent histogram count data

Hey! This is perhaps a broad question not limited specifically to PyMC, but I figured this is a good place to ask My input data is standard (as in, some vector of of covariates for each sample), but my observations are less so. In fact, my observations are a matrix (a vector for each sample) with histogram count data. Specifically, for each sample, the target is a vector of counts for each bin, which are of course somewhat correlated: index | X | bin_0 | bin_1 | bin_2 | ... | …

There’s a lot to unpack in that example notebook - feel free to call out specific parts that you’d like explained.

[s-morten]

Thank you for your help!
That is a great resource that I did not find before.
Much appreciated