# Approximate variance for 99.5th percentile for normal distribution

My question is similar to this one: Approximate order statistics for normal random variables

I am looking to find a formula for the variability of an arbitrary percentile of a normal distribution. The question cited does not quite solve it: - that is concerned with min and max only - I don't understand what n is

In particular, I am looking for the variability of the 99.5th percentile of a sample from a normal distribution.

Please, could you point me to a formula (with a reference) or a referred paper.

 Edits per Glen:


a) How I am collecting the samples seems irrelevant: the stdDev sample mean from normal dist is $\sigma / \sqrt{n}$ ... irrespective of anything. Should the variance of a percentile be independent of method of collection.

b) A direct derivation might be acceptable - given how quickly/frequently URLs change, I don't know whether it is sensible to cite stats.stackexchange as a source.

• Can you explain (in your question, via an edit) more about (i) how you're obtaining the estimate of your 99.5 percentile (is it via an order statistic/some average of order statistics, or in some other manner?) and (ii) roughly what sort of sample sizes you're looking at? (since that may affect the suitability of some asymptotic calculation or approximation for example; the properties at $n\approx 200$ may be somewhat different from the situation at very large $n$) – Glen_b Sep 16 '17 at 22:47
• Could you also clarify why you need a reference or a refereed paper? (Why would that be more convincing than a direct derivation, for example? For typical situations these things often yield formulas that are very simple to derive; in some cases no reliable stats journal would even publish very simply derived formulas - they'd be student exercises. Some simple formulas are found in old papers but quite a few are not) – Glen_b Sep 16 '17 at 22:50
• If you want to cite Cross Validated, you should use the cite button below the post, that will give a citation including the exact version visited. – kjetil b halvorsen Jan 21 '18 at 13:40
• It sounds like you are looking for the variation of the 99.5th percentile of a standard normal distribution as a function of sqrt(n). I would use a computer to get the fit, not algebra, but I'm not a stats major. Is this what you are asking? A variation on the extreme value distribution might be informative. Sample size is a required input. – EngrStudent Jan 21 '18 at 13:54
• Because for any fixed sample size a percentile is computed as an order statistic (or as a linear combination of adjacent order statistics), all our threads on assessing the variability of order statistics are relevant: see stats.stackexchange.com/search?q=order+statistic+variance. – whuber Jun 24 '18 at 14:15

I will try to answer without formulas many formulas because when we talk about percentile (you can google ORDER STATISTICS) the pdfs become pretty messy. I just to give you the main concepts.

1. We are talking about estimators of a quantile, in your case the 99.5th percentile.

2. Estimators are Random Variables and hence have moments.

3. You want to find the Variance of the sample 99.5th percentile from a Normal RV

The most rigorous approach I think is to evaluate an integral that is:

Suppose we call $T$ the estimator of 99.5th percentile from a Normal RV:

$\sigma^2(T) = {\displaystyle \int_{-\infty}^{+\infty} } \big(T-E(T)\big)^2f_T(t)dx$

where $E(T) = \mu(T) = {\displaystyle \int_{-\infty}^{+\infty} } Tf_T(t)dx$

As I said before $f_T(t)$ is pretty messy and you won't be able to find a close formula for the integral. Consequently you are going to evaluate the integral numerically. Just to give you an idea of the generic pdf for the $k$th ORDER STATISTIC here is what you get:

$f_{T_{(k)}}(t) =\frac{n!}{(k-1)!(n-k)!}[F_T(t)]^{k-1}[1-F_T(t)]^{n-k} f_T(t)$

So what should you do? In you question talk about approximation.

The easiest way to go about this is bootstrap. The steps are simple if we want a non sophisticated way to get some results:

1. From your ORIGINAL SAMPLE of size n calculate $\hat{\mu}$ the sample mean and ** $\hat{\sigma^2}$ the sample variance**.

2. Calculate the 99.5th percentile from the original sample.

3. Resample as many times as you want a sample of size n from a Normal distribution with mean $\hat{\mu}$ and variance $\hat{\sigma^2}$.

4. For each resample calculate the sample 99.5th percentile and store it in a vector.

5. Calculate the sample variance of this vector.

This is your approximate variance for the 99.5th percentile.

• As this page discusses, bootstrapping is unreliable at estimating extreme quantiles of distributions. The bootstrap results then tend to represent the vagaries of the sample at hand rather than the properties of the underlying distribution. The OP, however, seems more interested in theoretical results for a normal distribution rather than analysis of a particular sample. In that case, repeated random sampling from a normal distribution would provide an approach similar to what you suggest. – EdM Jun 24 '18 at 14:43

How I (non-stats-guy) would do this.

1) State the question I'm engaging

Given a standard normal distribution, track how variance in 99.5th percentile changes with sample size.

2) Describe basic approach

I'm going to simulate it, and make a graph.

If I were ambitious I would fit a generalized extreme value distribution at each n, and compute the variance of the distribution, and try to relate how parameters change with increasing n.

I feel like the true number of points informing a variance should be greater than 2. While the number of repeats selected is ~300, when $n < 100$ the percentile is informed by two points. The 99.5th percentile gets to be something other than interpolation between the last two samples in the tail after $n=100$.

3) execute

Here is my code:

#range of sample size "n"
n_min <- 100
n_max <- 2000

#number of samples at each value of "n"
N_samp <- 1000

#stage for loop
p995 <- numeric(length = N_samp)
summ <- numeric(length = n_max-n_min)

#big loop
for(i in n_min:n_max){

for(j in 1:N_samp){

#draw
y <- rnorm(n = i)

#compute 99.5th percentile
p995[j] <- quantile(x=y, probs = 0.995)

}
summ[1+i-n_min] <- var(p995)
}

x <- 1/sqrt(n_min:n_max)
y <- summ

est_1 <- lowess(y~x)
est_2 <- lm(y~x)

plot(x,y)
lines(est_1$x,est_1$y, col="Red", lwd=2)
lines(x,predict(est_2), col="Blue", lwd=2)
grid()

summary(est_2)


Here is the summary:

> summary(est_2)

Call:
lm(formula = y ~ x)

Residuals:
Min         1Q     Median         3Q        Max
-0.0127940 -0.0011302  0.0000608  0.0010720  0.0130319

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.0203914  0.0001284  -158.8   <2e-16 ***
x            1.3395199  0.0032323   414.4   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.002192 on 1899 degrees of freedom
Multiple R-squared:  0.9891,    Adjusted R-squared:  0.9891
F-statistic: 1.717e+05 on 1 and 1899 DF,  p-value: < 2.2e-16


Here is the plot:

4) look back on execution

The relationship isn't linear. It is close, but clearly not.

You have a linear approximation in the parameters of the linear fit, but it is not the "truth", whatever that is.

Update:
I tried a cubic instead of a linear, and after about n=400, there was a pretty good fit.

updated lines in code:

#range of sample size "n"
n_min <- 400
n_max <- 10000

#number of samples at each value of "n"
N_samp <- 1000


and

est_1 <- lowess(y~x, f = 0.05)  #really tight span value
est_2 <- lm(y~1 + x +I(x^2) +I(x^3))

plot(x,y, pch=19, cex=0.7, col="Green")
lines(est_1$x,est_1$y, col="Red", lwd=2)
lines(x,predict(est_2), col="Blue", lwd=2)
grid()

summary(est_2)


yields:

> summary(est_2)

Call:
lm(formula = y ~ 1 + x + I(x^2) + I(x^3))

Residuals:
Min         1Q     Median         3Q        Max
-0.0056997 -0.0001487 -0.0000011  0.0001404  0.0073250

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.979e-04  8.409e-05  -3.543 0.000398 ***
x            3.235e-02  1.195e-02   2.708 0.006788 **
I(x^2)       2.419e+01  5.042e-01  47.971  < 2e-16 ***
I(x^3)      -1.191e+02  6.390e+00 -18.641  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.0005175 on 9597 degrees of freedom
Multiple R-squared:  0.9951,    Adjusted R-squared:  0.9951
F-statistic: 6.474e+05 on 3 and 9597 DF,  p-value: < 2.2e-16


and this:

This gets up up around 10k samples. Given the decreasing width of for decreasing 1/sqrt(n) it seems to converge. This weakly suggests that the formula might extrapolate well.