Variance of sample median of normal distribution

Question

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.

Answer 1

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.

Answer 2

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).

Answer 3

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

Start with the simplest case, that of iid uniform random variables.

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.