What is the variance of the mean, conditional on being between two order statistics or quantiles? 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 
 A: 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.

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
set.seed(17)
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))

