6

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.

enter image description here

Stephan Kolassa
  • 95,027
  • 13
  • 197
  • 357
Jonas
  • 1,578
  • 1
  • 13
  • 16
  • 1
    (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
  • 1
    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
  • @whuber: Unfortunately, I cannot exclude small values of $\mu$ yet before I have been able to estimate it on a few thousand independent cases to get a handle on how the instrument works. – Jonas Dec 19 '11 at 19:42

2 Answers2

3

Some of the comments have pointed out that there may be better ways to solve your problem that do not involve finding the distribution of the sample variance of the Skellam distribution. Here I will answer the title question irrespective of those other issues.

The exact distribution of the sample variance from a Skellam distribution is complicated, since it is a quadratic function of Poisson random variables. Using the large-sample approximation with the chi-squared distribution (see e.g., O'Neill 2004) you can get the approximate distribution:

$$S_n^2 \sim 2 \mu \cdot \frac{\text{Chi-Sq}(DF_n)}{DF_n} \quad \quad \quad DF_n = \frac{2n \mu^2}{3 \mu^2 + (n-3) / (n-1)}.$$

This should get you a reasonable approximation to the distribution of the sample variance, so long as $n$ is not too small.

Ben
  • 91,027
  • 3
  • 150
  • 376
  • 1
    That's what I love about this community - even after several years, it is possible to get a helpful answer. Thanks a lot! – Jonas May 07 '18 at 08:44
  • I'm only six or so years late! (Hope it is still relevant.) – Ben May 07 '18 at 09:48
-2

I'm not sure what the problem is: sample variance is a random variable. I would guess it is distributed as a (re-scaled) Chi-square for this distribution. For a fixed number of samples, the standard error of the variance statistic will be proportional to the population variance you are trying to estimate. Thus I would expect these histograms to get 'wider' as the population variance increases.

shabbychef
  • 10,388
  • 7
  • 50
  • 93
  • 3
    Yes, it is a random variable. What is the distribution? I know that the variance of the Normal distribution follows a chi-square distribution (http://en.wikipedia.org/wiki/Variance#Distribution_of_the_sample_variance), but does this apply regardless of the underlying distribution? – Jonas Dec 18 '11 at 13:37
  • 1
    In general, the SE of the sample variance depends on the *fourth* moment of the underlying distribution, not just the second moment. – whuber Dec 18 '11 at 16:45