Finding the joint distribution from Poisson conditionals Suppose that for two discrete random variables $X_1$ and $X_2$, we know their conditional distributions. Namely
$$X_1~|~X_2 = x_2 \sim \mathrm{Poisson}(\lambda_1 + ax_2),$$
$$X_2~|~X_1 = x_1 \sim \mathrm{Poisson}(\lambda_2 + bx_1).$$
We want to calculate their joint distribution, $p(X_1, X_2)$.
My own idea is to divide the equations above and find $\frac{p_{X_1}(x_1)}{p_{X_2}(x_2)}$ and then sum over all values of $X_1$ to find $1/p_{X_2}(x_2)$. But for this, I have to compute a very bad series. 
Do you have any idea for this problem?
P.S. The series I have to compute is of form
$$\sum_{n=0}^\infty \frac{c^n}{n!(a+nb)^k},$$
which I hardly believe to have a closed form.
 A: Unless I've made a mistake, this doesn't seem to be possible.
One thing you can do is look at the joint MGF.  Conditioning on $X_1$,
$$ \eqalign{M(s,t) &= \mathbb E\left[ \exp(s X_1 + t X_2) \right] = 
    \mathbb E \left[ \exp(s X_1) \mathbb E[\exp(t X_2) | X_1] \right] \cr
&= \mathbb E \left[ \exp(s X_1) \exp((\lambda_2 + b X_1)(e^t-1))\right]\cr
&= \exp(-\lambda_2 (e^t-1)) M(s+ b (e^t-1),0)}$$
Similarly, conditioning on $X_2$,
$$M(s,t) =\exp(-\lambda_1 (e^s-1)) M(0, t + a (e^s-1)) $$
Combining these and replacing $s$ by $s - b (e^t-1)$:
$$ M(s,0) = \exp(\lambda_1 (1 - e^{s + b (1-e^t)})) + \lambda_2(e^t-1)) M(0, a\exp(s + b(1-e^t))-a+t)$$
and similarly
$$ M(0,t) =  \exp(\lambda_2 (1 - e^{t + a (1-e^s)})) + \lambda_1 (e^s-1)) M(b\exp(t + a(1-e^s))-b+s, 0) $$
But replacing $t$ by $a\exp(s + b(1-e^t))-a+t$ here gives us a rather complicated functional equation for $M(s,0)$: 
$$ M \left( s,0 \right) =M \left( b{{\rm e}^{a{{\rm e}^{-b{{\rm e}^{t}}+b
+s}}-a{{\rm e}^{s}}+t}}-b+s,0 \right) {{\rm e}^{-\lambda_{{2}}{{\rm e}
^{a{{\rm e}^{-b{{\rm e}^{t}}+b+s}}-a{{\rm e}^{s}}+t}}-\lambda_{{1}}{
{\rm e}^{-b{{\rm e}^{t}}+b+s}}+\lambda_{{2}}{{\rm e}^{t}}+{{\rm e}^{s}
}\lambda_{{1}}}}
$$
Now taking the derivative with respect to $t$ and evaluating at $t=0$ gives us a differential equation for $M(s,0)$:
$$\dfrac{\partial M}{\partial s}(s,0) = -{\frac {M \left( s,0 \right)  \left( -\lambda_{{2}} \left( -ab{
{\rm e}^{s}}+1 \right) +\lambda_{{1}}{{\rm e}^{s}}b+\lambda_{{2}}
 \right) }{b \left( -ab{{\rm e}^{s}}+1 \right) }}
$$
which we can solve (with initial condition $M(0,0)=1$):
$$M(s,0) =  \left( \dfrac{abe^s-1}{ab -1}  \right)^{\frac{\lambda_2}{b} + \frac{\lambda_1}{ab}}
$$
Unfortunately, this does not satisfy the functional equation.  So it seems there is no solution in general.
