Variance of sample median of normal distribution

Let $$\mathbf{R} \ni X_1, \dots, X_n \sim \mathcal{N}(0, \sigma^2)$$. I can show that $$\mathbb{E} ~ \text{Med} \{X_1, \dots, X_n\} = 0$$ and want to compute the variance of the same sample median, i.e. $$\mathbb{V}ar ~ \text{Med} \{X_1, \dots, X_n\}$$

I have a guess that it should behave like sample mean ($$\asymp \frac {\sigma^2} n$$), since everything is symmetric and stuff, but no idea how to show it rigorously.

• How do you define sample median for an even-sized sample? Oct 31 '19 at 15:27
• The Theorem referred to here is Thm 4.1 in Statistics and Data Analysis for Financial Engineering and gives you the $/n$ relationship. Oct 31 '19 at 18:40

For definiteness, consider odds $$n$$, say $$2m+1$$.

This is an additional comment more than an answer (the second comment by @LinAlg refers to a Theorem that gives the complete answer anyway), to grasp intuition about about the 1/n factor

Then the ordered statistics of such random variables are well known to be beta random variables, and the median itself will be Beta($$m$$, $$m+1$$) if I am not mistaken, the variance of which (check wikipedia) is

$$\frac{m(m+1)}{(2m+1)^2 \times 2(m+1)}= \frac{m}{2(2m+1)^2} \sim \frac{1}{8m}$$

Now you can map your iid uniform to iid Gaussian using the inverse distribution function (in fact, since the Beta concentrates around 1/2, we will only need the derivative of the function at this point). This way, you can even work out the constant in front of $$1/m$$, but doing this rigorously is parhaps not that easy I assume.

This is not an answer, but it's too long for a comment.

Clearly, the answer must scale proportionally to $$\sigma^2$$, so let's just assume $$\sigma=1$$.

Sampling $$x_1, \ldots, x_n$$ i.i.d. from a standard Gaussian is equivalent to sampling a random Gaussian vector $$\vec{x} \in \mathbb{R}^n$$ with mean $$\vec{0}$$ and identity covariance matrix. The median coordinate $$\mathrm{Med}(\vec{x})$$ satisfies $$\mathrm{Med}(\vec{x}) = ||\vec{x}|| \mathrm{Med}(\hat{x})$$ where $$\hat{x} = \frac{\vec{x}}{||\vec{x}||}$$ is the corresponding unit vector.

For Gaussians, the radius and unit vector are independent, so the calculation factors: $$\mathbb{E}_{\vec{x} \sim N(\vec{0}, I)}[\mathrm{Med}(\vec{x})^2] = \mathbb{E}_{r,\ \hat{x} \sim \mathrm{unif}(S^{n-1})}[r^2 \mathrm{Med}(\hat{x})] = \mathbb{E}_r[r^2] \mathbb{E}_{\hat{x} \sim \mathrm{unif}(S^{n-1})}[\mathrm{Med}(\hat{x})^2].$$

Here $$r$$ is just the norm of $$\vec{x}$$ so $$\mathbb{E}_r[r^2] = \mathbb{E}_{\vec{x}}[||\vec{x}||^2] = n$$. The distribution of $$\hat{x}$$ is uniform on the sphere by symmetry.

I don't know how to calculate this integral over the sphere though. You can do obvious things like dividing by $$n!$$ and restricting to the sector $$U_0 = \{x_1 \leq x_2 \leq \cdots \leq x_n\} \subset S^{n-1}$$: $$\mathbb{E}_{\hat{x} \sim \mathrm{unif}(S^{n-1})}[\mathrm{Med}(\hat{x})^2] = n! \int_{\hat{x} \in U_0} (x_{(n+1)/2})^2 d\mu(\hat{x}),$$ but already for $$n=5$$ Mathematica doesn't manage to compute the integral.

WLOG, assume that $$\sigma = 1$$.

We have a random sample $$X_1, X_2, \cdots, X_n$$ of size $$n$$ from $$X\sim N(0, 1)$$. The probability density function of $$X$$ is $$f(x) = \frac{1}{\sqrt{2\pi}}\mathrm{exp}(-\frac{x^2}{2})$$. The cumulative distribution function of $$X$$ is $$\Phi(x) = \int_{-\infty}^x \frac{1}{\sqrt{2\pi}}\mathrm{exp}(-\frac{t^2}{2}) dt$$.

The order statistics $$X_{(1)} \le X_{(2)} \le \cdots \le X_{(n)}$$ are obtained by ordering the sample $$X_1, X_2, \cdots, X_n$$ in ascending order.

The probability density function of the $$k$$-th order statistic $$X_{(k)}$$ is $$f_k(x) = \frac{n!}{(k-1)!(n-k)!}[\Phi(x)]^{k-1}[1-\Phi(x)]^{n-k}f(x), \quad -\infty < x < \infty.$$ The joint probability density function of $$k$$-th and $$(k+1)$$-th order statistics $$X_{(k)}$$ and $$X_{(k+1)}$$ is $$f_{k, k+1}(x, y) = \frac{n!}{(k-1)!(n-k-1)!}[\Phi(x)]^{k-1}[1-\Phi(y)]^{n-k-1} f(x) f(y), \quad x \le y.$$

If $$n$$ is odd, we have $$\mathrm{Med}(X_1, X_2, \cdots, X_n) = X_{(n+1)/2}$$ and hence \begin{align} &\mathrm{E}[\mathrm{Med}(X_1, X_2, \cdots, X_n)^2]\\ =\ & \int_{-\infty}^\infty x^2 f_{(n+1)/2}(x) dx \\ =\ & \int_{-\infty}^\infty x^2 \frac{n!}{(\frac{n-1}{2})!^2}[\Phi(x)- \Phi(x)^2]^{(n-1)/2}\frac{1}{\sqrt{2\pi}}\mathrm{exp}(-\frac{x^2}{2}) dx. \end{align}

If $$n$$ is even, we have $$\mathrm{Med}(X_1, X_2, \cdots, X_n) = \frac{1}{2}(X_{n/2} + X_{n/2+1})$$ and hence \begin{align} &\mathrm{E}[\mathrm{Med}(X_1, X_2, \cdots, X_n)^2]\\ =\ & \int_{-\infty}^\infty \int_{-\infty}^\infty \frac{1}{4}(x+y)^2 f_{n/2, n/2+1}(x, y)\ 1_{x < y}\ dx dy\\ =\ & \int_{-\infty}^\infty \int_{-\infty}^\infty \frac{1}{4}(x+y)^2 \frac{n!}{(n/2-1)!^2}[\Phi(x)]^{n/2-1}[1-\Phi(y)]^{n/2-1} f(x) f(y) 1_{x < y} dx dy \end{align} where $$1_{x < y}$$ is the indicator function.

For both cases, since $$\mathrm{E}[\mathrm{Med}(X_1, X_2, \cdots, X_n)] = 0$$, we obtain $$\mathrm{Var}[\mathrm{Med}(X_1, X_2, \cdots, X_n)] = \mathrm{E}[\mathrm{Med}(X_1, X_2, \cdots, X_n)^2].$$

Numerically verified: I used Maple software to calculate the integrations. I also used Matlab to do Monte Carlo simulation (i.e., generate many group of normal distribution data, calculate median of each group, average them).