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.
