A mostly traditional question on X validated about the “best” [minimum variance] unbiased estimator of θ² from a Poisson P(θ) sample leads to the Rao-Blackwell solution
![\mathbb{E}[X_1X_2|\underbrace{\sum_{i=1}^n X_i}_S=s] = -\frac{s}{n^2}+\frac{s^2}{n^2}=\frac{s(s-1)}{n^2}](https://s0.wp.com/latex.php?latex=%5Cmathbb%7BE%7D%5BX_1X_2%7C%5Cunderbrace%7B%5Csum_%7Bi%3D1%7D%5En+X_i%7D_S%3Ds%5D+%3D+-%5Cfrac%7Bs%7D%7Bn%5E2%7D%2B%5Cfrac%7Bs%5E2%7D%7Bn%5E2%7D%3D%5Cfrac%7Bs%28s-1%29%7D%7Bn%5E2%7D&bg=000000&fg=B0B0B0&s=0&c=20201002)
and a similar estimator could be constructed for θ³, θ⁴, … With the interesting limitation that this procedure stops at the power equal to the number of observations (minus one?). But, since the expectation of a power of the sufficient statistics S [with distribution P(nθ)] is a polynomial in θ, there is de facto no limitation. More interestingly, there is no unbiased estimator of negative powers of θ in this context, while this neat comparison on Wikipedia (borrowed from the great book of counter-examples by Romano and Siegel, 1986, selling for a mere $180 on amazon!) shows why looking for an unbiased estimator of exp(-2θ) is particularly foolish: the only solution is (-1) to the power S [for a single observation]. (There is however a first way to circumvent the difficulty if having access to an arbitrary number of generations from the Poisson, since the Forsythe – von Neuman algorithm allows for an unbiased estimation of exp(-F(x)). And, as a second way, as remarked by Juho Kokkala below, a sample of at least two Poisson observations leads to a more coherent best unbiased estimator.)
