# What is the distribution of the sample variance of the Skellam distribution?

I want to estimate the parameter $\mu$ using the difference between two Poisson distributions with the same parameter, i.e. a Skellam distribution with $\mu_1=\mu_2 = \mu$.

I can calculate the variance of the Skellam distribution as the average variance over multiple samples, however, depending on sample size and some good/bad luck, I have quite a bit of error in my estimate (red crosses in the figure below); the variance should be $2\mu$. If I knew what distribution I was looking for, I could fit that instead to hopefully get a more robust estimate.

Why am I doing this if I appear to know $\mu$? In my application, I only know $\hat{\mu}=\beta_1\mu+\beta_0$, so I need a second estimate of $\mu$ to find out the betas.

-
(1) From what data have you estimated $\beta_1\mu+\beta_0$? (2) Isn't this question really asking how to estimate $\mu$ from iid observations from a Skellam distribution? (Implicitly, it assumes that a method of moments estimator $\hat{\mu}=s^2/2$ ($s^2$ is the sample variance) would be good; but perhaps there are other approaches.) –  whuber Dec 18 '11 at 16:21
Do you need to deal with the possibility of small $\mu$ (less than $2$ or so)? If not, then estimates developed from a Normal approximation might work well (skewness is always zero and excess kurtosis is only $1/(2\mu)$). For largish $\mu$ (around $10$ or larger) it becomes worthwhile considering using order statistics, such as the IQR, as the basis for estimates of $\mu$. The relatively heavy tails of the Skellam contraindicate estimates based on higher moments; one should be inclined to look towards order statistics or M-estimators. –  whuber Dec 18 '11 at 16:49
@whuber: (1) I would like to estimate $\mu$. I can get an estimate for $\hat{\mu}$ from my data, but in order to find the $\beta$s, I need a second estimate of $f(\mu)$ that depends on the $\beta$s in a different way. –  Jonas Dec 19 '11 at 19:28
@whuber: The data comes actually from images, where every every pixel $pix(x,y)~f(poisson(\mu(x,y)) + some signal)$, where I assume $f()$ to be linear. The first estimate of $\hat{\mu}(x,y)$ is obtained by filtering the image, The Skellam distribution is obtained from taking the difference between two images, which eliminates the signal; each sample of the distribution is obtained by calculating the local variance (eg. in a 5x5 window) of the difference image. –  Jonas Dec 19 '11 at 19:36
@whuber: Since I know that both f() and the noise generation process are the same for every pixel, I can lump together the local variances from different pixels with the same $\hat{\mu}$ to get an estimate for the variance (which should be 2$\mu$). However, when I quickly tested this (see figure), I noticed that taking the average of the distribution does not seem to be very robust. In general, I like using order statistics, though the wiki entry for the Skellan distribution wasn't helpful (and I don't know how to calculate the IQR myself). –  Jonas Dec 19 '11 at 19:41
show 1 more comment