# How to obtain a confidence interval for a percentile?

I have a bunch of raw data values that are dollar amounts and I want to find a confidence interval for a percentile of that data. Is there a formula for such a confidence interval?

## 3 Answers

This question, which covers a common situation, deserves a simple, non-approximate answer. Fortunately, there is one.

Suppose $X_1, \ldots, X_n$ are independent values from an unknown distribution $F$ whose $q^\text{th}$ quantile I will write $F^{-1}(q)$. This means each $X_i$ has a chance of (at least) $q$ of being less than or equal to $F^{-1}(q)$. Consequently the number of $X_i$ less than or equal to $F^{-1}(q)$ has a Binomial$(n,q)$ distribution.

Motivated by this simple consideration, Gerald Hahn and William Meeker in their handbook Statistical Intervals (Wiley 1991) write

A two-sided distribution-free conservative $100(1-\alpha)\%$ confidence interval for $F^{-1}(q)$ is obtained ... as $[X_{(l)}, X_{(u)}]$

where $X_{(1)}\le X_{(2)}\le \cdots \le X_{(n)}$ are the order statistics of the sample. They proceed to say

One can choose integers $0 \le l \le u \le n$ symmetrically (or nearly symmetrically) around $q(n+1)$ and as close together as possible subject to the requirements that $$B(u-1;n,q) - B(l-1;n,q) \ge 1-\alpha.\tag{1}$$

The expression at the left is the chance that a Binomial$(n,q)$ variable has one of the values $\{l, l+1, \ldots, u-1\}$. Evidently, this is the chance that the number of data values $X_i$ falling within the lower $100q\%$ of the distribution is neither too small (less than $l$) nor too large ($u$ or greater).

Hahn and Meeker follow with some useful remarks, which I will quote.

The preceding interval is conservative because the actual confidence level, given by the left-hand side of Equation $(1)$, is greater than the specified value $1-\alpha$. ...

It is sometimes impossible to construct a distribution-free statistical interval that has at least the desired confidence level. This problem is particularly acute when estimating percentiles in the tail of a distribution from a small sample. ... In some cases, the analyst can cope with this problem by choosing $l$ and $u$ nonsymmetrically. Another alternative may be to use a reduced confidence level.

Let's work through an example (also provided by Hahn & Meeker). They supply an ordered set of $n=100$ "measurements of a compound from a chemical process" and ask for a $100(1-\alpha)=95\%$ confidence interval for the $q=0.90$ percentile. They claim $l=85$ and $u=97$ will work.

The total probability of this interval, as shown by the blue bars in the figure, is $95.3\%$: that's as close as one can get to $95\%$, yet still be above it, by choosing two cutoffs and eliminating all chances in the left tail and the right tail that are beyond those cutoffs.

Here are the data, shown in order, leaving out $81$ of the values from the middle:

$$\matrix{ 1.49&1.66&2.05&\ldots&\mathbf {24.33}&24.72&25.46&25.67&25.77&26.64\\ 28.28&28.28&29.07&29.16&31.14&31.83&\mathbf{33.24}&37.32&53.43&58.11}$$

The $85^\text{th}$ largest is $24.33$ and the $97^\text{th}$ largest is $33.24$. The interval therefore is $[24.33, 33.24]$.

Let's re-interpret that. This procedure was supposed to have at least a $95\%$ chance of covering the $90^\text{th}$ percentile. If that percentile actually exceeds $33.24$, that means we will have observed $97$ or more out of $100$ values in our sample that are below the $90^\text{th}$ percentile. That's too many. If that percentile is less than $24.33$, that means we will have observed $84$ or fewer values in our sample that are below the $90^\text{th}$ percentile. That's too few. In either case--exactly as indicated by the red bars in the figure--it would be evidence against the $90^\text{th}$ percentile lying within this interval.

One way to find good choices of $l$ and $u$ is to search according to your needs. Here is a method that starts with a symmetric approximate interval and then searches by varying both $l$ and $u$ by up to $2$ in order to find an interval with good coverage (if possible). It is illustrated with R code. It is set up to check the coverage in the preceding example for a Normal distribution. Its output is

Simulation mean coverage was 0.9503; expected coverage is 0.9523

The agreement between simulation and expectation is excellent.

#
# Near-symmetric distribution-free confidence interval for a quantile q.
# Returns indexes into the order statistics.
#
quantile.CI <- function(n, q, alpha=0.05) {
#
# Search over a small range of upper and lower order statistics for the
# closest coverage to 1-alpha (but not less than it, if possible).
#
u <- qbinom(1-alpha/2, n, q) + (-2:2) + 1
l <- qbinom(alpha/2, n, q) + (-2:2)
u[u > n] <- Inf
l[l < 0] <- -Inf
coverage <- outer(l, u, function(a,b) pbinom(b-1,n,q) - pbinom(a-1,n,q))
if (max(coverage) < 1-alpha) i <- which(coverage==max(coverage)) else
i <- which(coverage == min(coverage[coverage >= 1-alpha]))
i <- i[1]
#
# Return the order statistics and the actual coverage.
#
u <- rep(u, each=5)[i]
l <- rep(l, 5)[i]
return(list(Interval=c(l,u), Coverage=coverage[i]))
}
#
# Example: test coverage via simulation.
#
n <- 100      # Sample size
q <- 0.90     # Percentile
#
# You only have to compute the order statistics once for any given (n,q).
#
lu <- quantile.CI(n, q)$Interval # # Generate many random samples from a known distribution and compute # CIs from those samples. # set.seed(17) n.sim <- 1e4 index <- function(x, i) ifelse(i==Inf, Inf, ifelse(i==-Inf, -Inf, x[i])) sim <- replicate(n.sim, index(sort(rnorm(n)), lu)) # # Compute the proportion of those intervals that cover the percentile. # F.q <- qnorm(q) covers <- sim[1, ] <= F.q & F.q <= sim[2, ] # # Report the result. # message("Simulation mean coverage was ", signif(mean(covers), 4), "; expected coverage is ", signif(quantile.CI(n,q)$Coverage, 4))

• $u-1$ would be in the range of plausible values for number of $X_i$ less than or equal to $F^{-1}(q)$, whereas $u$ is not. Why do we take $X_{(u)}$ as the upper limit and not $X_{(u-1)}$. Would the CI not be conservative otherwise? Mar 8, 2021 at 8:10
• @retodomax I'm afraid I don't understand your argument, perhaps because what you mean by "plausible values" is unclear.
– whuber
Mar 8, 2021 at 14:35
• As I understand, ${l, l+1, \dots, u-1}$ are all "plausible values" for a random variable following a Binomial(n,q) distribution. $u$ is not included and therefore, it is unlikely to observe $u$ observations which are less than or equal to the quantile. Still, we use $X_{u}$ as our upper bound and not $X_{u-1}$? Mar 8, 2021 at 16:12
• @whuber, I second the point by retodomax: You're contradicting yourself in the explanation and the example. Explanation: "The expression at the left is the chance that a Binomial(𝑛,𝑞) variable has one of the values {𝑙,𝑙+1,…,𝑢−1}". Example: you're given 𝑙=85 and 𝑢=97. But then the values included in your interval are {𝑙,𝑙+1,…,𝑢} instead of {𝑙,𝑙+1,…,𝑢−1}. Mar 17, 2021 at 12:46
• @jollycat It's easy to make one-off mistakes like that, but I recall being careful not to and careful with the example and the software to match the book. I cannot find the explanation to which you refer. A close look at the plot shows it colors the cases 85 through 96, inclusive, in cyan; and my explanation explicitly characterizes the upper tail as "97 or more." This doesn't appear to match what you or retodomax seem to be saying.
– whuber
Mar 17, 2021 at 13:21

Derivation

The $\tau$-quantile $q_\tau$ (this is the more general concept than percentile) of a random variable $X$ is given by $F_X^{-1}(\tau)$. The sample counterpart can be written as $\hat{q}_\tau = \hat{F}^{-1}(\tau)$ -- this is just the sample quantile. We are interested in the distribution of:

$\sqrt{n}(\hat{q}_\tau - q_\tau)$

First, we need the asymptotic distribution of the empirical cdf.

Since $\hat{F}(x) = \frac{1}{n} \sum 1\{X_i < x\}$, you can use the central limit theorem. $1\{X_i < x\}$ is a bernoulli random variable, so the mean is $P(X_i < x) = F(x)$ and the variance is $F(x)(1-F(x))$.

$\sqrt{n}(\hat{F}(x) - F(x)) \rightarrow N(0, F(x)(1-F(x))) \qquad (1)$

Now, because inverse is a continuous function, we can use the delta method.

[**The delta method says that if $\sqrt{n}(\overline{y} - \mu_y) \rightarrow N(0,\sigma^2)$, and $g(\cdot)$ is a continuous function, then $\sqrt{n}(g(\overline{y}) - g(\mu_y)) \rightarrow N(0, \sigma^2 (g'(\mu_y))^2)$ **]

In the left hand side of (1), take $x=q_\tau$, and $g(\cdot) = F^{-1}(\cdot)$

$\sqrt{n}(F^{-1}(\hat{F}(q_\tau)) - F^{-1}(F(q_\tau))) = \sqrt{n}(\hat{q}_\tau - q_\tau)$

[** note that there is a bit of a slight of hand in the last step because $F^{-1}(\hat{F}(q_\tau)) \neq \hat{F}^{-1}(\hat{F}(q_\tau)) = \hat{q}_\tau$, but they are the asymptotically equal if tedious to show **]

Now, apply the delta method mentioned above.

Since $\frac{\textrm{d}}{\textrm{d}x} F^{-1}(x) = \frac{1}{f(F^{-1}(x))}$ (inverse function theorem)

$\sqrt{n}(\hat{q}_\tau - q_\tau) \rightarrow N\left(0, \frac{F(q_\tau)(1-F(q_\tau))}{f(F^{-1}(F(q_\tau)))^2}\right) = N\left(0, \frac{F(q_\tau)(1-F(q_\tau))}{f(q_\tau)^2}\right)$

Then, to construct the confidence interval, we need to calculate the standard error by plugging in sample counterparts of each of the terms in the variance above:

Result

So $se(\hat{q}_\tau) = \sqrt{\frac{\hat{F}(\hat{q}_\tau)(1-\hat{F}(\hat{q}_\tau))}{n \hat{f}(\hat{q}_\tau)^2}} =$ $\sqrt{\frac{\tau (1 - \tau)}{n \hat{f}(\hat{q}_\tau)^2}}$

And $CI_{0.95}(\hat{q}_\tau) = \hat{q}_\tau \pm 1.96 se(\hat{q}_\tau)$

This will require you to estimate the density of $X$, but this should be pretty straightforward. Alternatively, you could bootstrap the CI pretty easily too.

• Could you expand your answer with contents from the linked article? Links may not work forever and then this answer would become less useful
– Andy
May 23, 2014 at 17:29
• What is the advantage of this asymptotic result based on density estimates compared to the distribution free c.i.based on the binomial distribution? May 23, 2014 at 21:27
• Is this still based on the article you linked originally? May 24, 2014 at 17:25
• Yes, should I add that link back in? I think this is a well known result. I've seen it in class before and it is not hard to find by google. In a case like this, is it better to link to it or type it up, or both? May 25, 2014 at 17:18
• I'd say both, and that you should edit it back in if this is quoted / derived entirely from it for the sake of proper attribution. Otherwise it might not matter whether you edit it, but in general, Stack Exchange policy is to discourage link-only answers to avoid link rot and as a matter of principle (the idea is to be an independent repository, not a link index – but I'm not sure how much of that scenario is more than an imaginary "slippery slope"). May 26, 2014 at 7:16

A brute-force computing intensive solution is to use the bootstrap resampling method. The following function returns the bootstrap confidence intervals of a quantile.

quantile.CI.via.bootstrap <- function(x, p, alpha = 0.1) {
## Purpose:
##   Calculate a two-sided confidence interval with confidence level of (1 - alpha) for
##   a quantile, based on the (computing intensive) bootstrap resampling method.
##
## Arguments:
##   - x: a vector of values, representing a data sample.
##   - p: probability cutpoint for the quantile (between 0 and 1).
##   - alpha: type I error level (default to 0.1 so a 90% CI is calculated)
##
## Return:
##   - CI: the lower and upper limits of the two-sided CI.

q <- quantile(x, probs = p)
message("Quantile Point Estimate = ", q, " (Probability Cutpoint = ", p, ")\n")

## Bootstrap resampling with 2000 replications
library(boot)
set.seed(1)
b <- boot(x, function(x, i) quantile(x[i], probs = p), R = 2000)

boot.ci(b, conf = 1 - alpha, type = c("norm", "basic", "perc", "bca"))

}
if (F) {                                # Unit Test
x <- 1:100
p <- 0.9
alpha <- 0.05
quantile.CI.via.bootstrap(x, p, alpha)
## Intervals :
## Level      Normal              Basic
## 95%   (84.50, 96.34 )   (85.10, 97.10 )

## Level     Percentile            BCa
## 95%   (83.1, 95.1 )   (83.3, 95.1 )
}


You may note the "Basic" bootstrap method returns an interval [85, 97] that aligns with the analytical method (binomial distribution) in the previous post.

• A little analysis will relate these bootstrap intervals to the symmetrical intervals given by Hahn & Meeker. After all, the result does not depend on the actual values in the dataset (continuing to assume all are distinct), but only on how many there are and on the quantile. Ergo, it's unnecessary to carry out the bootstrap calculations. Just use H&M's approach directly.
– whuber
15 hours ago