# Bootstrap confidence interval in R using 'replicate' and 'quantile'

I have around one thousand measurements (numbers). All these measurements are my observations and I have calculated a 95% confidence interval for the mean and for the variance by using the normal formulas (without software).

Then, I have used the replicate function in R with around one hundred thousand simulations for my measurements (observations) and the parameter "replace" has been set to "true".

After that, I used the apply function with the three parameters: my measurements as the data, 2 as the margin and mean as the function.

Then to get a 95% confidence interval this way, I used the quantile function containing the variable for the apply function I used, and then 0.025 and 0.975 combined as the second parameter for the quantile function.

In that way, I got almost exactly the same 95% confidens interval as calculated with the normal formula (without software).

So now I wanted to do exactly the same thing for the variance, i.e. use replicate, apply and quantile to get a 95% confidence interval for the variance. So I just changed the third parameter "mean" to "var" in the apply function. I then noticed that the outputted 95% confidence interval for the variance (from the quantile function) is a little bit different (it is wider) than the one I calculated by the normal formula without software.

So my question is:

Did I use the replicate, apply and quantile function correctly for the variance confidence interval? I know I did use the functions right for the mean, since I got almost exactly the same result there as calculated by normal formula.

If I did use the functions correctly for the variance, why is the confidence interval a little bit different? Is it because as sample size increases, there might be more observations far away from the mean resulting in a bigger variance?

• Hello and welcome to math.stackexchange. It looks like you are interested in comparing bootstrap and formula confidence intervals for the mean and the variance. From your description, it looks to me that you are doing the right thing. As to the difference between the bootstrap and formula confidence intervals for the variance, this is most likely due to the fact that your distribution is not normal. The formula confidence interval for the variance is not as robust to such deviations as the formula confidence for the mean. That is, the bootstrap CI is more accurate. – Hans Engler Aug 6 '17 at 16:15
• Hi. Thank you for your comment. But the distribution of my data (my observations) is definitely normal. I have made both a histogram of it (and it forms a bell curve) and a Q-Q Normal plot containing an almost straight line through the dots. – John A Aug 6 '17 at 16:33
• Try the same thing with synthetic data from $N(0,1)$ and see if you get the same effects. – Hans Engler Aug 6 '17 at 18:20

From your description, I am not sure exactly what you are doing. But maybe I can give you some helpful ideas.

1) For normal data, the standard CI for the population mean $\mu$ uses Student's t distribution which is symmetrical. The 95% CI is of the form $\bar X \pm t^*S/\sqrt{n},$ where $t^*$ cuts 2.5% from the upper tail of $\mathsf{T}(\nu= n-1)$.

By contrast, the standard CI for the population variance $\sigma^2$ uses the chi-squared distribution which is not symmetrical. The 95% CI is of the form $\left(\frac{(n-1)S^2}{U}, \frac{(n-2)S^2}{L}\right),$ where $L$ and $U$ cut 2.5% of the probability from the lower and upper tails, respectively, of $\mathsf{Chisq}(\nu = n-1).$

Unless your bootstrapping procedure corrects for the bias in the case of the variance (with its skewed distribition), you cannot expect an accurate result.

2) It is not exactly 'fair' to compare CIs (using Student's t and chi-squared distributions) based on the assumption data are normal with results of nonparametric bootstrap CIs. The assumption that data are normal provides 'information' for the t and chi-squared intervals that is not used in nonparametric bootstrapping.

3) In making bootstrap CIs for variation, I have found it better to find CIs for $\sigma$ rather than $\sigma^2,$ possibly because the former has the same units as the data. Also, for scale parameters, I have found that it often works better to bootstrap ratios rather than differences.

4) Finally, along lines of @HansEngler's suggestion, I will generate $n = 1000$ observations from $\mathsf{Norm}(\mu = 100, \sigma=15),$ find the chi-sqared CI for $\sigma$ and compare it with a bias-corrected nonparametric bootstrap.

First, here are my fake data and the corresponding standard CI for $\sigma,$ which turns out to be $(9.55, 10.43).$ [I have provided set.seed statements, so that you can replicate the exact simulations I have used.]

set.seed(1234)
n = 1000;  mu = 100;  sg = 10;  x = rnorm(n, mu, sg)
a = mean(x); s = sd(x);  a;  s
## 99.73403    # sample mean
## 9.973377    # sample SD
# CI for sg
UL = qchisq(c(.975,.025), 999); UL
## 1088.487  913.301
sqrt((n-1)*var(x)/UL)
##  9.554619 10.430810  # CI for sg (includes 9.973, as it must)


Now for the (bias corrected) nonparametric bootstrap CI: I will bootstrap the ratio $R = S/\sigma.$ If I knew the distribution of $R,$ then I could find $L$ and $U$ with $P(L \le R = S/\sigma \le U) = 0.95$ so that $P(S/U \le \sigma \le S/L) = .95$ and a 95% CI for $\sigma$ would be of the form $(S/U,\, S/L).$ By bootstrapping $R,\,$ I can estimate $U$ and $L.$ I use the observed $s = 9.97$ temporarily as a proxy for unknown $\sigma.$ Suffixes .re indicate bootstrapped quantities. [On this site with relatively few experienced R users, I try to use only the most fundamental R functions; I will leave it to you to write more elegant R code, which is obviously possible.]

set.seed(1235)
B = 10^5;  r = numeric(B)
for(i in 1:B) {
s.re = sd(sample(x,n,repl=T))
r[i] = s.re/s }
L.re = quantile(r, .025);  U.re = quantile(r, .975)
c(s/U.re, s/L.re)
##    97.5%      2.5%
## 9.537866 10.462202


So the nonparametric bootstrap CI is $(9.54, 10.46),$ which is not a bad match for the standard chi-squared CI $(9.55, 10.43)$ obtained above. [The large sample size ($n$ = 1000) has mostly obviated my comment in (2) above, which still stands for smaller $n$. Even so, it is possible that the information that the data are normal accounts for the fact that the normal-based CI is a bit shorter.]

Note: If you do a parametric bootsrap CI, you may get closer to the result for the chi-squared CI. In parametric bootstrapping, one re-samples from a parametric distribution (here normal) with parameters suggested by the data (rather than from the data themselves). In particular, in the code above one would substitute the code s.re = sd(rnorm(n, a, s)) for the line with sample. With this change and bootstrap seed 1066, I got the 95% parametric bootstrap CI $(9.55, 10.43).$