Confidence Interval of p-Quantile from Empirical CDF

Question

I am trying to provide an interval estimate for the 0.8-quantile of some numeric data, which is assumed to be an IID sample from some unknown, continuous distribution. I constructed an Empirical CDF for my data (black line in the image below). I used the information contained in this Wikipedia article and this review to construct point-wise 95% confidence intervals for the true CDF at many values of x. (These intervals show up as the red lines in the image below. Note that I did NOT find simultaneous confidence bands.)

How do you turn bounds on the CDF into bounds on quantiles? I found where the bounds on the CDF reach 0.8. Can I take the x-values of these intersections as bounds for a 95% confidence interval on the 0.8-quantile? In other words, is [28 000, 36 000] a valid 95% CI for the 0.8-quantile? Is it valid to invert the bounds on the CDF in this manner in order to obtain a confidence interval for a quantile?

EDIT: @BruceET The green lines show the method @BruceET used in his answer.

Using @BruceET's method, the CI for the 0.8-quantile would be [28 500, 36 500]. I think that my approach in the blue lines finds one-sided confidence intervals at different values of x and then pieces those together to find a two-sided CI for the 0.8-quantile. @BruceET's method seems more straightforward. Would one be preferred to the other? More importantly for the purposes of this post, what makes it valid to transform a CI for the CDF into a CI for the p-quantile of the population?

Answer 1

Data. Because I do not have your dataset at hand, I will use a sample of size $100$ from $\mathsf{Norm}(\mu = 50, \sigma = 5),$ which has 80th percentile 54.21, as an example. Sampling in R:

set.seed(2020)
x = rnorm(100, 50, 5)
quantile(x, .8)
     80% 
54.62839 
qnorm(.8, 50, 5)
[1] 54.20811

The 80th quantile $54.62$ of the sample of size 100 is close to the 80th percentile of the population. If we did not know the population percentile, the question would be how to find a 95% CI for the population percentile based on the sample.

Binomial method. Roughly speaking, I believe the method shown in the figure of your question is to find a 95% CI for the true proportion $p$ of population values below 54.21. A Wald CI gives $(0.72, 0.88).$ [Various styles of binomial CIs are possible.]

pm = c(-1,1)
ci.p = .8 + pm*1.96*sqrt(.8*.2/100);  ci.p
[1] 0.7216 0.8784

The find the sample values with quantiles 0.72 and 0.88:

quantile(x, ci.p)
  72.16%   87.84% 
53.21117 56.21680 

So that a 95% CI for the 80th percentile of the population is $(53.21, 56.22).$ [Because we know the true 80th population percentile to be 54.21, we know that this CI covers the true value.]

Bootstrap method. Another method is to use a bootstrap. By taking $B = 5000$ 're-samples' from the 100 observations, we can find their 80th percentiles ad use them to make a CI for the population value. [A 're-sample' is a sample of size $n=100$ taken with replacement from the sample x.]

q.8 = replicate(5000, quantile(sample(x,100,rep=T),.8))
mean(q.8)
[1] 54.56447
quantile(q.8, c(.025,.975))
2.5%    97.5% 
53.03377 56.03073 

The average of these bootstrapped 80th re-sample percentiles is 54.56 and the middle 95% of them span an 95% bootstrap CI $(53.0, 56.0).$ Various styles of bootstraps are possible. The CI obtained here is not much different from the CI $(53.21, 56.22),$ obtained using the binomial method above.

Addendum: With the same 100 observations x as above, I might have guessed that the data are nearly normal. Looking at a normal probability plot (normal Q-Q plot), not show, points fall nearly along a straight line. A Shapiro-Wilk test finds the data are consistent with normal. So a parametric bootstrat might use re-samples repeatedly to estimate $\mu$ and $\sigma$ and then find the 80th percentile of a normal distribution with those estimated parameters. The R code below assumes data x are already available in R. Thus a 95% parametric bootstrap CI for the 90th percentile of the population is $(53.95, 56.50),$ not much different from CIs found by other methods.

set.seed(1234)
B = 5000;  q.80 = numeric(B)
for (i in 1:B) {
  x.re = sample(x, 100, rep=T)
  q.80[i] = qnorm(.8, mean(x.re), sd(x.re)) }
quantile(q.80, c(.025, .975))
    2.5%    97.5% 
53.95346 56.50207