R Script to Illustrate Direct Calculation of Bootstrap CIs in R ======================================================= 2/5/2021 > library(MASS) > gal.dat = galaxies gal.dat[78] = 26960 ## fixing a typo in the dataset: see help(galaxies) gal.dat=gal.dat/1000 ## n = 82 hist(gal.dat, prob=T, nclass=15) [ ## ordinary t interval > gal.bar = mean(gal.dat) gal.t = gal.bar + c(-1,1)*qt(0.975,81)*sd(gal.dat)/sqrt(82) [1] 19.82773 21.83519 ### Creating bootstrap intervals for the mean, WITH b=1000: > set.seed(4407) gal.boot = array(sample(gal.dat,1000*82, replace=T), c(1000,82)) ## FIRST THE BASIC PERCENTILE INTERVAL gal.bootmn = apply(gal.boot,1,mean) wquantB = quantile(gal.bootmn, prob=c(.025, .975)) gal.basic = 2*gal.bar - wquantB[2:1] 97.5% 2.5% 19.76921 21.81592 ## NEXT THE STUDENTIZED PERCENTILE INTERVAL wquantS = quantile((gal.bootmn - gal.bar)/apply(gal.boot,1,sd), prob=c(.025,.975)) gal.studpct = gal.bar - sd(gal.dat)*wquantS[2:1] 97.5% 2.5% 19.69052 21.80775 ## FINALLY THE "OTHER PERCENTILE" METHOD INTERVAL gal.pctOther = wquantB 2.5% 97.5% 19.84701 21.89372 Thus our example leads to the CI comparison > round( rbind(normal = gal.t, BasicBoot = as.vector(gal.basic), StudBoot = as.vector(gal.studpct), Pctile = as.vector(gal.pctOther)), 3) [,1] [,2] normal 19.828 21.835 BasicBoot 19.769 21.816 StudBoot 19.691 21.808 Pctile 19.847 21.894 ##--------------------------------------------------------- # THE INTERVALS CALCULATED HERE FROM FIRST PRINCIPLES CAN ALSO BE DONE # USING THE R PACKAGE boot: THE FUNCTION boot WOULD BE USED TO GENERATE # THE BOOTSTRAPPED DATA IN AN OBJECT OF CLASS "boot", TO SERVE AS # INPUT TO THE FUNCTION boot.ci ## boot.ci and Wikipedia refer to what I called Basic Percentile # method as "Basic Bootstrap Method" ## and to what I (and P. Hall) called "Other Percentile" or # "pivotal" method as "bootstrap percentile" method