The median of a jittered Poisson distribution

TL;DR

This paper analyzes the median of a jittered Poisson distribution, showing it approximates λ+1/3 for large λ, and introduces a simple, robust estimator for λ that is effective even with very large datasets.

Contribution

It establishes the asymptotic behavior of the median of a jittered Poisson distribution and proposes a new robust estimator for the Poisson parameter λ.

Findings

Median of N_λ + U is close to λ + 1/3 as λ→∞

Proposed estimator is consistent and asymptotically normal

Estimator is suitable for large datasets (up to 10^9 observations)

Abstract

Let $N_\lambda$ and $U$ be two independent random variables respectively distributed as a Poisson distribution with parameter $\lambda >0$ and a uniform distribution on $(0,1)$. This paper establishes that the median, say $M$, of $N_\lambda+U$ is close to $\lambda +1/3$ and more precisely that $M-\lambda-1/3=o(\lambda^{-1})$ as $\lambda\to \infty$. This result is used to construt a very simple robust estimator of $\lambda$ which is consistent and asymptotically normal. Compared to known robust estimates, this one can still be used with large datasets ($n\simeq 10^9$).

Tables

Table 1. Table 1: Average computational time in seconds (based on 10 replications) to evaluate each estimate based on a sample of size n 𝑛 n from a Poisson distribution with parameter λ = π 𝜆 𝜋 \lambda=\pi . The notation 0 ⋆ superscript 0 ⋆ 0^{\star} means that the average time is smaller than 0.05 sec.

	Sample size $n$
	${10}^4$	${10}^5$	${10}^6$	${10}^7$	${10}^8$	${10}^9$
${\hat{\lambda }}^{\mathrm{MLE}}$	$0^{\star }$	$0^{\star }$	$0^{\star }$	$0^{\star }$	0.1	1.2
$\widehat{\mathrm{Me}}\left({𝐍}_{\lambda }\right)$	$0^{\star }$	$0^{\star }$	$0^{\star }$	0.2	2.8	20.1
${\hat{\lambda }}^J$	$0^{\star }$	$0^{\star }$	0.1	0.5	5.5	48.2
${\hat{\lambda }}^{\mathrm{Tukey}}$	$0^{\star }$	0.1	0.8	8.2	98.2	NA

Equations

$$f_{Z_{\lambda }}(t)=\mathrm{P}(N_{\lambda }=\lfloor t\rfloor ),t\ge 0.$$

$$\sqrt{n}\left(\widehat{\mathrm{Me}}({\mathbf{Z}}_{\lambda })-{\mathrm{Me}}_{Z_{\lambda }}\right)\to N(0,{\sigma }_{\lambda }^2)$$

$$-\log 2\le {\mathrm{Me}}_{N_{\lambda }}-\lambda \le \frac{1}{3}.$$

$$\underset{\lambda \to \infty }{\mathrm{liminf}}{\mathrm{Me}}_{N_{\lambda }}-\lambda =-\frac{2}{3}\text{ and }\underset{\lambda \to \infty }{\mathrm{limsup}}{\mathrm{Me}}_{N_{\lambda }}-\lambda =\frac{1}{3}.$$

$$-\log (2)\le {\mathrm{Me}}_{Z_{\lambda }}-\lambda \le \frac{4}{3}.$$

$$\underset{\lambda \to \infty }{\mathrm{liminf}}\lambda {\delta }_{\lambda }=-\frac{8}{405}\text{ and }\underset{\lambda \to \infty }{\mathrm{limsup}}\lambda {\delta }_{\lambda }=\frac{4}{135}.$$

$$F_{Z_{\lambda }}(t)=\mathrm{P}(Z_{\lambda }\le t)=\sum _{k=0}^{\lfloor t\rfloor }\mathrm{P}(N_{\lambda }=k)+(t-\lfloor t\rfloor )\mathrm{P}(N_{\lambda }=\lfloor t\rfloor ),$$

$$M{e}_{Z_{\lambda }}=\lambda +\frac{1}{3}+\frac{\mathcal{H}(\lambda -\lfloor \lambda \rfloor )}{\lambda }+o\left(\frac{1}{\lambda }\right)$$

$$\mathcal{H}(x)=\left\{\begin{array}[]{ll}\frac{x^{2}(x-1)}{3}+\frac{4}{135},&\mbox{ if }x\in[0,2/3]\\ \frac{x}{3}(x^{2}-4x+5)-\frac{86}{135},&\text{ if }x\in[2/3,1].\end{array}\right.$$

$$\underset{x\in [0,1]}{\inf }\mathcal{H}(x)=-\frac{8}{405}\text{ and }\underset{x\in [0,1]}{\sup }\mathcal{H}(x)=\frac{4}{135}.$$

$$\frac{\mathcal{H}(x)-\epsilon }{n+x}\le {\mathrm{Me}}_{Z_{n+x}}-(n+x)-\frac{1}{3}\le \frac{\mathcal{H}(x)+\epsilon }{n+x}.$$

$$V_n\stackrel{d}{=}\sum _{i=1}^n(N_i-1)+\left(N_x+U-x-1/3-\frac{k}{n+x}\right)$$

$$w_n(x,k)\le w_{n+1}(x,k)\le \frac{1}{2}\iff {\mathrm{Me}}_{Z_{n+x}}\ge n+x+\frac{1}{3}+\frac{k}{n+x}$$

$$w_n(x,k)\ge w_{n+1}(x,k)\ge \frac{1}{2}\iff {\mathrm{Me}}_{Z_{n+x}}\le n+x+\frac{1}{3}+\frac{k}{n+x}.$$

$$w_n(x,k)=\mathrm{P}\left(N_{n+x}\le n\right)+\left(x-\frac{2}{3}+\frac{k}{n+x}\right)\mathrm{P}\left(N_{n+x}=n\right).$$

$$w_n(x,k)=\mathrm{P}\left(N_{n+x}\le n\right)+\left(x-\frac{2}{3}+\frac{k}{n+x}\right)\mathrm{P}(N_{n+x}=n+1).$$

$$\Delta_{n}(x,k)=\frac{3}{2(n+1+x)^{2}}\big{(}\mathcal{H}(x)-k\big{)}\;+\;o\left(\frac{1}{n^{2}}\right).$$

$${\hat{\lambda }}^{\mathrm{J}}-\lambda \approx N\left(0,\frac{{\sigma }_{\lambda }^2}{n}\right),$$

$${\Delta }_{\lambda }=2\sqrt{n}({\hat{\lambda }}^{\mathrm{J}}-\lambda )\mathrm{P}\left(N_{{\hat{\lambda }}^{\mathrm{J}}}=\lfloor {\hat{\lambda }}^{\mathrm{J}}+1/3\rfloor \right)$$

$${\psi }_{k,a}(y,\lambda )=\left(\frac{y-\lambda }{\sqrt{\lambda }}-a\right){\left(k^2-{\left(\frac{y-\lambda }{\sqrt{\lambda }}\right)}^2\right)}^2{\mathrm{I}}_{[-k,k]}\left(\frac{y-\lambda }{\sqrt{\lambda }}-a\right)$$

$${\tilde{N}}_{i,\lambda }=N_{i,\lambda }+(1-{\epsilon }_i)\sqrt{h},i=1,\dots ,n,$$

$$w_n(x,k)=\mathrm{P}\left(N_{n+x}\le n\right)+\left(x-\frac{2}{3}+\frac{k}{n+x}\right)\mathrm{P}\left(N_{n+x}=n\right).$$

$$w_n(x,k)=\mathrm{P}\left(N_{n+x}\le n\right)+\left(x-\frac{2}{3}+\frac{k}{n+x}\right)\mathrm{P}(N_{n+x}=n+1).$$

$$\Delta_{n}(x,k)=\frac{(n+1)!}{g_{n+1}(n+1+x)}\bigg{(}w_{n+1}(x,k)-w_{n}(x,k)\bigg{)}.$$

$${\Delta }_n(x,k)=$$

$$+r_{n+1}(x,k)-\frac{n+1}{n+x}{c}_n(0,x)r_n(x,k)$$

$$c_n(v,x)={\left(\frac{n+v+x}{n+1+x}\right)}^{n+1}\exp (1-v).$$

$${\Delta }_n(x,k)=$$

$$+\frac{n+1+x}{n+2}{r}_{n+1}(x,k)-c_n(0,x)r_n(x,k).$$

$$w_{n+1}(x,k)-w_n(x,k)$$

$$+\left(x-\frac{2}{3}+r_{n+1}(x,k)\right)\mathrm{P}\left(N_{n+1+x}=n+1\right)$$

$$-\left(x-\frac{2}{3}+r_n(x,k)\right)\mathrm{P}\left(N_{n+x}=n\right),$$

Full text

Keywords:: Robust estimate Poisson distribution Quantile

1 Introduction and position of the problem

The Poisson distribution is commonly used for modeling count data. Let $\mathbf{N}_{\lambda}=(N_{1,\lambda},\dots,N_{n,\lambda})$ be a sample of $n\geq 1$ independent and identically distributed random variables distributed as $N_{\lambda}$ a Poisson distribution with parameter $\lambda>0$. Different strategies exist to make the maximum likelihood estimator of $\lambda$ more robust to outliers. For example, specific M-estimators (such as the modified Tukey’s type estimate) for $\lambda$ have been investigated deeply by [5]. The authors also investigate weighted likelihood type estimators and trimmed-mean estimators. In the present paper, we focus our attention on the simplest robust alternative to the maximum likelihood estimator, which is the sample median of $\mathbf{N}_{\lambda}$ and actually on a theoretical problem induced by the use of such an estimator.

To introduce our contribution, let us consider the following standard notation. For a random variable $Y$, we denote by $F_{Y}(\cdot)$ its cumulative distribution function (cdf), by $F_{Y}^{-1}(p)$ its quantile of order $p\in(0,1)$ and by $\mathrm{Me}_{Y}=F_{Y}^{-1}(1/2)$ its theoretical median. Based on a sample $\mathbf{Y}=(Y_{1},\dots,Y_{n})$ of $n$ identically distributed random variables we denote by $\widehat{F}(\cdot;\mathbf{Y})$ the empirical cdf, by $\widehat{F}^{-1}(p;\mathbf{Y})$ the sample quantile of order $p$ given by $\widehat{F}^{-1}(p;\mathbf{Y})=\inf\{x\in\mathbb{R}:p\leq\widehat{F}(x;\mathbf{Y})\}.$ The sample median is simply denoted by $\widehat{\mathrm{Me}}(\mathbf{Y})=\widehat{F}^{-1}(1/2;\mathbf{Y})$. Finally, the density of $Y$ when it exists, is denoted by $f_{Y}$.

Due to the discrete nature of the Poisson distribution, the limiting distribution of $\widehat{\mathrm{Me}}(\mathbf{N}_{\lambda})$ does not follow from standard theory, see e.g. [9] or [7], since it is required that the model possesses a positive density at the true median. To circumvent this problem, one classical strategy introduced by [8], applied to count data by [6] and to the estimation of the intensity of a homogeneous spatial point process by [4], consists in artificially imposing smoothness in the problem through jittering: i.e. we add to each count variable $N_{i,\lambda}$ a random variable $U_{i}\sim\mathcal{U}((0,1))$. Let $\mathbf{Z}_{\lambda}=(Z_{1,\lambda},\dots,Z_{n,\lambda})$, where $Z_{i,\lambda}=N_{i,\lambda}+U_{i}$ for $i=1,\dots,n$, be the sample of independent random variables distributed as $Z_{\lambda}=N_{\lambda}+U$ where $U\sim\mathcal{U}((0,1))$ is independent of $N_{\lambda}$. It can be shown that $Z_{\lambda}$ admits a density almost everywhere, which is given by

[TABLE]

Standard asymptotic theory (e.g. [7]) is now valid: as $n\to\infty$

[TABLE]

in distribution, where $\sigma^{-2}_{\lambda}=4f_{Z_{\lambda}}(\mathrm{Me}_{Z_{\lambda}})^{2}$.

Equation (2) is the source of motivation for the present paper since we are clearly invited to understand how far $\mathrm{Me}_{Z_{\lambda}}$ is from $\lambda$. The study of the median for Poisson and Gamma distributions has a long story, see [3] and the references therein. We can even go back to an old and outstanding formula by Ramanujan, see [3, Equation (3)]. Among several results, [3] proves a conjecture proposed by [2] which is that for every $\lambda>0$

[TABLE]

It is worth mentioning that these bounds are optimal, in the sense that there exists at least one value of $\lambda$ for which the lower-bound or upper-bound is reached. [1] complete this work and prove that asymptotically as $\lambda\to~{}\infty$

[TABLE]

Going back to $\mathrm{Me}_{Z_{\lambda}}$, using these results one easily deduces that

[TABLE]

Such a result is definitely pessimistic since the contribution of this paper is to show that we have the surprising and unexpected following result: $\mathrm{Me}_{Z_{\lambda}}$ is actually very close to $\lambda+1/3$. Even more, our main result implies that, by denoting $\delta_{\lambda}=\mathrm{Me}_{Z_{\lambda}}-\lambda-1/3$

[TABLE]

The latter results suggests us to propose $\hat{\lambda}^{\mathrm{J}}=\mathrm{Me}(\mathbf{Z}_{\lambda})-1/3$ as a new estimator for $\lambda$.

The rest of the paper is organized as follows. Section 2 presents our main result and provides a sketch of the proof while Section 3 illustrates this result. We investigate statistical properties of $\hat{\lambda}^{\mathrm{J}}$ and compare its performances with the maximum likelihood estimator and the Tukey’s modified estimator proposed by [5]. Finally, we show that $\hat{\lambda}^{\mathrm{J}}$ does not suffer from computational problems and can still be used with very large datasets. The proof of our main result relies upon simple technical lemmas which are postponed to Appendix.

2 Main result

We consider the notation introduced in the previous section. Let us first mention that the cumulative distribution function $F_{Z_{\lambda}}$ is given for any $t\geq 0$ by

[TABLE]

whereby it can be checked that $Z_{\lambda}$ indeed admits a density almost everywhere and that this density is given by (1). Our main result is based on the empirical finding depicted in Figure 1. Figure 1 (a) illustrates that for any $\lambda>0$, $\mathrm{Me}_{N_{\lambda}}-\lambda\in[-\log(2),1/3]$, while Figure 1 (b)-(c) illustrate that indeed $\mathrm{Me}_{Z_{\lambda}}\approx\lambda+1/3$. Note that to evaluate $\mathrm{Me}_{Z_{\lambda}}$ we use root-finding algorithm for the function $|F_{Z_{\lambda}}(\cdot)-1/2|$. We now present our main result.

Theorem 2.1**.**

Let $Z_{\lambda}=N_{\lambda}+U$ where $N_{\lambda}$ and $U$ are two independent random variables respectively distributed as a Poisson distribution with parameter $\lambda>0$ and a uniform distribution on $(0,1)$. Then, as $\lambda\to\infty$, the median, $\mathrm{Me}_{Z_{\lambda}}$ of $Z_{\lambda}$ satisfies

[TABLE]

where $\mathcal{H}:[0,1]\to\mathbb{R}$ is the continous function given by

[TABLE]

Equation (3) is easily deduced since (see also Figure 1(d)) we can check that

[TABLE]

Theorem 1 ensues from the following proposition for which we provide a sketch of the proof.

Proposition 2.1**.**

Let $n\geq 1$ and $x\in[0,1)$. For any $\varepsilon>0$, there exists $n_{0}\in\mathbb{N}$ such that for all $n\geq n_{0}$

[TABLE]

Proof.

Let $k\in\mathbb{R}$, and $(w_{n}(x,k))_{n\geq 1}$ be the sequence given by $w_{n}(x,k)=\mathrm{P}(Z_{n+x}\leq n+x+1/3+k/(n+x))$ and $V_{n}=Z_{n+x}-(n+x+1/3+k/(n+x))$. First,

[TABLE]

where $\stackrel{{\scriptstyle d}}{{=}}$ stands for equality in distribution, $N_{1},\ldots,N_{n},N_{x}$ and $U$ are independent random variables, $N_{i}\sim P(1)$ for $i=1,n$, $N_{x}\sim P(x)$ and $U\sim\mathcal{U}([0,1])$. Central limit theorem and Slutsky’s lemma show that $V_{n}/\sqrt{n}\to N(0,1)$ in distribution as $n\to\infty$, whereby we deduce that for any $x\in[0,1)$ and $k\in\mathbb{R}$, $w_{n}(x,k)\to 1/2$ as $n\to\infty$.

Now, assume $x$ and $k$ are such that for $n$ sufficiently large $(w_{n}(x,k))_{n\geq 1}$ is increasing or decreasing then obviously

[TABLE]

or

[TABLE]

So, the rest of the proof simply consists in proving that the sequence $(w_{n}(x,k))_{n\geq 1}$ is monotonic for $n$ sufficiently large. We start by noting that the discontinuity at $x=2/3$ of the function $\mathcal{H}$ comes from the definition of $w_{n}(x,k)$. Lemma A.1 shows in particular that for $n$ sufficiently large

• •

if $x\in(0,2/3)$ or $x=2/3$ and $k<0$,

[TABLE]

• •

if $x\in(2/3,1)$ or $x=2/3$ and $k\geq 0$

[TABLE]

Then, we define $\Delta_{n}(x,k)=(n+1)!/g_{n}(n+1+x)\{w_{n+1}(x,k)-w_{n}(x,k)\}$ where for any $n\geq 1$ and $u\in\mathbb{R}$, $g_{n}(u)=\mathrm{e}^{-u}u^{n}$. Lemma A.3, which is based on simple but lengthy Taylor expansions, shows that for any $x\in[0,1)$

[TABLE]

So, if we set $k=\mathcal{H}(x)+\varepsilon$, $(w_{n}(x,k))_{n\geq 1}$ is a decreasing sequence for $n$ sufficiently large which, from (9) leads to the upper-bound of (7). In the same way, if we set $k=\mathcal{H}(x)-\varepsilon$, $(w_{n}(x,k))_{n\geq 1}$ is an increasing sequence for $n$ sufficiently large which, from (8) leads to the lower-bound of (7). ∎

Our main result has a simple statistical application.We suggest to estimate $\lambda$ by $\hat{\lambda}^{\mathrm{J}}=\widehat{\mathrm{Me}}(\mathbf{Z}_{\lambda})-1/3$: $\hat{\lambda}^{\mathrm{J}}$ is almost an unbiased estimator of $\lambda$, and we can use the approximation

[TABLE]

where $1/\sigma_{\lambda}=2\mathrm{P}(N_{\lambda}=\lfloor\mathrm{Me}_{Z_{\lambda}}\rfloor)$. Note that $1/\sigma_{\lambda}$ can simply be estimated by $1/\hat{\sigma}_{\lambda}=2\mathrm{P}(N_{\hat{\lambda}^{\mathrm{J}}}=\lfloor\hat{\lambda}^{\mathrm{J}}+1/3\rfloor)$. When $\lambda$ is large, we can even use Stirling’s formula to approximate $\sigma_{\lambda}\sim\sqrt{\pi\lambda/2}$ which is then simply estimated by $\sqrt{\pi\hat{\lambda}^{\mathrm{J}}/2}$. Therefore, $\sqrt{\pi/2}$ represents the ratio of asymptotic standard deviations of $\hat{\lambda}^{\mathrm{J}}$ (when $\lambda$ is large) and the maximum likelihood estimator. We can wonder where this $\sqrt{\pi/2}$ comes from: actually this ratio is also the ratio of standard deviations of the sample median to the sample mean when we consider a sample of i.i.d. Gaussian random variables with mean 0 and variance 1.

We end this section by stressing on the simplicity of the estimator $\hat{\lambda}^{\mathrm{J}}$. We do not resist to provide the R instruction to evaluate it based on a sample stored in a vector y:

median(y+runif(length(y)))-1/3

3 Numerical results

3.1 Performances of $\hat{\lambda}^{\mathrm{J}}$ without outliers

For 100 values of $\lambda$ between 1 and 10, we generate 10000 replications of samples of Poisson distribution of size $n$ with parameter $\lambda$. We consider $n=50$ and $n=200$. For each sample $\mathbf{N}_{\lambda}$, we evaluate the maximum likelihood estimate that we denote in the sequel by $\hat{\lambda}^{\mathrm{MLE}}$, $\hat{\lambda}^{\mathrm{J}}$ and $\widehat{\mathrm{Me}}(\mathbf{N}_{\lambda})$. Figure 2 reports empirical biases in terms of $\lambda$. As expected, $\hat{\lambda}^{\mathrm{MLE}}$ and $\hat{\lambda}^{\mathrm{J}}$ are almost unbiased while $\widehat{\mathrm{Me}}(\mathbf{N}_{\lambda})$ has some bias which doesn’t disappear with large $\lambda$ or large $n$. Figure 3 shows $\mathrm{RMSE}(\hat{\lambda}^{\mathrm{J}})/\mathrm{RMSE}(\hat{\lambda}^{\mathrm{MLE}})$ which is the ratio of the root mean squared error (RMSE) of $\hat{\lambda}^{\mathrm{J}}$ to the one of $\hat{\lambda}^{\mathrm{MLE}}$. Obviously the MLE outperforms $\hat{\lambda}^{\mathrm{J}}$ and we observe that the ratio of RMSE is close to $\sqrt{\pi/2}$ when $\lambda$ gets large. Finally, to confirm the estimation of the standard error of $\hat{\lambda}^{\mathrm{J}}$ and its asymptotic normality, we investigate the random variable

[TABLE]

for which its distribution should be close to standard normal distribution. For each value of $\lambda$ considered, Figure 4 depicts for $n=50$ and $n=200$ the 100 normal probability plots. Actually, we only represent the fitted linear regression models and the expected theoretical line $y=x$. We conclude that for every value of $\lambda$, $\Delta_{\lambda}$ seems indeed well-approximated by a $N(0,1)$ distribution.

3.2 Simulation study in presence of outliers

In this section, we compare three estimators of $\lambda$: $\hat{\lambda}^{\mathrm{MLE}}$ which serves as a baseline, $\hat{\lambda}^{\mathrm{J}}$ and the Tukey’s modifed estimator proposed by [5]. This estimator, denoted by $\hat{\lambda}^{\mathrm{Tukey}}$ in the sequel, is an $M$-estimator, see e.g. [9], with objective function

[TABLE]

where $a=a(\lambda,k)$ is such that $\mathrm{E}(\psi_{k,a}(Y,\lambda))=0.$ [5] considered several other estimators (other version of $M$-type estimators, weighted likelihood type estimators, etc), made an extensive simulation study and concluded that in many situations ${\hat{\lambda}}^{\mathrm{Tukey}}$ was the best one. To tune the constant $k$ and thus the corrective term $a=a(\lambda,k)$, we follow the suggestion by [5] and set the constant $k=6$. Given a first estimate of $\lambda$, a first corrective term $a$ is found which serves as a first $M$-estimation and so on. The algorithm is stopped when the difference between two successive estimates of $\lambda$ does not exceed $10^{-4}$.

The simulation model we consider is an additive outliers type model where we assume to observe $\widetilde{N}_{i,\lambda}$ given by

[TABLE]

where $P(\varepsilon_{i}=1)=1-\pi$ and where $\pi$ corresponds to the proportion of outliers and $h$ is a constant. For a given $\pi$, we consider the signal-to-noise ratio defined in decibels by $SNR=10\log_{10}(\frac{\lambda}{h\pi(1-\pi)})$ and set $\sqrt{h}$ (as an integer value) such that $\mathrm{SNR}=-10$ (db).

Figure 5 reports empirical biases and RMSE in terms of $\pi\in\{0\%,1\%,5\%,10\%,20\%\}$ for estimators of $\lambda=5$. These Monte-Carlo results are based on 10000 replications from the model (12) with $n=50$ and $n=200$. As expected, the MLE gets quickly biased as soon as $\pi>0$ which makes its RMSE very high. When $\pi$ is not too large, $\hat{\lambda}^{\mathrm{J}}$ seems less biased than $\hat{\lambda}^{\mathrm{Tukey}}$. However, the latter one is more efficient which explains why the RMSE of the $\hat{\lambda}^{\mathrm{Tukey}}$ is smaller than the one of $\hat{\lambda}^{\mathrm{J}}$. It is to be noticed that this difference tends to decrease when $n$ increases. When $\pi=10\%$ or $20\%$, $\hat{\lambda}^{\mathrm{Tukey}}$ is much less biased than $\hat{\lambda}^{\mathrm{J}}$ and still has a smaller variance.

As a conclusion of this simulation study, it turns out that the naive and very simple estimator $\hat{\lambda}^{\mathrm{J}}$ behaves nicely compared to very efficient estimators such as the Tukey’s modified estimator, when the proportion of outliers is low. When this proportion increases, the performances of the median-based estimator $\hat{\lambda}^{\mathrm{J}}$ degrade, whereas $\hat{\lambda}^{\mathrm{Tukey}}$ still remains efficient.

3.3 Computational time

Table 1 reports average computational time of four estimators of $\lambda$ ($\hat{\lambda}^{\mathrm{MLE}}$, $\widehat{\mathrm{Me}}(\mathbf{N}_{\lambda})$, $\hat{\lambda}^{\mathrm{J}}$ and $\hat{\lambda}^{\mathrm{Tukey}}$) in terms of the sample size. The estimates are implemented in R on a laptop using a 2,9 GHz Intel Core i5 process. The MLE is obviously the cheapest one. Evaluating $\hat{\lambda}^{\mathrm{J}}$ is approximately twice more expensive than evaluating $\widehat{\mathrm{Me}}(\mathbf{N}_{\lambda})$. This factor is due to the generation of uniform distributions. However, the computational time is very reasonable even for very large datasets compared to the Tukey’s modified estimate which is unfeasible for $n=10^{9}$ due to memory storage.

Acknowledgements

The authors are sincerely grateful to H. Elsaied and R. Fried for discussions and for sharing the R code implementing the Tukey’s modified $M$-estimator. The research of J.-F. Coeurjollly is supported by the Natural Sciences and Engineering Research Council of Canada.

Appendix A Technical lemmas

Appendix gathers technical lemmas, used in the proof of Proposition 1.

Lemma A.1**.**

*Let $k\in\mathbb{R}$, $x\in[0,1)$ and $r_{n}(x,k)=k/(n+x)$. Then the sequence $w_{n}(x,k)$ reads as follows:

(i) If $x+r_{n}(x,k)\in[-1/3,2/3)$*

[TABLE]

(ii) If $x+r_{n}(x,k)\in[2/3,5/3)$

[TABLE]

Let $k\in\mathbb{R}$, then for $n$ sufficiently large, if $x\in[0,2/3)$ or $x=2/3$ and $k<0$, then $w_{n}(x,k)$ reads as in (13). In the same way, if $x\in(2/3,1)$ or $x=2/3$ and $k\geq 0$, then $w_{n}(x,k)$ reads as in (14). The proof of Lemma A.1 is omitted as it derives easily from (4).

Lemma A.2**.**

Let $k\in\mathbb{R}$, $x\in[0,1)$ and $r_{n}(x,k)=k/(n+x)$. Let $\Delta_{n}(x,k)$ be given by

[TABLE]

*There exists $n_{0}\in\mathbb{N}$ such that for all $n\geq n_{0}$, we have the two following cases:

(i) If $x+r_{n}(x,k)\in[-1/3,2/3)$*

[TABLE]

where $c_{n}(\cdot,x):[0,1]\to\mathbb{R}$ is defined by

[TABLE]

(ii) If $x+r_{n}(x,k)\in[2/3,5/3)$

[TABLE]

Proof.

(i) Using the Poisson-Gamma relation $\mathrm{P}(N_{\lambda}\leq n)=\frac{1}{n!}\int_{\lambda}^{\infty}g_{n}(u)\mathrm{d}u$ with $g_{n}(u)=\mathrm{e}^{-u}u^{n}$ and Lemma A.1 (i), we can rearrange the difference $w_{n+1}(x,k)-w_{n}(x,k)$ as

[TABLE]

which leads to the result after little algebra by noticing that

[TABLE]

(ii) Using the Poisson-Gamma relation and Lemma A.1 (ii), we can rearrange the difference $w_{n+1}(x,k)-w_{n}(x,k)$ as

[TABLE]

which leads to the result after little algebra by noticing that

[TABLE]

and

[TABLE]

∎

Lemma A.3**.**

Let $k\in\mathbb{R}$, $x\in[0,1)$ and $r_{n}(x,k)=k/(n+x)$, for any $k\in\mathbb{R}$ and $x\in[0,1)$, then for $n$ sufficiently large we have

[TABLE]

where $\mathcal{H}$ is the function given by (6).

Proof.

Let $k\in\mathbb{R}$, then for any $x\in[0,1)$, there exists $n_{0}\in\mathbb{N}$ such that for $n\geq n_{0}$, either (i) $x+r_{n}(x,k)\in[-1/3,2/3)$ or (ii) $x+r_{n}(x,k)\in[2/3,5/3)$. In the sequel, we consider both cases and expand the expression of $\Delta_{n}(x,k)$ given by Lemma A.2. The following expansions extensively make use of auxiliary results gathered in Lemma A.4.

(i) Case $x+r_{n}(x,k)\in[-1/3,2/3)$.

[TABLE]

(ii) Case $x+r_{n}(x,k)\in[2/3,5/3)$.

[TABLE]

∎

Lemma A.4**.**

*Let $v,x\in[0,1]$, then we have the following expansions as $n\to\infty$:

(i)*

[TABLE]

(ii)

[TABLE]

(iii)

[TABLE]

(iv)

[TABLE]

(v)

[TABLE]

(vi)

[TABLE]

(vii)

[TABLE]

Proof.

(i) Using Taylor expansions, we have

[TABLE]

Let $\rho_{n}(x)=\frac{n+1}{n+1+x}=1-\frac{x}{n+1+x}$, then

[TABLE]

(ii) It can be easily deduced from Equation (17).

(iii) Starting from Equation (17), since $\int_{0}^{1}(1-v)^{k}dv=\frac{1}{k+1}$,we have

[TABLE]

(iv) Using Taylor expansions, we get

[TABLE]

We expand $c_{n}(0,x)\rho^{\prime}_{n}(x)$ using Equation (19).

[TABLE]

whereby the result is deduced.

(v) Let $x+r_{n}(x,k)\in[-1/3,2/3)$, then

[TABLE]

(vi) We use the Taylor expansion of $\frac{n+1+x}{n+2}$, which is $1+\frac{x-1}{n+1+x}+\frac{(x-1)^{2}}{(n+1+x)^{2}}+o\left(\frac{1}{n^{2}}\right)$ to deduce that

[TABLE]

(vii) Now, let $x+r_{n}(x,k)\in[2/3,5/3)$. We use the Taylor expansion of $\frac{n+1+x}{n+2}$, which is $1+\frac{x-1}{n+1+x}+o\left(\frac{1}{n^{2}}\right)$ to deduce that

[TABLE]

∎