Distribution of the ratio of Poisson $\lambda$'s in closed form

Being the evaluation of ratio of rates (to which the ratio of $\lambda$'s is related) an important issue in Physics, it is worth trying to get an analytic expression for its pdf. This can be done extending to the continuum Eq. (),¹³that is replacing the sums by integrals, and applying the constraint between the two variables by a Dirac delta [13]:

$$f(\rho_\lambda\,\vert\,x_1,x_2) = \int_0^\infty\!\!\!\int_0^\infty \!\delta\left(\rho_\lambda-\frac... ...\lambda_2}\right)\cdot f(\lambda_1\,\vert\,x_1)\cdot f(\lambda_2\,\vert\,x_2)\, \lambda_1 \lambda_2\,.$$ (9)

Making use of the properties of the $\delta()$, we can rewrite it as

$$\delta\left(\rho_\lambda-\frac{\lambda_1}{\lambda_2}\right) = \frac{\delta(\lambda_1-\lambda_1^*)} {\left\vert\frac{\mbox{d}}{\... ...lambda-\frac{\lambda_1}{\lambda_2}\right) \right\vert _{\lambda_1=\lambda_1^*}}$$ (10)

$$= \lambda_2\cdot \delta(\lambda_1-\lambda_1^*)\,,$$ (11)

with $\lambda_1^*$ root of the equation $\rho_\lambda-\lambda_1/\lambda_2=0$, and therefore equal to $\rho_\lambda\cdot\lambda_2$. Equation () becomes then

$$f(\rho_\lambda\,\vert\,x_1,x_2) = \int_0^\infty\!\!\!\int_0^\infty \!\lambda_2\cdot \delta(\lambda_... ...a\cdot\lambda_2) \cdot f(\lambda_1\,\vert\,x_1)\cdot f(\lambda_2\,\vert\,x_2)\, \lambda_1 \lambda_2$$ (12)

$$= \int_0^\infty\!\lambda_2\cdot \frac{(\rho_\lambda\cdot \lambda_2)... ...cdot\lambda_2}} {x_1!}\cdot \frac{\lambda_2^{x_2}\cdot e^{-\lambda_2}}{x_2!} \, \lambda_2$$ (13)

$$= \frac{\rho_\lambda^{x_1}}{x_1!\,x_2!}\cdot \int_0^\infty\!\lambda_2^{x_1+x_2+1}\cdot e^{-(\rho_\lambda+1)\cdot\lambda_2} \, \lambda_2\,.$$ (14)

Once more we recognize in the integrand something related to the Gamma distribution. In fact, identifying the power of $\lambda_2$ with `$\alpha-1$' of a Gamma pdf, and ` $(1+\rho_\lambda)$' at the exponent with the `rate parameter' $\beta$, that is

$$\alpha-1 = x_1+x_2+1$$ (15)

$$\beta = \rho_\lambda+1\,,$$ (16)

the integrand in Eq. () can be rewritten as

$$\lambda_2^{\alpha-1}\cdot e^{-\beta\,\lambda_2} = \frac{\Gamma(\alpha)}{\beta^{\,\alpha}}\cdot \left(\frac{\beta^{\... ...}}{\Gamma(\alpha)}\cdot \lambda_2^{\alpha-1}\cdot e^{-\beta\,\lambda_2} \right)$$ (17)

in order to recognize within parentheses a Gamma pdf in the variable $\lambda_2$, whose integral over $\lambda_2$ is then equal to one because of normalization. We get then

$$f(\rho_\lambda\,\vert\,x_1,x_2) = \frac{\rho_\lambda^{\,x_1}}{x_1!\,x_2!}\cdot \frac{\Gamma(\alpha)... ...,\alpha}}{\Gamma(\alpha)}\cdot \lambda_2^{\alpha-1}\cdot e^{-\beta\lambda_2} \, \lambda_2$$ (18)

$$= \frac{\rho_\lambda^{x_1}}{x_1!\,x_2!}\cdot \frac{\Gamma(\alpha)}{\beta^{\,\alpha}}$$ (19)

$$= \frac{\Gamma(x_1+x_2+2)}{\Gamma(x_1+1)\,\Gamma(x_2+1)}\cdot \frac{\rho_\lambda^{\,x_1}}{(\rho_\lambda+1)^{\,x_1+x_2+2}}$$ (20)

$$= \frac{(x_1+x_2+1)!}{x_1!\,x_2!}\cdot \rho_\lambda^{\,x_1}\cdot (\rho_\lambda+1)^{\,-(x_1+x_2+2)}$$ (21)

The mode of the distribution can be easily obtained finding the maximum (of the log) of the pdf, thus getting

$$(\rho_\lambda) = \frac{x_1}{x_2+2}\,,$$ (22)

in agreement with what we have got in Figs.  and by Monte Carlo (indeed, done there in a fast and rather rough way - see Appendix B.2).

In order to get expected value and standard deviation, we need to evaluate the relevant integrals¹⁴

• First we check that $f(\rho_\lambda\,\vert\,x_1,x_2)$ is properly normalized. Indeed the integral $\int_0^\infty\! f(\rho_\lambda\,\vert\,x_1,x_2)\,$d$\rho_\lambda$ is equal to unity for `all possible' $x_1$ and $x_2$.¹⁵
• The expected value is equal to
  

  $$(\rho_\lambda\,\vert\,x_1,x_2) = \frac{x_1+1}{x_2} \hspace{0.7cm}(\mathbf{x_2>0})\,,$$ (23)

  
  in perfect agreement with what we can read from the Monte Carlo results of Figs.  and .
• The expected value of $\rho_\lambda^2$ is given by
  

  $$(\rho_\lambda^2\,\vert\,x_1,x_2) = \frac{(x_1+1)\cdot (x_1+2)}{x_2\cdot (x_2-1)} \hspace{3.0cm}(\mathbf{x_2>1})\,,$$ (24)

  
  from which we evaluate (subtracting to it the square of the expected value)
  

  $$(\rho_\lambda\,\vert\,x_1,x_2) = \frac{x_1+1}{x_2}\cdot \left(\frac{x_1+2}{x_2-1} - \frac{x_1+1}{x_2}\right) \hspace{0.8cm}(\mathbf{x_2>1})\,,$$ (25)

  
  from which the standard deviation follows, that we rewrite in a more compact form as
  

  $$\sigma(\rho_\lambda) = \sqrt{\mu_{\rho_\lambda}\cdot \left(\frac{x_1+2}{x_2-1} - \mu_{\rho_\lambda}\right)} \hspace{4.0cm}(\mathbf{x_2>1})\,,$$ (26)

  
  having indicated by $\mu_{\rho_\lambda}$ the expected value of $\rho_\lambda$. For the values $x_1$ and $x_2$ used in Figs.  and , we get, starting from $\mathbf{x_1=x_2=2}$ in increasing order, the following standard deviations: 1.936, 1.247, 0.507, 0.269 and 0.143, in agreement with the Monte Carlo results (or the other way around).

The detailed comparison between closed expression of the pdf and the Monte Carlo outcome is shown in Fig.  for the toughest case we have met, that is $x_1=x_2=1$.

Figure: Comparison of the distribution of $\rho_\lambda=\lambda_1/\lambda_2$ obtained by the closed expression () with that estimated by Monte Carlo (same as top plot of Fig. ). The vertical lines indicate mode and expected value, evaluated using Eqs. () and (), equal to 1/3 and 2, respectively. (Note that none of these values is close to 1, that is what one would naively expect for the ratio of the rates - indeed, only for $x_1$ and $x_2$ above ${\cal O}(100)$ mode, expected value and ratio of the observed counts become approximately equal).