# G test with small observed counts (~zero)

I have an experiments that counts how many events happened during a certain period of time. The experiments is repeated many times in different "runs" (I have 65 runs).

Here for example the list of observed events for every run:

obs = [48, 88, 9, 38, 84, 84, 8, 14, 47, 6, 202, 156, 69, 460, 75, 145, 243, 390, 3, 541, 613, 57, 353, 231, 5, 360, 484, 705, 745, 61, 516, 437, 71, 96, 183, 1121, 405, 66, 1000, 995, 271, 93, 368, 10, 182, 509, 1044, 0, 208, 732, 23, 817, 415, 238, 111, 1765, 188, 1129, 283, 89, 1550, 1974, 471, 281, 825]


and here the duration of every run (the number of events is proportional to the time of observation):

times = [6.25, 11.91, 0.66, 4.22, 9.00, 10.87, 0.68, 1.60, 6.65, 0.74, 21.82, 19.96, 9.96, 55.58, 7.32, 17.10, 32.03, 52.81, 0.35, 68.68, 75.22, 6.77, 48.46, 30.07, 1.02, 45.11, 66.18, 89.72, 95.19, 8.43, 69.69, 59.28, 7.45, 12.42, 24.83, 151.14, 45.64, 7.16, 132.74, 126.59, 35.19, 12.45, 46.30, 0.85, 24.19, 64.07, 136.96, 10.55, 26.52, 93.41, 2.97, 108.62, 57.15, 29.44, 13.52, 235.91, 24.08, 146.92, 39.54, 10.21, 202.20, 249.35, 57.50, 35.37, 104.50]


There is no constraint on the total number of events. There is no correlation between runs.

I want to test the null hypothesis "the rate is the same for every run (two side test: the incompatibility could be for excess or deficit)".

From the data I now:

mean_obs_rate = sum(obs) / sum(times) = 0.0352


So I can compute the expected number of events under the null hypothesis, using the observed mean rate exp = mean_obs_rate * times:

exp = [48.31, 92.00, 5.12, 32.61, 69.55, 83.99, 5.22, 12.33, 51.35, 5.74, 168.58, 154.16, 76.92, 429.37, 56.57, 132.13, 247.46, 407.94, 2.69, 530.54, 581.06, 52.33, 374.35, 232.29, 7.90, 348.46, 511.26, 693.06, 735.33, 65.08, 538.34, 457.92, 57.55, 95.91, 191.79, 1167.55, 352.56, 55.28, 1025.39, 977.87, 271.81, 96.17, 357.68, 6.60, 186.90, 494.95, 1058.01, 81.54, 204.90, 721.59, 22.92, 839.10, 441.52, 227.41, 104.41, 1822.38, 186.02, 1134.93, 305.44, 78.85, 1562.01, 1926.26, 444.23, 273.23, 807.27]


Then I have tried with the likelihood ratio using as a model a product of Poissonian distribution:

$$-2 \log \frac{\sup_\lambda \text{Pois}(n_1|\lambda t_1)\cdots\text{Pois}(n_{65}|\lambda t_{65})}{\sup_{\{\lambda_i\}} \text{Pois}(n_1|\lambda_1 t_1)\cdots\text{Pois}(n_{65}|\lambda_{65}t_{65})}$$

The two hypotheses are nested (you can write $\lambda_i = \lambda + \delta\lambda_i$ and the numerator is the special case with $\delta\lambda_i=0$)

The denominator is maximized with $\lambda_i t_i = n_i$, while the numerator is maximized with $\lambda t_i = t_i\times$ mean_obs_rate.

After some simplifications I have found that this is equal to:

$$2 \sum_i n_i \log\frac{n_i}{e_i}$$ with $e_i = \lambda t_i$. This is exactly the G-test. That should be asymptotically distributed as a $\chi^2$ distribution with 65 - 1 dof (even if generating pseudo experiments under the null-hypothesis I found that there is a small bias). Small number of the G-test should be good compatibility, while large should be bad compatibility with the null-hypothesis. Is this correct, taking into account the fact that I want a two-sided test?

I was trying to see what is the contribution to the G-test of every run ($n_i \log\frac{n_i}{e_i}$), and I looked to one where the number of observed event is 0 and expected is 81.54. It should give a large contribution to the G-test quantity, since its contribution should be proportional to the probability to observe 0 events when the expected is 81.45, but actually its contribution is 0 (just taking the limit for $n_i\to 0$ in the $G$ definition). That is completely the opposite compared to a simpler $\chi^2$-test $(n_i - e_i)^2 / e_i = (0 - 81.54)^2 / 81.54 = 81.54$.

So the question: is this ok? How should I implement runs where I observe 0 events? Why runs with 0 observed events count nothing in the final value of the G-test? Should they contribute in the number of degree of freedom of the asymptotic $\chi^2$ distribution?

The likelihood ratio

$$-2\log\Lambda = -2 \log \frac{\sup_\lambda \text{Pois}(n_1|\lambda t_1)\cdots\text{Pois}(n_{65}|\lambda t_{65})}{\sup_{\{\lambda_i\}} \text{Pois}(n_1|\lambda_1 t_1)\cdots\text{Pois}(n_{65}|\lambda_{65}t_{65})}$$

can be simplified to:

$$-2\log\Lambda = 2\sum (e_i - n_i + n_i\log(n_i/e_i))$$

then using $\sum e_i=\sum n_i$, that hold in this case, the expression can be further simplified to:

$$-2\log\Lambda = 2\sum n_i\log(n_i/e_i)$$

The point is that the last simplification holds when doing the sum, not when considering each single term of the summation. So, if I want to see the contribution of each term to the G-test I have to consider the quantity $2 (e_i - n_i + n_i\log(n_i/e_i))$, which has a minimum for $e_i=n_i$.

In the figure the dotted line is the expression above, while the dotted lines are the chi-square $(e_i - n_i)^2 / e_i$ for different values of $n$ ($O$ in the figure) as a function of $e_i$ ($E$ in the figure).