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