# Moments of the sample median of a normal distribution

Given a sample $X_1,X_2,X_3,........,X_n$ from a normal distribution, that is $X_i$~$N(\mu,\sigma^2)$, the objective is to find the moments or as of now, the expectation of the sample median say $Y_{med}$. Assuming that the number of observations is odd i.e. $n->2n+1$.

By formulating the distribution of the $k\:th$ order statistic, the distribution of the median is given by: $$f_{Y_{med}}(y)=\dfrac{(2n+1)!}{(n!)^2}[F_{X_i}(y)(1-F_{X_i}(y))]^{n}f_{X_i}(y)$$

where $f_{X_i}(y)$ is a probability distribution function of a normal random variable and $F_{X_i}(y)$ is the cumulative distribution function.

The only hint I came across is that it seems like a beta distribution as $0\leq F_{X_i}(y)\leq1$, also $\dfrac{(2n+1)!}{(n!)^2}$ seems like $\dfrac{\tau(2n+2)}{\tau(n+1)\tau(n+1)}$.

Proving $\int f_{Y_{med}}(y)=1$ is relatively easy, taking $F_{X_i}(y)=t$ the integral does convert into a beta distribution, but what about $E(Y_{med})$, I am really concerned about the extra variable $y$ involved: $\int \dfrac{(2n+1)!}{(n!)^2}y[F_{X_i}(y)(1-F_{X_i}(y))]^{n}f_{X_i}(y)dy$.

Can someone help ?

• Approximations are easy to come by and numerical integration is straightforward. If you need theoretical expressions for the moments, you might be disappointed. – whuber Sep 14 '17 at 16:18
• One can show, for a standard Normal parent, that all the (positive) odd order moments of the sample median are zero. For a general Normal parent, that should generalise at least to E[sample median] = $\mu$ – wolfies Sep 14 '17 at 16:28
• So there's no way to tackle $\int \dfrac{(2n+1)!}{(n!)^2}y[F_{X_i}(y)(1-F_{X_i}(y))]^{n}f_{X_i}(y)dy$ ? @whuber – User9523 Sep 14 '17 at 17:12
• Some of those moments can be found by symmetry arguments. For instance, when $\mu=0$ obviously the odd moments are zero. Otherwise, I believe these integrals cannot be expressed in closed form using any of the standard functions of math, stats, and physics. – whuber Sep 14 '17 at 17:27
• To OP, as per my comment above, one can formally show that the sample median from a $N(\mu,\sigma^2)$ parent is an unbiased estimator of $\mu$, i.e. that E[sample median] = $\mu$, given your odd-sized sample. – wolfies Sep 16 '17 at 17:23

The following is not an exact answer, but it does provide some help in the form of a heuristic for the central moment.

As mentioned in the comments by others, a theoretical expression is difficult (except for the odd moments which are zero), if not unobtainable. Yet, for the 2nd moment of the median of samples with 3 and 5 variables, it is possible to find expressions using integration by parts and tables of integrals.

# Experimental

I started plotting some computations in order to get an idea in what direction the integral might go. And the result is close to a power-law relationship in for the ratio of the moment of the Gaussian distribution and the moment of the median is of 2n+1 Gaussian distributed variables.

For the n-th central moment, $m$, of the Gaussian distribution we have:

$E((X-\mu_X)^m) = \begin{cases} \sigma^m(m-1)!! & \text{for$m$even,} \\ 0 & \text{for$m$odd.} \end{cases}$

with $n!!$ the double factorial.

For the n-th central moment of the median of a sample with 2n+1 Gaussian distributed variables, $X_{2n+1}$ we find roughly a powerlaw:

$\frac{E((X_{2n+1}-\mu_{X_{2n+1}})^m)}{E((X-\mu_X)^m)} \sim (2n+1)^{-m/2}$

In the plot with the computations I have added a line that is based on the following fit:

$\frac{E((X_{2n+1}-\mu_{X_{2n+1}})^m)}{E((X-\mu_X)^m)} \sim (1+0.6355323(n-1))^{-m/2}$

Note that the parameters in this fit only change a little bit if different moments are calculated.

Legend for the image below:

• 2-nd moment of the median as function of number of variables in the sample
• black dots are computed 2-nd central moments
• green crosses are theoretic calculations
• red line is a power law relationship based on an experimental linearized model fit

# Theoretic

The image above contains two points that have been calculated exactly. These points were found using integration by parts.

First we can eliminate this $\dfrac{(2n+1)!}{(n!)^2}$ part by using a ratio:

$\frac{\int \dfrac{(2n+1)!}{(n!)^2}y^m[F_{X_i}(y)(1-F_{X_i}(y))]^{n}f_{X_i}(y)dy}{\int \dfrac{(2n+1)!}{(n!)^2}[F_{X_i}(y)(1-F_{X_i}(y))]^{n}f_{X_i}(y)dy}$

The integration by parts follows a pattern like (I won't write down the entire procedure):

$\int y^a f^b (F-0.5)^c = \frac{(a-1)}{b}\int y^{a-2} f^b (F-0.5)^c + \frac{c}{b} \int y^{a-1} f^{b+1} (F-0.5)^{c-1}$

And in this way you eliminate the $y$-term such that you are left with only products of $f^b F^c$ in the integral. The cases $f^b$, $f F^c$, $f^b F$, are easy to calculate. However, if both $b>2$ and $c>2$, then it becomes difficult. For the case of 5 variables we do get an integral with a $f^3F^2$ term. This can be solved using the tables from Murray Geller 'A table of integrals of the Error functions' in which we find:

$\int_0^\infty \mathrm{erf}(ax) \, \mathrm{erf}(bx) \, e^{-c^2x^2} dx = \frac{1}{c\sqrt(\pi)} \mathrm{tan}^{-1} \left( \frac{ab}{c\sqrt{a^2+b^2+c^2}} \right)$

In cases of more than 5 variables we end up with larger power terms for which there are not standard integrals published. You will have to play around with more than simple integration by parts.

So we get for the cases $n=3$ and $n=5$:

$E((X_{3}-\mu_{X_{3}})^2) \sim E((X-\mu_X)^2) \frac{\frac{1}{6}-\frac{\sqrt{3}}{\sqrt{2 \pi}^2}}{\frac{1}{6}}$

and

$E((X_{5}-\mu_{X_{5}})^2) \sim E((X-\mu_X)^2) \frac{\frac{1}{30}-\frac{1}{4 \sqrt{3} \pi}+\frac{6}{\pi^2 4 \sqrt{3}} \mathrm{tan}^{-1}\left( \frac{1}{\sqrt{15}} \right)}{\frac{1}{30}}$

In which I excuse myself for not simplifying the expression and showing every detail of the derivation. I leave that open. The code below shows that these expressions work.

Code for the image:

#function for double factorial
odfactorial <- function(x) {
l <- (x+1)/2
return(factorial(2*l)/factorial(l)/2^l)
}

#define some stuff
mom = 2
sigm = 1
x = 10*c(-1000:1000)/1000 #x variable
k <- c(0:500) # loop over odd numbers of observations
n <- 1+2*k # number of observations

#calculate moments
moments_p <- lapply(k,
FUN <- function(kj) sum(x^mom*dnorm(x, 0, sigm)*(pnorm(x, 0, sigm) - pnorm(x, 0, sigm)^2)^kj)/
sum(      dnorm(x, 0, sigm)*(pnorm(x, 0, sigm) - pnorm(x, 0, sigm)^2)^kj))
moments_p <- as.numeric(moments_p) # convert array to vector

# plot expected moment versus number of observations
plot(n, moments_p, xlab="observations", ylab=paste0(mom, "-th moment"), log="xy")

#lines(n, sigm^(mom) * (odfactorial(mom-1))/n^(mom/2))

#draw heuristic line
m_t <- sigm^(mom) * (odfactorial(mom-1)) * (1+0.6355323*(n-1))^(-mom/2)
lines(n,m_t,col=2)

#draw calculated points
m1_predict = sigm^(mom)*(odfactorial(mom-1))
points(1,m1_predict,col=3,pch=3,cex=1.45)
m2_predict = sigm^(mom)*(odfactorial(mom-1))*(1/6-sqrt(1/3)/sqrt(2*pi)^2)/(1/6)
points(3,m2_predict,col=3,pch=3,cex=1.45)
m3_predict = sigm^(mom)*(odfactorial(mom-1))*((1/30)-0.5/sqrt(3)/2/pi+6*(atan(1/sqrt(3)/sqrt(5))/pi^2/sqrt(3)/4))/(1/30)
points(5,m3_predict,col=3,pch=3,cex=1.45)