Find estimator for $\lambda$ for $X\sim \operatorname{Poisson}(\lambda)$ using the 2nd method of moment

parameter estimationpoisson distributionstatistical-inferencestatistics

So an estimator for for $\lambda$ using the first method of moment is simply the sample mean:

$$\mu_1 = E(X) = \lambda = \bar{X} = \hat{\mu_2}$$

For the second moment,

$$\operatorname{Var}(X) = E(X^2) – E(X)^2$$
$$\lambda = E(X^2) – \lambda^2$$
$$\mu_2 = E(X^2) = \lambda + \lambda^2$$
$$s^2 = \frac{1}{n}\sum_{i=1}^n (X_i – \bar{X})^2$$
$$s^2 = \frac{1}{n}\sum_{i=1}^n X_i^2 – \bar{X}^2$$
$$\frac{1}{n}\sum_{i=1}^n X_i^2 = s^2 + \bar{X}^2$$
$$\hat{\mu_2} = s^2 + \bar{X}^2$$
Making $\mu_2 = \hat{\mu_2}$, we get:

$$\mu_2 = \lambda + \lambda^2 = s^2 + \bar{X}^2 = \hat{\mu_2}$$
$$ \lambda + \lambda^2 = s^2 + \bar{X}^2$$

From here I've no idea how to factorise it so I have $\lambda$ all on one side.

Best Answer

So you have three potential estimators:

  • The first, $\bar X$ is unbiased and has the smallest variance among unbiased estimators.
  • The second, $S^2$ is also unbiased, but has much larger variance than $\bar X.$
  • The third is given in @JimB's comment. It is unbiased and its variance is not much larger than the variance of $\bar X.$

in the R simulation below, SDs of estimators are given instead of their variances. With $m = 100\,000$ iterations we can expect about two significant digits of accuracy.

set.seed(2020)
m = 10^5;  n = 10;  lam = 3
x = rpois(m*n, lam)
MAT = matrix(x, nrow=m)
a = rowMeans(MAT)
v = apply(MAT, 1, var)
est.1 = a
mean(est.1);  sd(est.1)
[1] 3.0008      # aprx 3 (unbiased)
[1] 0.5484283   # smallest variability
est.2 = v
mean(est.2);  sd(est.2)
[1] 2.99409     # aprx 3 (unbiased)
[1] 1.518167    # largest variability of the 3
est.3 = .5*(sqrt(4*v+4*a^2+1)-1)
mean(est.3);  sd(est.3)
[1] 2.994443    # aprx 3 (unbiased)
[1] 0.581836    # slightly more variable than mean

enter image description here

par(mfrow=c(3,1))
mn = min(est.1, est.2, est.3)
mx = max(est.1, est.2, est.3)
 hist(est.1, prob=T, xlim=c(mn,mx), col="skyblue2")
 hist(est.2, prob=T, br=60, xlim=c(mn,mx), col="skyblue2")
 hist(est.3, prob=T, xlim=c(mn,mx), col="skyblue2")
par(mfrow=c(1,1))
Related Question