Ratio of Poisson to sum of Poisson

Question

Let $X$ and $Y$ be independent random variables which follows Poisson distribution with rate $\lambda_1$ and $\lambda_2$ respectively. I would like to know and find (if it is possible) the distribution for $R=\frac{X}{X+Y}$. Does the distribution of $R$ exist and if yes, how to find it? A part of me suggests that the distribution does not exist for both $x,y=0$, simultaneously.

Answer 1

Special cases

Suppose $R$ is defined to be $r_0$ when $X=0$ and $Y=0$. That happens with probability $e^{-\lambda_1-\lambda_2}$.

The probability that $R=1$ can be found by knowing that $X>0$ and $Y=0$:

$$\text{Pr}(R=1)=\sum _{s=1}^{\infty } \frac{\exp (-\lambda_1) \exp (-\lambda_2) \lambda_2^s}{s!}=\left(e^{\lambda_1}-1\right) e^{-\lambda_1-\lambda_2}$$

(This, of course, is ignoring the possibility that $r_0$ has been defined to be 1.)

The probability that $R=1/2$ is found when $X=Y>0$:

$$\sum _{s=1}^{\infty } \frac{\exp (-\lambda_1) \exp (-\lambda_2) \lambda_1^s \lambda_2^s}{s! s!}=e^{-\lambda_1-\lambda_2} \left(I_0\left(2 \sqrt{\lambda_1} \sqrt{\lambda_2}\right)-1\right)$$

Where $I_0$ is the modified Bessel function of the first kind. (Here a closed-form solution might very well be in the eye of the beholder.)

A similar approach for $R=1/3$ can be constructed:

$$\sum _{s=1}^{\infty } \frac{\exp (-\lambda_1) \exp (-\lambda_2) \lambda_1^s \lambda_2^{2 s}}{s! (2 s)!}=e^{-\lambda_1-\lambda_2} \left(\, _0F_2\left(;\frac{1}{2},1;\frac{\lambda_1 \lambda_2^2}{4}\right)-1\right)$$

where $_0F_2$ is the generalized hypergeometric function.

General case

How one goes about this depends on what software is going to be used to implement specific values of $\lambda_1$ and $\lambda_2$. Here I'm using Mathematica. A brute-force formula for any specific value of $R$ (ignoring the value set for $r_0$) is as follows:

pr[r_, λ1_, λ2_, n_] := Sum[Sum[Boole[x1/s == r] Binomial[s, x1] (λ1/(λ1 + λ2))^x1*
  (1 - λ1/(λ1 + λ2))^(s - x1), {x1, 0, s}]*Exp[-(λ1 + λ2)] (λ1 + λ2)^s/s!, {s, 1, n}]

As a partial check consider values of $R$ being $0$, $1$, and $1/2$ from the formulas above:

$\text{Pr}(R=0)$

E^(-λ1 - λ2) (-1 + E^λ1) /. {λ1 -> 2, λ2 -> 5}
(* (-1 + E^5)/E^7 *)
pr[0, 2, 5, Infinity]
(* (-1 + E^5)/E^7 *)

$\text{Pr}(R=1)$

E^(-λ1 - λ2) (-1 + E^λ1)
(* (-1 + E^2)/E^7 *)
pr[1, 2, 5, Infinity]    
(* (-1 + E^2)/E^7 *)

$\text{Pr}(R=1/2)$

Sum[Exp[-λ1] λ1^s/s! Exp[-λ2] λ2^s/s!, {s, 1, Infinity}] /. {λ1 -> 2, λ2 -> 5} // N
(* 0.0815915 *)
pr[1/2, 2, 5, Infinity] // N
(* 0.0815915 *)

So how to determine a set of ratios such that most of the distribution is covered? We can determine a value for the total ($X+Y$) that is relatively extreme and figure out all of the unique values of $R$ given that threshold for the total and then calculate all of the probabilities. The quality of that job will be if the sum of the resulting probabilities is close to 1.

With[{λ1 = 2, λ2 = 5},
  nmax = Ceiling[λ1 + λ2 + 6 Sqrt[λ1 + λ2]];
  t = Flatten@Table[i/j, {j, 1, nmax}, {i, 0, j}] // DeleteDuplicates //
Sort;
  dist = {#, pr[#, λ1, λ2, nmax]} & /@ t;]
Total[dist[[All, 2]]] // N
(* 0.999088 *)

This approach covers 0.999088 of the total probability for this particular case.

For this example here is a plot of the values of $R$ and the associated probabilities: