13
$\begingroup$

So, I have a random process generating log-normally distributed random variables $X$. Here is the corresponding probability density function:

Figure representing a lognormal probability density function

I wanted to estimate the distribution of a few moments of that original distribution, let's say the 1st moment: the arithmetic mean. To do so, I drew 100 random variables 10000 times so that I could calculate 10000 estimate of the arithmetic mean.

There are two different ways to estimate that mean (at least, that's what I understood: I could be wrong):

  1. by plainly calculating the arithmetic mean the usual way: $$\bar{X} = \sum_{i=1}^N \frac{X_i}{N}.$$
  2. or by first estimating $\sigma$ and $\mu$ from the underlying normal distribution: $$\mu = \sum_{i=1}^N \frac{\log (X_i)}{N} \quad \sigma^2 = \sum_{i=1}^N \frac{\left(\log (X_i) - \mu\right)^2}{N}$$ and then the mean as $$\bar{X} = \exp(\mu + \frac{1}{2}\sigma^2).$$

The problem is that the distributions corresponding to each these estimates are systematically different:

The two estimators give different distributions as shown on the picture.

The "plain" mean (represented as the red dashed line) provides generally lower values than the one derived from the exponential form (green plain line). Though both means are calculated on the exact same dataset. Please note that this difference is systematic.

Why aren't these distributions equal?

$\endgroup$
  • $\begingroup$ what are your true parameters for $\mu$ and $\sigma$? $\endgroup$ – Christoph Hanck Jun 30 '16 at 13:58
  • $\begingroup$ $\mu = 3$ and $\sigma = 1.5$, but please note that I am interested in estimating these parameters, hence the Monte-Carlo approach instead of calculating the thing from these raw numbers. $\endgroup$ – JohnW Jun 30 '16 at 14:30
  • $\begingroup$ sure, this is for replication your results. $\endgroup$ – Christoph Hanck Jun 30 '16 at 14:33
  • 4
    $\begingroup$ Interestingly, this phenomenon has nothing to do with lognormality. Given positive numbers $x_i$ with logarithms $y_i$, it is well known their arithmetic mean (AM) $\sum x_i/n$ is never less than their geometric mean (GM) $\exp(\sum y_i/n)$. In the other direction, the AM is never greater than the GM multiplied by $\exp(s_y^2/2)$ where $s_y^2$ is the variance of the $y_i$. Thus, the dotted red curve must lie to the left of the solid green curve for any parent distribution (describing positive random numbers). $\endgroup$ – whuber Jun 30 '16 at 15:22
  • $\begingroup$ If much of the mean comes from a tiny probability of huge numbers, a finite sample arithmetic mean may underestimate the population mean with high probability. (In expectation it's unbiased, but there's a large probability of a small underestimate and a small probability of a large over estimate.) This question may also relate to this one: stats.stackexchange.com/questions/214733/… $\endgroup$ – Matthew Gunn Jun 30 '16 at 16:16
12
$\begingroup$

The two estimators you are comparing are the method of moments estimator (1.) and the MLE (2.), see here. Both are consistent (so for large $N$, they are in a certain sense likely to be close to the true value $\exp[\mu+1/2\sigma^2]$).

For the MM estimator, this is a direct consequence of the Law of large numbers, which says that $\bar X\to_pE(X_i)$. For the MLE, the continuous mapping theorem implies that $$ \exp[\hat\mu+1/2\hat\sigma^2]\to_p\exp[\mu+1/2\sigma^2],$$ as $\hat\mu\to_p\mu$ and $\hat\sigma^2\to_p\sigma^2$.

The MLE is, however, not unbiased.

In fact, Jensen's inequality tells us that, for $N$ small, the MLE is to be expected to be biased upwards (see also the simulation below): $\hat\mu$ and $\hat\sigma^2$ are (in the latter case, almost, but with a negligible bias for $N=100$, as the unbiased estimator divides by $N-1$) well known to be unbiased estimators of the parameters of a normal distribution $\mu$ and $\sigma^2$ (I use hats to indicate estimators).

Hence, $E(\hat\mu+1/2\hat\sigma^2)\approx\mu+1/2\sigma^2$. Since the exponential is a convex function, this implies that $$E[\exp(\hat\mu+1/2\hat\sigma^2)]>\exp[E(\hat\mu+1/2\hat\sigma^2)]\approx \exp[\mu+1/2\sigma^2]$$

Try increasing $N=100$ to a larger number, which should center both distributions around the true value.

See this Monte Carlo illustration for $N=1000$ in R:

enter image description here

Created with:

N <- 1000
reps <- 10000

mu <- 3
sigma <- 1.5
mm <- mle <- rep(NA,reps)

for (i in 1:reps){
  X <- rlnorm(N, meanlog = mu, sdlog = sigma)
  mm[i] <- mean(X)

  normmean <- mean(log(X))
  normvar <- (N-1)/N*var(log(X))
  mle[i] <- exp(normmean+normvar/2)
}
plot(density(mm),col="green",lwd=2)
truemean <- exp(mu+1/2*sigma^2)
abline(v=truemean,lty=2)
lines(density(mle),col="red",lwd=2,lty=2)

> truemean
[1] 61.86781

> mean(mm)
[1] 61.97504

> mean(mle)
[1] 61.98256

We note that while both distributions are now (more or less) centered around the true value $\exp(\mu+\sigma^2/2)$, the MLE, as is often the case, is more efficient.

One can indeed show explicitly that this must be so by comparing the asymptotic variances. This very nice CV answer tells us that the asymptotic variance of the MLE is $$V_t = (\sigma^2 + \sigma^4/2)\cdot \exp\left\{2(\mu + \frac 12\sigma^2)\right\},$$ while that of the MM estimator, by a direct application of the CLT applied to samples averages is that of the variance of the log-normal distribution, $$ \exp\left\{2(\mu + \frac 12\sigma^2)\right\}(\exp\{\sigma^2\}-1) $$ The second is larger than the first because $$ \exp\{\sigma^2\}>1+\sigma^2 + \sigma^4/2, $$ as $\exp(x)=\sum_{i=0}^\infty x^i/i!$ and $\sigma^2>0$.

To see that the MLE is indeed biased for small $N$, I repeat the simulation for N <- c(50,100,200,500,1000,2000,3000,5000) and 50,000 replications and obtain a simulated bias as follows:

enter image description here

We see that the MLE is indeed seriously biased for small $N$. I am a little surprised about the somewhat erratic behavior of the bias of the MM estimator as a function of $N$. The simulated bias for small $N=50$ for MM is likely caused by outliers that affect the non-logged MM estimator more heavily than the MLE. In one simulation run, the largest estimates turned out to be

> tail(sort(mm))
[1] 336.7619 356.6176 369.3869 385.8879 413.1249 784.6867
> tail(sort(mle))
[1] 187.7215 205.1379 216.0167 222.8078 229.6142 259.8727 
$\endgroup$
  • $\begingroup$ Ah okay. It really did not occur to me that one method could be more efficient than the other given the same data. So I could say that the MLE solution converges faster with respect to $N$ than the other method if I understood correctly. Thanks! $\endgroup$ – JohnW Jun 30 '16 at 14:36
  • 1
    $\begingroup$ I made a little edit about the bias. For $N=100$ the bias is indeed negative for the MM estimator, but that does not seem like a general result, see the plot for the bias as a function of $N$. $\endgroup$ – Christoph Hanck Jun 30 '16 at 14:44
  • 2
    $\begingroup$ Well, I am surprised too that there is such a large difference between the two methods, however this example is absolutely perfect to demonstrate why "just averaging stuff" can be awful! $\endgroup$ – JohnW Jun 30 '16 at 14:44
  • 1
    $\begingroup$ @JohnW, I added a little analytical explanation of why the MLE has smaller variance. $\endgroup$ – Christoph Hanck Jul 1 '16 at 11:33
  • 1
    $\begingroup$ The discrepancy stems from the fact that the bias is a finite sample problem, i.e., it vanishes as $N$ goes off to infinity. The asymptotic variance (as the name says) comparison only shows what happens in the limit, as $N\to\infty$. $\endgroup$ – Christoph Hanck Jul 1 '16 at 11:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.