# Median and mean of the sample mean of i.i.d. log-normal

Let $y:=\frac1n\sum_{i=1}^n x_i$, where $\{x_i\}_{i=1}^n$ is a set of i.i.d. random variables, and every $x_i$ has a lognormal distribution $x_i \sim\text{Lognormal}(\mu,\sigma^2)$. Let $\text{Med}[y]$ be the median of $y$. Is the following inequality true $\forall (n,\mu,\sigma)$? $$\text{Med}[y]<\mathbf E[y]$$

Motivation: I am computing the sample mean of the lognormal random variables via Monte Carlo. The sample mean seems tend to concentrate below the mean for large $\sigma$. I am wondering whether this is true for all cases. It is true for a sigle sample. However, there is no explicit formula for the distribution of the mean of finite number of --- not even two -- samples i.i.d. lognormal variables. I have no idea how to prove it.

• A discrete parallel to consider: if the x's are Bernoulli with probability 40%, then Med[x1] < E[x1], Med[x2] < E[x2], but Med[x1+x2] > E[x1+x2]. – Matt F. Jun 6 '17 at 9:59
• @MattF. is there a general name for such comparison problem between mean and median, I thought there is such a name but cannot remember... – Henry.L Jun 6 '17 at 13:21
• @HenryL, I know of no good name for this. If you say "positive Pearson's second skewness", some people might understand you; that skewness is defined as 3(mean-median) / (standard deviation). But I don't recommend that, since I haven't seen any good citation of a text of Pearson where he proposed this particular measure, and I prefer quartile skewness in any case. mathworld.wolfram.com/PearsonsSkewnessCoefficients.html – Matt F. Jun 9 '17 at 14:45

I would like to provide an asymptotic solution, as the sample size $n\rightarrow\infty$.

Using the Fenton-Wilkson empirical asymptotic result in [1], we know that the sample sum of i.i.d. lognormals $$X_i\sim \text{LN}(\mu,\sigma^2)$$ can be approximated by lognormal with different mean and variance $Y=\sum_{i=1}^nX_i$ is distributed as $$Y\sim \text{LN}(\mu_n,\sigma^2_n)$$ where $\sigma_{n}^{2}=\log\left(\frac{e^{\sigma^{2}}-1}{n}+1\right)$ and $\mu_{n}=\log\left(n\cdot e^{\mu}\right)+\frac{1}{2}\left(\sigma_{n}^{2}-\sigma^{2}\right)$. We know from properties of lognormal distributions that $\text{Median}(Y)=e^{\mu_n}$ and $\text{Mean}(Y)=\exp\big(\mu_n+\frac{\sigma^2_n}{2}\big)$and therefore $\text{Median}(Y)\le \text{Mean}(Y)$ is asymptotically proved since $\sigma^2_n\rightarrow 0$. However this argument holds only when the Fenton-Wilkson approximation works within the order of accuracy of $O\Big(e^{\frac{\sigma^2_n}{2}}\Big)$, which is not always the case. If we know the exact distribution of $Y$, this argument can be modified to see if such an argument holds in general.

[1] Cobb, Barry R., R. Rumi, and Antonio Salmerón. "Approximating the distribution of a sum of log-normal random variables." Statistics and Computing 16.3 (2012): 293-308. http://leo.ugr.es/pgm2012/submissions/pgm2012_submission_6.pdf

• Is the error Fenton-Wilkson approximation proved or a numerical empirical result tested for some large number of examples? Would you mind specifying the definition of the "accuracy" that you say the Fenton-Wilkson works within? Thank you for your answer, Henry L. – Hans Jun 6 '17 at 18:22
• @Hans AFAIK, it is a numerical empirical result tested for a wide range of applications, but it may fail, you should look for expertise on the topic. For the accuracy I think they match the first two moments, so it depends on how good a second order approximation is around the point you are to estimate. – Henry.L Jun 6 '17 at 20:52
• For comparing mean and median via moment-matching, I would want an approximation that matched the first three moments. After all, if we matched the first two moments with a normal distribution, we would get the wrong answer to this problem. – Matt F. Jun 6 '17 at 21:11
• @MattF. Sorry but why first 3 moments? Not sensical to me. – Henry.L Jun 6 '17 at 22:25
• Henry.L: I think @MattF. was suggesting that the third moment could point towards skewness in a way the first two might not. Stating that a log-normal approximation to the sum of i.i.d. log-normals is better than a CLT-type normal approximation to the sum for finite $n$ needs some justification which could come from this third moment. Otherwise, as $n$ increases without limit the log-normal approximation will converge in distribution to the normal approximation (in a standardised sense). – Henry Jun 21 '17 at 17:53

Here is a positive result for the average of two iid lognormals. Recall that the average of a lognormal distributed as $LN(\mu,\sigma)$ is $\exp(\mu+\sigma^2/2)$.

Let $y_2$ be the average of two iid lognormals each distributed as $LN(\mu,\sigma)$. Then:

For small $\sigma$, $Med[y_2] \approx \exp(\mu+\sigma^2/4)$, which is less than $E[y_2]$.

For large $\sigma$, $Med[y_2] \approx \exp(\mu+.545\sigma)$, which is less than $E[y_2]$.

From those we can establish that $Med[y_2] < E[y_2]$ for all positive $\sigma$.

The approximations are valid to second order in $\sigma$, in the sense that

$$P\left[y_2 < \exp(\sigma^2/4)\right] = \frac{1}{2} + O(\sigma^3)$$

$$P\left[y_2 < \exp(0.545\sigma)\right] = \frac{1}{2} + O(\frac{1}{\sigma^3})$$

where .545 is an abbreviation for $\ln(F_1^{-1}(\frac{1}{\sqrt{2}}))$, and $F_{\sigma}$ is the CDF of $LN[0,\sigma]$, and we take $\mu=0$ for simplicity.

Proof:

We prove the first approximation from the general convolution formula

\begin{align} P[y_2 < m] &= \int\int_{\frac{x_1+x_2}2<m}dF(x_1)dF(x_2)\ \ \text{for arbitrary i.i.d. with CDF }F \\ &= \int_{x_1} F(2m-x_1)\, dF(x_1)\ \ \text{by integrating over }x_2 \\ &= \int_{0}^{1} F\big(2m - F^{-1}(p)\big)dp\ \ \text{where } p=F(x_1) \\ \end{align}

which in this particular case gives $$P[y_2 < m] = \int_{0}^{F_\sigma(2m)} F_\sigma\big(2m - F_\sigma^{-1}(p)\big)dp.$$

$F_\sigma(2m)$ is the largest possible $p$ such that $(x_1+x_2)/2<m$ for some positive $x_2$.

The graph shows the line $1-p$ and the integrand for $\sigma=1/2, \ m=\exp(\sigma^2/4)$. The integrand approaches the line as $\sigma$ goes to 0. The median is $m$ iff the areas under the two curves are equal.

We use Mathematica to differentiate the formula for $P[y_2 < m]$ under the integral sign and at the limit of integration, and this computes that the coefficients in the power series are $1/2, 0, 0$.

F[u_] = Simplify[CDF[LogNormalDistribution[0, sigma], u], Assumptions -> u > 0]
x1 = u /. Solve[p == F[u], u][[1]]
x2 = u /. Solve[q == F[u], u][[1]]
qsol = q /. Solve[(x1 + x2)/2 == Exp[sigma^2/4], q][[1]]
plim = F[2 Exp[sigma^2/4]]
Limit[plim, sigma -> 0]
Limit[D[plim, sigma], sigma -> 0]
Limit[D[plim, {sigma, 2}], sigma -> 0]
Integrate[Limit[qsol, sigma -> 0], {p, 0, 1}]
Integrate[Limit[D[qsol, sigma], sigma -> 0], {p, 0, 1}]
Integrate[Limit[D[qsol, {sigma, 2}], sigma -> 0], {p, 0, 1}]
Plot[{qsol /. sigma -> 1/2, 1 - p, 1}, {p, 0, 1}, PlotRange -> {0, 1},
AxesLabel -> {"p=F[x1]", "q=F[x2]"}, PlotStyle -> {, , {Black, Thin}}]


We prove the second approximation from $$\frac{1}{2} < P{\Large[}y_2 < e^{0.545\sigma}{\Large]} < F_\sigma(2 e^{0.545\sigma})^2\ \text{ and } \lim_{\sigma\rightarrow\infty}F_\sigma(2 e^{0.545\sigma})^2=\frac{1}{2}$$

and the fact (not proved here) that $P[y_2 < e^{0.545\sigma}]$ has eventually monotonic derivatives.

The inequality comes from \begin{align} \frac{1}{2} &= F_\sigma(e^{0.545\sigma})^2 \\ &= P[x_1 < e^{0.545\sigma}\ \& \ x_2 < e^{0.545\sigma}] \\ &< P[\ \ \ \ \ \ \ x_1 + x_2 < 2e^{0.545\sigma} \ \ \ \ \ \ \ \ \, ] \\ &< P[x_1 < 2e^{0.545\sigma}\ \& \ x_2 < 2e^{0.545\sigma}] \\ &= F_\sigma(2 e^{0.545\sigma})^2 \\ \end{align}

We use Mathematica to show that the limit of the right-hand side is 1/2.

x3 = u /. Solve[1/Sqrt[2] == F[u], u][[1]]
Limit[F[2 x3]^2, sigma -> Infinity] // FullSimplify
N[x3]

• this is no more than a simulation result.... – Henry.L Jun 25 '17 at 14:38
• Here is my derivation for your formula. You can see if you would like to incorporate it into your answer. \begin{align}P[y_2 < m] &= \int_{y_2=\frac{x_1+x_2}2<m}dF_2(x_1,x_2) \quad\text{for arbitrary joint CDF }F_2 \\ &= \int\int_{\frac{x_1+x_2}2<m}dF(x_1)dF(x_2)\quad\text{for arbitrary i.i.d. with CDF }F \\ &= \int_{x_1} F(2m-x_1)\, dF(x_1)\quad\text{integrate over }x_2 \\ &= \int_{0}^{F_\sigma(2m)} F_\sigma\big(2m - F_\sigma^{-1}(p)\big)dp, \end{align} where $p=F(x_1)$ and $F = F_\sigma$ thus $x_1\in[0,2m]$ for a CDF $F_\sigma$ on a positive random variable $x_1$. – Hans Jun 27 '17 at 8:21