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|\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|\lambda)$.
This is quite nice, however it requires that I have samples directly from $\text{Poisson}(k|\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|\lambda)$ in place of $\lambda$, because that will not produce an unbiased estimate of $\text{Poisson}(k|\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|\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!
Best Answer
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|\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|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|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).