Looking for unbiased estimators for Poisson probabilities

Question

I am looking for unbiased estimators for Poisson probabilities. That is, some estimator $\hat{g}(k)$ such that

$$ E( \hat{g}(k) ) = \text{Poisson}(k\mid\lambda) $$

I discovered one in this old paper:

https://www.jstor.org/stable/1266576 ("Minimum Variance Unbiased Estimators for Poisson Probabilities"), which is

$$ \hat{g}(k) = \binom{x}{k} \left(\frac{1}{N} \right)^k \left(1 - \frac{1}{N}\right)^{x-k}$$

where $x = \sum_{i=1}^N k_i $, where the $k_i$ are samples from $\text{Poisson}(k\mid\lambda)$.

This is quite nice, however it requires that I have samples directly from $\text{Poisson}(k\mid\lambda)$, which I don't have, or at least I don't think I can get these. Instead, what I have are ways of estimating $\lambda$, in particular I can get this by Monte Carlo integration.

Actually, to be more specific, I can get $p$ by Monte Carlo integration, where then $\lambda = p M$ for some known factor $M$. So this estimate of $\lambda$ is unbiased, however I can't just plug it into $\text{Poisson}(k\mid\lambda)$ in place of $\lambda$, because that will not produce an unbiased estimate of $\text{Poisson}(k\mid\lambda)$, and I really need an unbiased estimate for my application. I'm not quite sure if the Monte Carlo estimate of $\lambda$ can be modified to instead sample from $\text{Poisson}(k\mid\lambda)$, but at least it doesn't look to me like I can do that. But perhaps something like that is a possibility?

I am not very familiar with the theory on constructing unbiased estimators of things, so I am not really sure how to proceed here. Any pointers would be much appreciated!

Answer 1

I'm not sure I understand your question, but I will guess that you mean you want an unbiased estimator of $\lambda^k e^{-\lambda}/k!$ that is a function of an unbiased estimator of $\lambda.$ If you have i.i.d. observations $X_1,\ldots,X_n$ from the Poisson distribution with expected value $\lambda,$ then $\overline X = (X_1+\cdots+X_n)/n$ is the best unbiased estimator of $\lambda.$

And $\overline X$ is a sufficient statistic for this family of distributions, i.e. the conditional distribution of $X_1,\ldots,X_n$ given $\overline X$ does not depend on $\lambda.$

The indicator random variable $I_k= \left. \begin{cases} 1 & \text{if } X_1=k \\ 0 & \text{if } X_1\ne k \end{cases} \right\} $ is an unbiased estimator of $\lambda^k e^{-\lambda}/k!.$

So we seek $\operatorname E\left( I_k \mid \overline X\right).$ \begin{align} & \operatorname E\left( I_k\mid \overline X = x/n \right) = \operatorname E\left( I_k \mid X_1+\cdots + X_n = x \right) \\[8pt] = {} & \Pr( X_1=k \mid X_1+\cdots + X_n = x) \\[12pt] = {} & \frac{\Pr(X_1=k\ \&\ X_1+\cdots + X_n = x)}{\Pr(X_1+\cdots+X_n=x)} \\[12pt] = {} & \frac{\Pr( X_1=k\ \&\ X_2+\cdots + X_n = x-k)}{\Pr(X_1+\cdots+ X_n=x)} \\[12pt] = {} & \frac{\Pr( X_1=k)\Pr(X_2+\cdots + X_n = x-k)}{\Pr(X_1+\cdots+ X_n=x)} \\[12pt] = {} & \frac{\Big(\lambda^k e^{-\lambda}/k!\Big) \Big( \big((n-1)\lambda\big)^{x-k} e^{-(n-1)\lambda}/(x-k)! \Big)}{(n\lambda)^x e^{-n\lambda}/x!} \\[12pt] = {} & \binom x k \left. \left( \frac{n-1}n \right)^x \right/ (n-1)^k. \end{align}

Letting $\widehat\lambda = \overline X,$ we then have $$ \operatorname E\left( I_k \mid \widehat\lambda \right) = \left. \binom{n\widehat\lambda} k \left( \frac{n-1}n \right)^{n\widehat\lambda} \right/ (n-1)^k. $$

That should be a fairly good unbiased estimator of $\Pr(X_1=k) $ $= \lambda^k e^{-\lambda}/k!.$ But I haven't looked at what its variance is or whether it is better of worse, in the sense of mean squared error, than the maximum-likelihood esimator, which is just $\widehat\lambda^{\,k} e^{-\widehat\lambda}/k!.$

Postscript: I just ran the following commands in R: \begin{align} & \texttt{x <- rpois(1000,6.3)} \\ & \texttt{6.3^5 * exp(-6.3)/120} \\[5pt] & \text{The value returned was }0.151868. \\[5pt] & \texttt{choose(6337,5) * (999/1000)^6337 / 999^5} \\[5pt] & \text{This is the value of the unbiased estimator.} \\ & \text{It was } 0.1507382. \\[5pt] & \texttt{mean(x)^5 * exp(-mean(x))/120} \\[5pt] & \text{This gives the value of the maximum-likelihood} \\ & \text{estimator. It was } 0.1507, \text{ with rounding error less} \\ & \text{than } 4\times10^{-8}. \text{ (somewhat surprising)} \end{align} In this case there is not much difference between the unbiased estimator and the MLE.

PS: The argument above assumes $n>1.$ When $n=1,$ the bottom line you see above is $0/0.$ But in that case $\operatorname E\left( I_n\mid \overline X\right)$ is still defined and is equal in every instance to $1$ or $0.$ But the MLE for $\lambda$ in that case is just $X_1,$ and that for $\Pr(X_1=k)$ is just $X_1^k e^{-X_1}/k!.$ Clearly the MLE performs better than the unbiased estimator in that case.

Answer 2

Ok well I made some progress on this, and I guess it technically answers my question, although I have some dissatisfaction with it, as I will explain.

The estimator I gave in the OP can indeed be improved to the cases I am interested in, as follows.

Suppose that the Poisson distribution we are interested in can be written as

$$\Pr(k\mid\lambda=Np) = \frac{e^{-Np}(Np)^k}{k!}$$

as in, this is the Poisson limit of N Bernoulli draws with success probability $p$. Now, $p$ is unknown, but we can estimate it via Monte Carlo methods. Suppose these can be describe as draws from another Poisson distribution

$$\Pr(x\mid np) = \frac{e^{-np}(np)^x}{x!}$$

where $x$ is the number of 'successes' in the Monte Carlo simulation. We can then use the MC result to obtain an unbiased estimator for $\Pr(k\mid Np)$ by following closely the proof in the paper I linked in the question. Essentially, we can expand $e^{-Np}$ in a power series, and then use the MC samples to estimate all powers of $p$ unbiasedly as $p^i \sim x!/(n^i (x-i)!)$ , which then gives the estimator we want. Here is the proof:

$\begin{aligned}\frac{e^{-Np}(Np)^k}{k!} &= \frac{(Np)^k}{k!}\sum_{\alpha=0}^\infty \frac{(-Np)^\alpha}{\alpha!} \\ &= \frac{(N)^k}{k!}\sum_{\alpha=0}^\infty \frac{(-N)^\alpha}{\alpha!} p^{\alpha + k} \\ &\sim \frac{(N)^k}{k!}\sum_{\alpha=0}^\infty \frac{(-N)^\alpha}{\alpha!} \left(\frac{1}{n}\right)^{\alpha+k}\frac{x!}{(x-k-\alpha)!} \\ &= \binom{x}{k} \left(\frac{N}{n}\right)^k \sum_{\alpha=0}^\infty \binom{x-k}{\alpha}\left(-\frac{N}{n}\right)^\alpha \\ &= \binom{x}{k} \left(\frac{N}{n}\right)^k \left(1 - \frac{N}{n}\right)^{x-k} \end{aligned}$

Which indeed seems to work out numerically when I test it. So technically that's my question answered. However, there are a couple of unsatisfying properties here which I wonder if they can be avoided somehow, or if not, I'd like to know why.

First, the estimator is zero or maybe undefined for $x<k$. This means that you need to simulate more successful events than are actually observed, which, if $p$ is small, can require quite a lot of Monte Carlo work. Secondly, it requires that $n>N$ in order for it to be non-negative (and maybe to make sense at all), so again, you need more trials in the simulation than the real experiment, which is a problem if $N$ is big. Neither of these requirements exist in order to use the maximum likelihood estimator from the simulation, which is just $\hat{p}=x/n$, however that estimator is biased.

Is there no way to have my cake and eat it too here? Is the cost of having an unbiased estimate really what this estimator requires, or is there some other estimator that has less strict requirements on the MC distribution? I think the estimator derived here is the minimum variance unbiased estimator given our MC samples, but actually I would be happier with an estimator with larger variance but which worked for $x<k$ and $n<N$ (while still being unbiased).