What are the properties of the distribution of medians from samples of a set of integers?

Given a set of integers $$S = \{1, 2, 3, ..., n\}$$ I draw uniformly at random without repetition $$m$$ integers a large number of times ($$1 \leq m \leq n$$). For each of these samples (each of size $$m$$) I compute the median. What is the expected distribution of these medians?

By intuition, we can say the sample space from which we can take the samples is the set of all permutations without repetition of the elements of $$S$$ of size $$m$$. Different permutations of this set will have the same median making certain median values more likely than others so the resulting distribution is not uniform for $$1 \lt m \lt n$$. For the extreme cases $$m = 1$$ we should observe a uniform distribution in the range $$[1, n]$$ and for $$m = n$$ all medians will be the median of $$S$$. I am unsure on what to state about the distribution of medians for all other values though.

• If you make $m$ odd you can keep all the medians as integers Dec 1 '21 at 2:31
• Building off the last comment, we can explicitly count the number of combinations leading to a median of k. These combinations need to contain k, and the contain (m-1)/2 elements < k and (m-1)/2 elements > k. How many combinations is this? Divide that number by the total number of combinations Dec 1 '21 at 2:36
• The formula gets a bit more complicated when m is even. In that case, one way to split up the world that may be fruitful: how many ways can you get k such that the middle two numbers are precisely x<y such that (x+y)/2 = k? Then add up for all x,y pairs satisfying this property Dec 1 '21 at 2:39
• @stats_model Your first comment suggests a combinatorial identity I had not seen before: $\sum\limits_{j=k}^{n-k} {j \choose k}{n-j \choose k} = {n+1 \choose 2k+1}$ Dec 1 '21 at 2:55

Let $$X_1,...,X_m \sim \text{SRSWOR}(n)$$ denote the sampled values. If you have an odd number of values in the sample then the median is just the middle order statistic. I am going to generalise from this and show how you can obtain the distribution of any order statistic $$X_{(k)}$$ over the order values $$k =1,...,m$$. (Getting the median in the case where you have an even number of values is a bit trickier, since it involves the joint distribution of two adjacent order statistics; I will leave this as an exercise.)

There are a couple of different ways you can derive the distribution for this order statistic. The simplest method is through a direct combinatorial argument. There are $${n \choose m}$$ possible samples that can be drawn from the population (with equal probability). To get the event $$X_{(k)} = x$$ we have $${x-1 \choose k-1}$$ ways of choosing the previous order statistics and $${n-x \choose m-k}$$ ways of choosing the subsequent order statistics. Consequently, we have:

$$\mathbb{P}(X_{(k)} = x) = \frac{{x-1 \choose k-1} {n-x \choose m-k}}{{n \choose m}} \quad \quad \quad \quad \quad \text{for } x = 1,...,n.$$

(To confirm that this is a valid probability mass function ---i.e., that it sums to one over the support--- note this related summation result for products of binomial coefficients.) It is also possible to obtain this result by using the relationship of the order statistics to the hypergeometric distribution, but this latter method is much more complicated (omitted here).

With a bit of algebra, it is possible to derive the general form for the factorial moments of this distribution, and through this obtain the mean and variance. These are given by:

$$\mathbb{E}(X_{(k)}) = \frac{n+1}{m+1} \cdot k \quad \quad \quad \quad \quad \mathbb{V}(X_{(k)}) = \frac{(n+1)(n-m)}{(m+1)^2 (m+2)} \cdot k(m-k+1).$$

It is simple to program this mass function in R using standard syntax for density functions. I will refer to this distribution as the finite-population-order-statistic (FPOS) distribution and program its density function as dFPOS. Appealing to standard notation in these problems, I will use N for the population size and n for the sample size. The density function is programmed using computation in log-space, in order to ensure computational stability.

#Probability mass function for FPOS distribution
dFPOS <- function(x, k, n, N, log = FALSE) {

#Check input types
if (!is.numeric(x))    stop('Input x must be numeric')
if (!is.numeric(k))    stop('Order k must be a positive integer')
if (!is.numeric(n))    stop('Sample size n must be a positive integer')
if (!is.numeric(N))    stop('Population size N must be a positive integer')
if (!is.logical(log))  stop('Input log must be a logical value')
ki <- as.integer(k)
ni <- as.integer(n)
Ni <- as.integer(N)
if (k != ki)           stop('Order k must be a positive integer')
if (n != ni)           stop('Sample size n must be a positive integer')
if (N != Ni)           stop('Population size N must be a positive integer')

#Check input lengths and values
if (length(k) != 1)    stop('Order k should be a single positive integer')
if (length(n) != 1)    stop('Sample size n should be a single positive integer')
if (length(N) != 1)    stop('Population size N should be a single positive integer')
if (length(log) != 1)  stop('Input log should be a single logical value')
if (k < 1)             stop('Order k must be positive')
if (n < 1)             stop('Sample size n must be positive')
if (N < 1)             stop('Population size N must be positive')
if (k > n)             stop('Order k cannot be greater than sample size n')
if (n > N)             stop('Sample size n cannot be greater than population size N')

#Compute log-probabilities
LOGPROBS <- rep(0, N)
for (t in 1:N) {
LOGPROBS[t] <- lchoose(t-1, k-1) + lchoose(N-t, n-k) - lchoose(N, n) }

#Generate output
L <- length(x)
OUT <- rep(-Inf, L)
for (i in 1:L) {
xx <- x[i]
if (xx %in% (1:N)) {
OUT[i] <- LOGPROBS[xx] } }

#Return output
if (log) { OUT } else { exp(OUT) } }

Using this function we can obtain the distribution of the median (for an odd sample size) by taking $$m=2k-1$$ and finding the distribution of the order statistic $$X_{(k)}$$. In the code below we compute and plot the mass function for the median for the case where $$n=16$$ and $$m=9$$.

#Set parameters
n <- 16
m <- 9

#Compute and plot the mass function
PROBS <- dFPOS(1:n, k=(m+1)/2, N = n, n = m)
names(PROBS) <- 1:n
barplot(PROBS, xlab = 'Value', ylab = 'Probability', col = 'blue') • Although your answer is a bit hard to follow for me (have no formal background on maths or stats), it is useful to get an idea of what is going on. The intuition I get is that the probability of observing any given median less than or equal to $x$ depends on the probability of getting $k + 1$ values greater than or equal to $x$ in $m$ tries (which is why the hypergeometric is used). Then to compute an exact probability for a specific $x$, you just compute the cumulative mass probability of $x$ and $x - 1$ then substract them to get the probability "in the middle". Am I getting it more or less? Dec 1 '21 at 8:36
• I think you are getting it more or less, though the probability of an individual value is not really "in the middle". Rather, you have $\mathbb{P}(\text{median} = x) = \mathbb{P}(\text{median} \leqslant x) - \mathbb{P}(\text{median} \leqslant x-1)$. In any case, the above R function should allow you to compute the distribution for any parameters you wish (within the computational ability of the machine you are using).
– Ben
Dec 1 '21 at 8:56

First, consider the situation where the sample size $$m=2k+1$$ is odd so that the median is uniquely defined.

For $$x\in S,$$ the chance the median $$M$$ equals $$x$$ is the chance that at $$k$$ of the sample values are from the set $$\{1,2,\ldots, x-1\}$$ and another $$k$$ of the sample values are from $$\{x+1,\ldots, n\},$$ with the remaining value equal to $$x$$ itself. Such samples are in one-to-one correspondence with ordered pairs $$(I,J)$$ where $$I$$ is a $$k$$-subset of $$x-1$$ values and $$J$$ is a $$k$$-subset of $$n-x$$ values. Since there are $$\binom{n}{2k+1}$$ equiprobable samples, this chance therefore is

$$\Pr(M=x) = \frac{\binom{x-1}{k}\binom{n-x}{k}}{\binom{n}{m}} ;\ x \in S, m=2k+1.$$

When $$m=2k$$ is even, let the median be the midpoint of the interval that splits the sample into two equal halves. When $$2 \le 2x \le 2n$$ is an integer, the median equals $$x$$ when the two closest numbers to $$x$$ in the sample are equally far from $$x$$--let's call them $$i\lt x$$ and (therefore) the other is $$2x-i \gt x$$--and the remaining $$2(k-1)$$ elements of the sample are split into two equal-size groups.

To restate this a little more precisely, there must be a (unique) integer $$i \lt x$$ for which (a) $$i$$ and $$2x-i$$ are in the sample, (b) exactly $$k-1$$ elements of the sample are in the set $$\{1,2,\ldots, i-1\}$$ (implying $$i\ge k$$) and (c) exactly $$k-1$$ elements are in the set $$\{2x-i+1, 2x-i+2, \ldots, n\}$$ (a set with $$n - (2x-i)$$ elements). Because $$i$$ is unique, the events so described are disjoint: this implies their probabilities sum to the chance that $$x$$ is the median. Thus

$$\Pr(M = x) = \frac{1}{\binom{n}{m}}\sum_{i=k}^{\lfloor x\rfloor} \binom{i-1}{k-1}\binom{n-2x+i}{k-1};\ m=2k, 2x\le n+1.$$

Because the distribution is symmetric about $$(n+1)/2,$$ use $$n+1-x$$ in place of $$x$$ in this formula when $$2x \gt n+1.$$

Here are some plots of the median distributions for various sample sizes from $$n=39$$ values: The alternating probabilities of the half-integral and integral medians for even sample sizes $$m$$ are apparent.

Simulations confirm the correctness of these formulas using a chi-squared test.

n <- 7
N <- 5e2
for (m in 1:n) {
X.2 <- replicate(N, 2*median(sample.int(n, m)))
p <- pmed(seq(1/2, n, by=1/2), n, m)
i <- p > 0
observed <- tabulate(X.2, nbins=2*n)[i]
p <- p[i]
Y <- rbind(observed, expected = N*p)
colnames(Y) <- which(i)/2

cat("m =", m, "; n =", n, "\n")
print(signif(Y, 3))

if (sum(i) > 1) print(chisq.test(observed, p=p))
}

Typical output tallies the simulated medians, compares them to the formulas, and displays the results of the chi-squared test.

m = 4 ; n = 7
2.5    3 3.5    4 4.5    5  5.5
observed 69.0 34.0 107 67.0 118 45.0 60.0
expected 57.1 42.9 114 71.4 114 42.9 57.1

Chi-squared test for given probabilities

data:  observed
X-squared = 5.4006, df = 6, p-value = 0.4936

It performs the calculations in R using this quick and dirty implementation of the two formulas. Using logarithms of the binomial coefficients makes the calculation feasible for very large $$n$$ (millions or more), whereas computing with the binomial coefficients themselves would overflow double-precision floats once $$n$$ exceeds $$1000$$ or so.

pmed <- Vectorize(function(x, n, m) {
k <- floor(m/2)
if (!isTRUE(m == 2*k)) {
if (x == round(x)) {
exp(lchoose(x-1, k) + lchoose(n-x, k) - lchoose(n, m))
} else {
0
}
} else {
x <- pmin(x, n+1-x)
x.0 <- floor(x-1/2)
if(x.0 >= k) {
i <- seq(k-1, x.0-1, by=1)
q <- lchoose(i, k-1) + lchoose(n+1-2*x+i, k-1)
sum(exp(q - lchoose(n, m)))
} else {
0
}
}
}, "x")