Statistical Efficiency of Median

I am interested in the method / formula which computes the values of the Relative Statistical Efficiency for sample sizes 2 to 16. That is, as to how we get the following values: "1.000" "0.743" "0.838" "0.697" "0.776" "0.679" "0.743" "0.669" "0.723" "0.663" "0.709" "0.659" "0.699" "0.656" "0.692" as shown in the graph from link 1.

There are two posts on the subject:

Asymptotic relative efficiency of median vs mean for Student t distribution

What is statistical efficiency of the median?

Your numbers are for Normal distributions. Other distributions yield other relative efficiencies. I will explain how to find these values -- but ultimately, for Normal distributions, you have to perform numerical integration.

Setting and definitions

We are contemplating taking a sample. This is modeled as $$n$$ identically distributed, independent random variables with common distribution $$F.$$ With this sample we will estimate the central tendency of $$F.$$ Two estimators are under consideration:

1. The sample mean is a random variable with variance $$\sigma^2/n$$ where $$\sigma^2$$ is the variance of $$F.$$ This follows from basic properties of the variance.

2. The sample median also is a random variable. When $$F$$ is a continuous distribution with density $$f$$ and $$n=2k-1$$ is an odd number, the density of the sample median is that of the $$k^\text{th}$$ order statistic given by $$f_{k}(x) = \binom{n}{k-1,1,n-k} F^{k-1}(x)\, f(x)\, (1-F(x))^{n-k}$$ The variance of the sample median is $$\mu_{2} - \mu_{1}^2$$ where, as usual, the moments are given as the integrals $$\mu_{m} = \int_{\mathbb{R}} x^m\, f_{k}(x)\,\mathrm{d}x.\tag{*}$$

The efficiency of an estimator is the reciprocal of its sampling variance. The relative efficiency of two estimators (hopefully of the same quantity, for otherwise this makes little sense) is the ratio of their efficiencies.

Why is this important? Because the variance of an estimator typically is directly proportional to the size of a sample needed to detect a given effect with a given confidence and power. Relative efficiency, therefore, is the relative cost of your sample when you use one estimator instead of another. For instance, when $$F$$ is Normal and $$n$$ is large, the sample mean's efficiency is approximately $$2/\pi\approx 0.637$$ relative to the sample median for estimating the central tendency. Thus, when your objective is to estimate this tendency precisely and you are certain $$F$$ really has a Normal distribution, then by using the sample mean you will need only $$63.7\%$$ as large a sample than by using the sample median for your estimator.

Application

When $$F$$ is a Normal distribution and $$n$$ exceeds $$1,$$ computing the moments $$(*)$$ requires numerical integration -- but this isn't hard to do, because $$F$$ is well understood and can be quickly and accurately computed. Moreover, we may assume with no loss of generality that $$F$$ is the standard Normal distribution because all Normal distributions are just $$F$$ but with a different unit of measurement for the variable and changing the unit won't change the relative efficiency.

To illustrate, I wrote code to perform the integrals in $$(*)$$ and plotted the relative efficiency of the sample mean compared to the sample median for three distributions: Normal (which could be any Normal distribution), Uniform (which could be any uniform distribution on an interval), and Student t with $$3$$ degrees of freedom.

(The calculations took no measurable time.)

• The values in the "Normal" plot at left are every other value quoted in the question, because I have analyzed the situation only for odd values of $$n.$$ Namely, these values are (from left to right) 1, 0.743, 0.697, 0.679, 0.669, 0.663, 0.659, 0.656.

• Notice how the sample mean is even more (relatively) efficient for the Uniform distribution (shown in the middle plot). This is because outliers cannot occur: uniform distributions are bounded between two values and samples tend to be spread evenly between those bounds.

• The sample mean is worse than the sample median for the Student t distribution (notice the values on the vertical axis in the right hand plot). This distribution has "heavy" tails that frequently produce huge outliers. The sample mean is sensitive to these outliers, much more so than the median.

This illustrates the risk you take in using the sample mean when you think you might have a Normal distribution but really you have something like a Student t distribution governing your data. It explains why the median might be preferred as an estimator even though (under ideal assumptions) it appears to be more efficient.

Remarks

Why focus on odd $$n$$? Only because the median for even $$n$$ is not uniquely defined: there is no definite middle value in an even sample. The usual rule of averaging the two values closest to the middle results in a more complicated expression than $$(*)$$ and I just didn't want to deal with it here, since it adds little insight.

Why perform numerical integration? Another solution is via simulation: generate a gazillion samples, compute the mean and median of each one, and compare the variances of those gazillion statistics. Although that's easily done, (1) the formula $$(*)$$ is amenable to mathematical analysis, which can reveal far more than any bunch of simulations, and (2) the numerical integrals can be computed far more quickly and with far greater accuracy than you will ever obtain through simulation.

Here is the R code that produced the figure. It was written to apply to any distribution you want--but it includes no safety checks to make sure the integration succeeded, so be careful!

#
# Use numerical integration to compute moments of the kth order statistic from a
# distribution with log cdf f, log pdf ff, and logsurvival function sf.
# Provide hints about the domain of integration (if possible) by supplying
# values for the limits lower and upper.
#
moment.order <- function(n, k, f, ff, sf, moments, lower, upper, ...) {
if (missing(moments)) moments <- 1
if (missing(lower)) lower <- -Inf
if (missing(upper)) upper <- Inf

cnst <- lfactorial(n) - lfactorial(k-1) - lfactorial(n-k)
g <- function(x) exp(cnst + f(x) + (k-1) * ff(x) + (n-k) * sf(x))

sapply(moments, function(m) {
integrate(function(x) x^m * g(x), lower, upper, ...)$$value }) } # # Create a list of distributions to supply to moment.order. # NB: Rather than asking the code to compute the variance of each distribution, # we simply stipulate it. # distros <- list( Normal = list(f = function(x) dnorm(x, log = TRUE), ff = function(x) pnorm(x, log = TRUE), sf = function(x) pnorm(x, log = TRUE, lower.tail = FALSE), lower = -10, upper = 10, variance = 1), Uniform = list(f = function(x) dunif(x, log = TRUE), ff = function(x) punif(x, log = TRUE), sf = function(x) punif(x, log = TRUE, lower.tail = FALSE), lower = 0, upper = 1, variance = 1/12), Student (3) = list(f = function(x) dt(x, 3, log = TRUE), ff = function(x) pt(x, 3, log = TRUE), sf = function(x) pt(x, 3, log = TRUE, lower.tail = FALSE), lower = qt(1e-4, 3), upper = qt(1-1e-4, 3), variance = 3) ) # # Process each of the distributions for a sequence of sample sizes n and a # given quantile q (which will be estimated by an order statistic). # n <- seq(1, 16, by=2) # Sample sizes q <- 0.50 # Quantile par(mfrow=c(1, length(distros))) for (s in names(distros)) { # # Compute the variance of the sample median. # NB: the subdivisions argument is passed to integrate to make it do its # work a little more precisely and accurately. # v <- with(distros[[s]], { sapply(n, function(n) { k <- ceiling(n * q) mu <- moment.order(n, k, f, ff, sf, 1:2, lower, upper, subdivisions=1e3) mu[2] - mu[1]^2 }) }) # # Compute the variance of the sample mean, using the variance of the # distribution as supplied by the caller. # phi <- distros[[s]]$$variance / n
#
# Plot relative efficiency.
#
plot(n, phi / v, pch=21, bg="Gray",
ylab="Relative Efficiency",
main=s)
}
par(mfrow=c(1,1))

• Excellent. One question, it is printing values in step size of 2. As a result the output is "1.000," 0.743," "0.697," "0.679" and so on for Normal. How do I get the full list of values on the Graph as shown on link 1. That is, "1.000" "0.743" "0.838" "0.697" "0.776" "0.679" "0.743" "0.669" "0.723" "0.663" "0.709" "0.659" "0.699" "0.656" "0.692". May 28 '20 at 9:30
• In the below link, last para, values by Dixon and Massey (1969) table is the same as "1.000" "0.743" "0.838" "0.697" "0.776" "0.679" "0.743" "0.669" "0.723" "0.663" "0.709" "0.659" "0.699" "0.656" "0.692". The graph output in Link 1. (books.google.co.in/…) May 28 '20 at 9:34
• As I tried to explain, the values for even $n=2k$ are based on the arithmetic means of the $k^\text{th}$ and $(k+1)^\text{st}$ order statistics. This requires a similar but more complicated calculation (there are two variables, so it's a double integral). It doesn't require any additional concepts and reveals little more about the overall behavior of the relative efficiencies.
– whuber
May 28 '20 at 13:51
• If you can share a basic R-code framework, for even n=2k, I will attempt it my end. Thanks, for the help. May 29 '20 at 18:18

I would suggest the following methodology using Monte Carlo methods, in a spreadsheet (note, Google produces free spreadsheet access, at least for their product line).

1. Use an inversion formula based on the quantile function to translate a uniformly simulated random deviate into a distribution of choice.

2. For a sample of size n, compute the observed absolute relative accuracy of the mean, median,..

3. Tabulate the statistics for each run from Step 2.

4. Repeat until you have obtained the desired level of accuracy or run into other constraints.