Predicting numbers of counts, their difference and their ratio

The Poisson distribution hardly needs any introduction, beside, perhaps, that it can be framed within the Poisson process, which has indeed its roots in the Bernoulli process. This picture makes the Poissonian related to other important distributions, as reminded in Appendix A, which can be seen as a technical preface to the paper.

Using the notation introduced there, the Poisson probability function is given by²

$$f(x\,\vert\,\lambda) \equiv P(X\!\!=\!x\,\vert\,\lambda) = \frac{\lambda^x}{x!}\cdot e^{-\lambda} \hspace{1.0 cm} \left\{ \b... ...}{l} 0 < \lambda < \infty \\ x = 0, 1, \ldots, \infty\\ \end{array} \right.\,$$ (1)

and it quantifies how much we believe that $x$ counts will occur, if we assume an exact value of the parameter $\lambda$.³As well known, expected value and standard deviation of $X$ are $\lambda$ and $\sqrt{\lambda}$. The most probable value of $X$ (`mode') is equal to the integer just below $\lambda$ (`floor$(\lambda)$') in the case $\lambda$ is not integer. Otherwise is equal to $\lambda$ itself, and also to $\lambda-1$ (remember that $\lambda$ cannot be null).

If we have two independent Poisson distributions characterized by $\lambda_1$ and $\lambda_2$, i.e.

$$X_1 \sim {\cal P}_{\lambda_1}$$

$$X_2 \sim {\cal P}_{\lambda_2}\,,$$

we `expect' their difference $D\!=\! X_1\!-\!X_2$ `to be' $(\lambda_1\!-\!\lambda_2) \pm \sqrt{\lambda_1+\lambda_2}$, as it results from well known theorems of probability theory⁴(hereafter, unless indicated otherwise, the notation ` $xxx \pm yyy$' stands for `expected value of the quantity $\pm$ its standard uncertainty' [9], that is the standard deviation of the associated probability distribution).

The probability distribution of $D$ can be obtained `from the inventory of the values of $X_1$ and $X_2$ that result in each possible value of $D$', that is

$$P(D\!=\!d\,\vert\,\lambda_1,\lambda_2) \equiv f(d\,\vert\,\lambda_1,\lambda_2) = \sum_{\begin{array}{c}x_1,x_2 \\ x_1\!-\!x_2\!=\!d\end{array}} \!\! f(x_1\,\vert\,\lambda_1)\cdot f(x_2\,\vert\,\lambda_2)\,.$$ (2)

For example, in the case of $\lambda_1=\lambda_2=1$, the most probable contributions to $D$ are shown in Tab. .

Table: Table of the most probable differences $D\!=\! X_1\!-\!X_2$ for $\lambda_1\!=\!\lambda_2\!=\!1$ (probability of each entry in the table within square brackets).

		$X_2$
		0	1	2	3	4	5
	0	0	-1	-2	-3	-4	-5
		[0.135335]	[0.135335]	[0.067668]	[0.022556]	[0.005639]	[0.001128]
	1	1	0	-1	-2	-3	-4
		[0.135335]	[0.135335]	[0.067668]	[0.022556]	[0.005639]	[0.001128]
	2	2	1	0	-1	-2	-3
$X_1$		[0.067668]	[0.067668]	[0.033834]	[0.011278]	[0.002819]	[0.000564]
	3	3	2	1	0	-1	-2
		[0.022556]	[0.022556]	[0.011278]	[0.003759]	[0.000940]	[0.000188]
	4	4	3	2	1	0	-1
		[0.005639]	[ 0.005639]	[0.002819]	[0.000940]	[0.000235]	[0.000047]
	5	5	4	3	2	1	0
		[0.001128]	[0.001128]	[0.000564]	[0.000188]	[0.000047]	[0.000009]

For instance, the probability to get $D=0$ sums up to 30.9%. The probability decreases symmetrically for larger absolute values of the difference. Without entering into the question of getting a closed form of $f(d\,\vert\,\lambda_1,\lambda_2)$,⁵ it can be instructive to implement Eq. (), although in an approximate and rather inefficient way, in a few lines of R code [11]:⁶

 
dPoisDiff <- function(d, lambda1, lambda2) {
   xmax = round(max(lambda1,lambda2)) + 20*sqrt(max(lambda1,lambda2))
   sum( dpois((0+d):xmax, lambda1) * dpois(0:(xmax-d), lambda2) )
}

This function is part of the code provided in Appendix B.1, which produces the plot of Fig. , evaluating also expected value and standard deviation (indeed approximated values, being xmax not too large).

Figure: Distribution of the difference of counts resulting from two Poisson distributions with $\lambda_1=\lambda_2=1$.

Moving to the ratio of counts, numerical problems might arise, as shown in Tab. , analogue of Tab. .

Table: Table of the most probable ratios $X_1/X_2$ for $\lambda_1\!=\!\lambda_2\!=\!1$. `NaN' and `Inf' are the R symbols for undefined (`not a number') and infinity, resulting from a vanishing denominator.

		$X_2$
		0	1	2	3	4	5
	0	red NaN	0	0	0	0	0
		[ red0.135335]	[0.135335]	[0.067668]	[0.022556]	[0.005639]	[0.001128]
	1	red Inf	1	1/2	1/3	1/4	1/5
		[red 0.135335]	[0.135335]	[0.067668]	[0.022556]	[0.005639]	[0.001128]
	2	red Inf	2	1	2/3	1/2	2/5
$X_1$		[red 0.067668]	[0.067668]	[0.033834]	[0.011278]	[0.002819]	[0.000564]
	3	red Inf	3	3/2	1	3/4	3/5
		[red 0.022556]	[0.022556]	[0.011278]	[0.003759]	[0.000940]	[0.000188]
	4	red Inf	4	2	4/3	1	4/5
		[red 0.005639]	[ 0.005639]	[0.002819]	[0.000940]	[0.000235]	[0.000047]
	5	red Inf	5	5/2	5/3	5/4	1
		[red 0.001128]	[0.001128]	[0.000564]	[0.000188]	[0.000047]	[0.000009]

In fact for rather small values of $\lambda_2$ there is high chance (exactly the probability of getting $X_2=0$) that the ratio results in an undefined form or an infinite, reported in the table using the R symbols NaN (`not a number') and Inf, respectively. As we can see, we have now quite a variety of possibilities and the probability distribution of the ratios is rather irregular. For this reason, in this case we evaluate it by Monte Carlo methods using R.⁷Figure 

Figure: Monte Carlo distribution of the ratio of counts resulting from two Poisson distributions with $\lambda_1=\lambda_2$.

shows the distributions of the ratio for $\lambda_1\!=\!\lambda_2 = 1,2,3$. The figure also reports the probability to get an infinite or an undefined expression, equal to $P(X_2\!=\!0\,\vert\,\lambda_i)$. When $\lambda_2$ is very large the probability to get $X_2=0$, and therefore of $X_1/X_2$ being equal to Inf or NaN, vanishes. But the distribution of the ratio remains quite `irregular', if looked into detail, even for $\lambda_1$ `not so small', as shown in Fig.  for the cases of $\lambda_1\!=\!5,10,20,50$ and $\lambda_2\!=\!1000$.

Figure: Monte Carlo distribution of the ratio of counts resulting from two Poisson distributions with $\lambda_1=5,10,20,50$ and $\lambda_2=1000$.

However we should not be worried about this kind of distributions, which are not more than entertaining curiosities, as long as physics questions are concerned. Why should we be interested in the ratio of counts that we might observe for different $\lambda$'s? If we want to get an idea of how much the counts could differ, we can just use the probability distribution of their possible differences, which has a regular behavior, without divergences or undefined forms.

The deep reason for speculating about “ratios of small numbers of events” and their “errors”[3] is due to a curious ideology at the basis of a school of Statistics which limits the applications of Probability Theory. Indeed, we, as physicists, are often interested in the ratio of the rates of Poisson processes, that is in $\rho=r_1/r_2$, being this quantity related to some physical quantities having a theoretical relevance. Therefore we aim to learn which values of $\rho$ are more or less probable in the light of the experimental observations. Stated in this terms, we are interested in evaluating `somehow' (not always in closed form) the probability density function (pdf) $f(\rho\,\vert\,x_1,T_1,x_2,T_2,I)$, given the observations of $x_1$ counts during $T_1$ and of $x_2$ counts during $T_2$ (and also conditioned on the background state of knowledge, generically indicated by $I$). But there are statisticians who maintain that we can only talk about the probability of $X$ counts, assuming $\lambda$, and not of the probability distribution of $\lambda$ having observed $X\!=\!x$, and even less of $\lambda_1/\lambda_2$ (same as $r_1/r_2$, if $T_1\!=\!T_2$) having observed $x_1$ and $x_2$ counts.⁸

If, instead, we follow an approach closer to that innate in the human mind, which naturally attaches the concept of probability to whatever is considered uncertain [16], there is no need to go through contorted reasoning. For example, if we observe $X=1$, then we tend to believe that this effect has been caused more likely by a value of $\lambda$ around 1 rather than around 10 or 20, or larger, although we cannot rule out with certainty such values. Similarly, sticking to the observation of $X=1$, we tend to believe to $\lambda\approx 1$ much more than to $\lambda \approx 10^{-2}$, or smaller. In particular, $\lambda=0$ is definitely ruled out, because it could only yield 0 counts (this is just a limit case, since the Poisson $\lambda$ is defined positive).

It is then clear that, as far as the ratio of $\lambda_1/\lambda_2$ is concerned, there are no divergences, no matter how small the numbers of counts might be, obviously with two exceptions. The first is when $X_2$ is observed to be exactly 0 (but in this case we could turn our interest in $\lambda_2/\lambda_1$, assuming $X_1>0$). The second is when $X_2$ is not zero, but there could be some background, such that $r_2=0$ is not excluded with certainty [13,17]. The effect of possible background is not going to be treated in detail in this paper, and only some hints on how to include it into the model will be given.