Suppose I have a simple random sample drawn from a population with a known distribution on some population characteristic like height or income, with probability density function (pdf) $f(x)$. Order the sample on the measured value. How would I calculate the variance of the mean of a subset of the sample between the order statistics?

I would accept either an analytic or a numerical methods approach as an answer. I’ll accept a resampling approach if there is not a known analytic or numerical approach for an arbitrary known distribution.

I believe this to be a simplified mathematical statement of the problem that I posted previously motivated with a more narrative statement, here: Statistician's quantiles versus policymaker's quantiles: calculating the variance correctly

  • $\begingroup$ Could you please more clearly articulate what "the variance of the mean of a subset of the sample between the order statistics" might possibly be? Are you saying, from a sample $x$ and two prespecified ranks $j\le k,$ consider the sampling distribution of the empirical mean of the order statistics between those ranks? A formula is $$\bar x_{j;k} := \frac{1}{k-j+1}\sum_{i=j}^k x_{[i]}.$$ $\endgroup$
    – whuber
    Sep 1, 2019 at 16:28

2 Answers 2


This answer takes seriously the assumption that the data-generating distribution is known. Whether a closed-form analytical solution exists depends on the distribution, but it's possible to write down a universal formula for continuous distributions that permits numerical calculation of the variances of the quantile means.

Let's begin by formulating the problem. A random sample of a continuous distribution $F$ (with density $f$) is collected. When sorted, the values may be written as order statistics

$$x_{(1)} \le x_{(2)} \le \cdots \le x_{(n)}.$$

The background for this question concerns partitioning these data into "quantiles." For instance, we might take the first fifth of the order statistics from $x_{(1)}$ through $x_{(n/5)}$ (rounding $n/5$ suitably), the next fifth, and so on, creating five subsets; and characterize each such subset in terms of its mean.

Generally, given two indexes $1\le i \le j \le n,$ define the subsample

$$X_{i;j} = (x_{(i)}, x_{(i+1)}, \ldots, x_{(j)})$$

and let

$$\bar x_{i;j} = \frac{1}{j-i+1} (x_{(i)} + x_{(i+1)} + \cdots + x_{(j)})$$

be its mean. These are statistics--functions of the sample--and as such have a sampling distribution. The question I will address concerns finding the first and second moments of $(\bar x_{i;j}, \bar x_{i^\prime;j^\prime}).$ These could be used to compute the variance-covariance matrix of the quantile means, for instance.

The analysis is straightforward upon observing that the distribution of $x_{(i)}$ has density

$$f_i(x) = \binom{n}{i-1,1,n-i} F(x)^{i-1}(1-F(x))^{n-i} f(x)$$

and the joint distribution of $(x_{(i)}, x_{(j)})$ has density

$$\eqalign{&f_{i,j}(x,y) = \\&\binom{n}{i-1,1,j-i-1,1,n-j} F(x)^{i-1} (F(y)-F(x))^{j-i-1} (1-F(y))^{n-j} f(x) f(y) \mathcal{I}(x\le y).}$$

(See https://stats.stackexchange.com/a/4684/919 for an intuitive explanation of how such a formula can be derived and remembered.) The normalizing constants are multinomial coefficients and can be evaluated as

$$\binom{n}{a,b,\ldots, z} = \exp\left(\log \Gamma(n+1)-\log\Gamma(a+1)-\log\Gamma(b+1)-\cdots - \log\Gamma(z+1)\right).$$

The definitions of moments imply

$$\mu_{i;k} = E[x_{(i)}^k] = \int_\mathbb{R} x^k f_i(x) \mathrm{d}x$$

for the univariate moments of order $k$ and

$$\mu_{i,j;k,l} = E[x_{(i)}^k x_{(j)}^l] = \iint_{\mathbb{R}^2} x^k y^l f_{i,j}(x,y) \mathrm{d}x \mathrm{d}y$$

for the bivariate moments of order $(k,l).$ For $j\gt i$ the covariance is

$$\operatorname{Cov}(x_{(i)}, x_{(j)}) = \mu_{i,j;1,1} - \mu_{i;1}\mu_{j;1}$$

and for $j=i$ it is

$$\operatorname{Cov}(x_{(i)}, x_{(i)})=\operatorname{Var}(x_{(i)}) = \mu_{i,2} - \mu_{i;1}^2.$$

In some special cases--notably, when $F$ is uniform--these integrals can be evaluated, but usually they need numerical integration.

Let, then, $\Sigma_F$ denote the variance-covariance matrix of all the order statistics,

$$\Sigma_F(i,j) = \operatorname{Cov}(x_{(i)}, x_{(j)}).$$

The quantile means $\bar x_{i,j}$ are linear combinations of the order statistics with coefficients $1/(j-i+1)$ for $x_{(k)}$ when $i\le k \le j$ and zero otherwise. Let $\xi_{i,j}$ be this vector of coefficients:

$$\xi_{i,j}(k) = \left\{\matrix{\frac{1}{j-i+1}& i\le k \le j \\ 0& \text{otherwise.}}\right.$$

Because covariance is a quadratic form,

$$\operatorname{Cov}(\bar x_{i,j}, \bar x_{i^\prime, j^\prime}) = \xi_{i,j}^t\,\Sigma_F\,\xi_{i^\prime, j^\prime}.$$

The R code below carries out this analysis for any distribution for which you can compute the values of its CDF $F$ and PDF $f.$ It illustrates it with the standard Normal distribution, then checks its results with a simulation of 10,000 independent samples. The output of that check is a histogram of the relative errors of the $5\times 5 = 25$ entries in the variance-covariance matrix of the five quintiles. That they are all small--not much larger than $1/\sqrt{10000}=0.01$ in size--attests to the correctness of the results.

Histogram of relative errors shows they lie between -0.01 and +0.04.

For large values of $n,$ the computational demands will grow too large. Further analysis of the integrals for any particular $F$ will help develop approximations or asymptotic formulas. Saddle-point approximations may be particularly attractive.

# Log multinomial coefficients.
# The arguments are the *bottom* numbers: they must sum to the top value `n`.
lmultinom <- function(...) {
  i <- unlist(list(...))
  lgamma(sum(i)+1) - sum(lgamma(i+1))
# Univariate moments of order statistics.
# `lpf` returns the log of the CDF.  It must recognize a "lower.tail=FALSE"
# argument to return the log of the complementary CDF.
# `ldf` returns the log of the PDF.
# `k` is the order of the moment.
e1 <- function(i, n, lpf, ldf, lower, upper, k=1, ...) {
  lcnst <- lmultinom(i-1, 1, n-i)
  f <- function(x) {
   x^k * exp(lcnst + (i-1)*lpf(x) + (n-i)*lpf(x, lower.tail=FALSE) + ldf(x))
  integrate(f, lower, upper, ...)
# Bivariate moments of order statistcs.
# `k` and `l` are the orders of the moments.
e2 <- function(i, j, n, lpf, ldf, lower, upper, k=1, l=1, ...) {
  lcnst <- lmultinom(i-1, 1, j-i-1, 1, n-j)
  f <- function(x, y) {
    x^k * y^l * exp(lcnst + (i-1)*lpf(x) + (n-j)*lpf(y, lower.tail=FALSE) +
                  (j-i-1)*log(exp(lpf(y)) - exp(lpf(x))) + ldf(x) + ldf(y))
  g <- Vectorize(function(y) {
    integrate(f, lower, y, y=y, ...)$value
  integrate(g, lower, upper, ...)
# EXAMPLE: The standard Normal distribution.
lpf <- function(x, ...) pnorm(x, log=TRUE, ...)
ldf <- function(x, ...) dnorm(x, log=TRUE, ...)
lower <- -8; upper <- 8 # Effective range

n <- 20 # The double integrals take a few seconds, but only have to be 
        # ..performed once and can be stored for efficiency.

#---Compute the means of the order statistics.
y <- sapply(1:n, function(i) e1(i, n, lpf, ldf, lower, upper)$value)
# plot(1:n, y, 
#      sub=bquote(n==.(n)),
#      main="Expected Order Statistics\nStandard Normal Distribution",
#      xlab="Order", ylab="Value")

#---Compute the second moments of the order statistics.
z <- matrix(0.0, n, n)
for (i in seq_len(n)) {
  for (j in seq_len(n-i)+i) {
    z[i,j] <- z[j,i] <- e2(i, j, n, lpf, ldf, lower, upper)$value
diag(z) <- sapply(1:n, function(i) e1(i, n, lpf, ldf, -8, 8, k=2)$value)

#---Obtain the variance-covariance matrix.
v <- z - outer(y,y)
# Compute the variance-covariance matrix for the quintile means.
quantiles <- round(n * 0:5/5)
Sigma <- matrix(NA_real_, length(quantiles)-1, length(quantiles)-1)
for (a in seq_len(dim(Sigma)[1])) {
  i <- quantiles[a]+1
  j <- quantiles[a+1]
  x <- rep(1/(j-i+1), j-i+1)
  for (b in seq_len(dim(Sigma)[2])) {
    i0 <- quantiles[b]+1
    j0 <- quantiles[b+1]
    x0 <- rep(1/(j0-i0+1), j0-i0+1)
    Sigma[a,b] <- Sigma[b,a] <- crossprod(x, v[i:j, i0:j0] %*% x0)
# Verify with a simulation.
# Each column of `x` is an ordered sample.
n.sim <- 1e4
x <- apply(matrix(rnorm(n*n.sim), n, n.sim), 2, sort)

#---Compute the quantile means of each sample.
xbar <- sapply(seq_len(dim(Sigma)[1]), function(a) {
  i <- quantiles[a]+1
  j <- quantiles[a+1]
  colMeans(x[i:j, ])

#---Estimate the variance-covariance matrix from the simulated quantile means.
Sigma.hat <- cov(xbar)

#---Compare the estimate to the numerical result.  We're looking for tiny
#   values, on the order of 1/sqrt(n.sim).
hist(Sigma.hat / Sigma - 1, main="Relative errors",
     xlab=expression(hat(Sigma) / Sigma - 1))
  • 1
    $\begingroup$ Wow. What a spectacularly great answer. This not only answered my stated question, but also the next two questions the first answer raised in my mind, and a third question that I hadn't the wit to even ask. Thanks & double thanks. $\endgroup$
    – andrewH
    Nov 3, 2019 at 3:50

Data known to be normal. As an easy example, suppose you have data from a population known to be normal, but for which the population mean $\mu$ and population standard deviation $\sigma$ are unknown. (Using R.)

mu = 59; sg = 3.5
x = rnorm(10^6, mu, sg)

Pretending not to know the parameter values, you can estimate them as the sample mean $A$ and the sample SD $S.$

a = mean(x);  s = sd(x);  a;  s
[1] 59.00138
[1] 3.502753

Then use the estimated parameters to estimate the 60th and 80th percentiles (3rd and 4th quintiles). [The exact 3rd and 4th quintiles of $\mathsf{Norm}(69, 3.5)$ are 59.88671 and 61.94567, but we have "forgotten" the actual values of $\mu$ and $\sigma.]$ The statistician's 4th quintile is 59.889.

qq.mle = qnorm(c(.6,.8), a, s);  qq.mle
[1] 59.88879 61.94937

If you didn't know this is a normal distribution, somewhat less accurate estimates of these quintles could be found directly from the data.

qq.dir = quantile(x, c(.6,.8));  qq.dir
     60%      80% 
59.89191 61.95077 

Then using the quintile estimates 59.889 and 61.949 based on $A$ and $S,$ we find the average value of the data values on the interval $(59.889, 61.949).$ The planner's mean of third quintile is 60.86,

mean(x[x > 59.889 & x < 61.949])
[1] 60.86442

If I have understood your Question correctly, the next step is to try to find the 3rd quintile of a population distribution of unknown type and then to find the mean of the observations between the 3rd and 4th quintiles.

Data for unknown population type. For a very large sample from a distribution of unknown type, it may be best to get the quintiles directly from the data and then to average the relevant data values. For smaller samples, results of that will be reflect the population less accurately. Then you may want to use bootstrapping to get 95% nonparametric CIs for the quintiles and the mean of the values between.

For example, y is a sample of size $n=1000$ from a skewed population with 3rd quintile 52.366, observed 3rd quintile 51.282, and other descriptive statistics given below.

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  11.40   34.85   46.48   50.16   62.95  163.91 
q3.obs = quantile(y, .6); q3.obs

The following program makes one style of 95% nonparametric bootstrap CI $(49.33,\, 53.05)$ for the third quintile.

d.re = replicate(2000, quantile(sample(y,1000,rep=T),.6) -q3.obs)
UL = quantile(d.re, c(.975,.025))
q3.obs - UL
   97.5%     2.5% 
49.33390 53.04923 

Note: Data vector y simulated in R as:

set.seed(2019);  y = rgamma(1000, 5, .1)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.