function: bootstrap(x,nboot,theta,…,func=NULL) x: vector containing data. nboot:number of bootstrap samples desired. theta:function to be bootstrapped. func: optional argument specifying the functional distribution.
sub.function
sub.sum.quantile=function(x){sum(quantile(x,c(.25,.75)))/2}
sub.median.var=function(x){((quantile(x,.75,names=F)-
quantile(x,.25,name=F))/1.349)^2}
sub.Q1Q3.var=function(x){(sum(abs(x-median(x)))/(length(x)*.6745))^2}
bootstrap.biase.function=function(x,B){
sample.mean=mean(x)
sample.median=median(x)
sample.average.quantile=sub.sum.quantile(x)
sample.variance=var(x)
sample.IQR=sub.median.var(x)
sample.average.quantile=sub.Q1Q3.var(x)
sample.mu=c(sample.mean,sample.median,sample.average.quantile,
sample.variance,sample.IQR,sample.average.quantile)
resample.mean=bootstrap(x,nboot=B,mean)$thetastar
resample.median=bootstrap(x,nboot=B,median)$thetastar
resample.average.quantile=bootstrap(x,nboot=B,sub.sum.quantile)$thetastar
resample.variance=bootstrap(x,nboot=B,var)$thetastar
resample.IQR=bootstrap(x,nboot=B,sub.median.var)$thetastar
resample.MAD=bootstrap(x,nboot=B,sub.Q1Q3.var)$thetastar
temp=cbind(resample.mean,resample.median,resample.average.quantile,
resample.variance,resample.IQR,resample.MAD)
resample.mu=apply(temp,2,mean)
resample.se=apply(temp,2,sd)
bias=resample.mu-sample.mu
MSE=bias^2+resample.se^2
output=cbind(bias,resample.se,MSE)
print(output)
}
Load Package
library(gtools)
library(boot)
##
## Attaching package: 'boot'
##
## The following object(s) are masked from 'package:gtools':
##
## inv.logit, logit
library(bootstrap)
Sample Data
x=rnorm(20)
B=10000
bootstrap.biase.function(x,B)
## bias resample.se MSE
## resample.mean -0.002881 0.1906 0.03632
## resample.median -0.048464 0.2185 0.05009
## resample.average.quantile -1.098944 0.2219 1.25690
## resample.variance -0.039891 0.2388 0.05861
## resample.IQR 0.074019 0.4256 0.18660
## resample.MAD -0.054082 0.3707 0.14034
Q1Q3.sub.function=function(x){sum(quantile(x,c(.25,.75)))/2}
simple.percentile.method=function(x,B,theta,alpha){
y=bootstrap(x,B,theta)$thetastar
y=sort(y)
v1=quantile(y,alpha/2,names=F)
v2=quantile(y,(1-alpha/2),names=F)
simple.percentile.CI=c(2*theta(x)-v2,2*theta(x)-v1)
names(simple.percentile.CI)=c("Lower CL","Upper CL")
return(simple.percentile.CI)
}
main.simple.percentile.method=function(x,B,alpha,choice){
if(choice==1){
theta=mean
print("MEAN")
}
else if(choice==2){
theta=median
print("MEDIAN")
}
else if(choice==3){
theta=Q1Q3.sub.function
print("(Q1+Q3)/2")
}
CI=simple.percentile.method(x,B,theta,alpha)
print(CI)
}
Choice of Statistic Choice 1=Mean, Choice 2=Median, Choice 3=(Q1+Q3)/2
Enter Parameters
x=rnorm(20)
B=10000
alpha=.05
choice=1
main.simple.percentile.method(x,B,alpha,choice)
## [1] "MEAN"
## Lower CL Upper CL
## -0.04516 0.64305
choice=2
main.simple.percentile.method(x,B,alpha,choice)
## [1] "MEDIAN"
## Lower CL Upper CL
## -0.2301 0.5708
choice=3
main.simple.percentile.method(x,B,alpha,choice)
## [1] "(Q1+Q3)/2"
## Lower CL Upper CL
## -0.2540 0.4567
Used when the bias is small
Q1Q3.sub.function=function(x){sum(quantile(x,c(.25,.75)))/2}
sub.bootstrap.percentile.method=function(x,B,alpha,theta){
thetastar=bootstrap(x,B,theta)$thetastar
probs=(1+c(-1,1)*(1-alpha))/2
CI=quantile(thetastar,probs=probs)
return(CI)
}
main.bootstrap.percentile.method=function(x,B,alpha,choice){
if(choice==1){
print("MEAN")
theta=mean
sub.bootstrap.percentile.method(x,B,alpha,theta)
}
else if(choice==2){
print("MEDIAN")
theta=median
sub.bootstrap.percentile.method(x,B,alpha,theta)
}
else if(choice==3){
print("(Q1+Q3)/2")
theta=Q1Q3.sub.function
sub.bootstrap.percentile.method(x,B,alpha,theta)
}
}
Set Parameters
x=rnorm(20)
B=10000
alpha=.05
main.bootstrap.percentile.method(x,B,alpha,1)
## [1] "MEAN"
## 2.5% 97.5%
## -0.4340 0.3186
main.bootstrap.percentile.method(x,B,alpha,2)
## [1] "MEDIAN"
## 2.5% 97.5%
## -0.6034 0.3871
main.bootstrap.percentile.method(x,B,alpha,3)
## [1] "(Q1+Q3)/2"
## 2.5% 97.5%
## -0.4885 0.4683
Special Case of BCa
Q1Q3.sub.function=function(x){sum(quantile(x,c(.25,.75)))/2}
BcInt=function(x,B,theta,alpha){
thetahat=theta(x)
thetastar=bootstrap(x,B,theta)$thetastar
z0=qnorm(sum(thetastar<thetahat)/B)
zalpha=qnorm(alpha)
tt=pnorm(z0+z0+zalpha)
ooo=trunc(tt*B)
confpoints=sort(thetastar)[ooo]
confpoints=cbind(alpha,confpoints)
dimnames(confpoints)[[2]]=c("alpha","BcInt point")
return(confpoints)
}
main.percentile.method=function(x,B,alpha,choice){
if(choice==1){
print("mean")
theta=mean
}
else if(choice==2){
theta=median
print("median")
}
else if(choice==3){
print("(Q1+Q3)/2")
theta=Q1Q3.sub.function
}
BcInt(x,B,theta,alpha)
}
Set Parameters
x=rnorm(20)
alpha=c(.05,.95)
B=10000
main.percentile.method(x,B,alpha,1)
## [1] "mean"
## alpha BcInt point
## [1,] 0.05 -0.4631
## [2,] 0.95 0.3915
main.percentile.method(x,B,alpha,2)
## [1] "median"
## alpha BcInt point
## [1,] 0.05 -0.3537
## [2,] 0.95 0.6327
main.percentile.method(x,B,alpha,3)
## [1] "(Q1+Q3)/2"
## alpha BcInt point
## [1,] 0.05 -0.2831
## [2,] 0.95 0.5379
Adjusts percentiles to correct bias and skewness
alpha=.05
x=rnorm(20)
theta1=function(x){mean(x)}
results1=bcanon(x,100,theta1)
theta2=function(x){median(x)}
results2=bcanon(x,100,theta2)
theta3=function(x){sum(quantile(x,c(.25,.75)))/2}
results3=bcanon(x,100,theta3)
a1=results1$confpoints[results1$confpoints[,1]==alpha][2]
b1=results1$confpoints[results1$confpoints[,1]==(1-alpha)][2]
bca1=c(a1,b1)
a2=results2$confpoints[results2$confpoints[,1]==alpha][2]
b2=results2$confpoints[results2$confpoints[,1]==(1-alpha)][2]
bca2=c(a2,b2)
a3=results3$confpoints[results3$confpoints[,1]==alpha][2]
b3=results3$confpoints[results3$confpoints[,1]==(1-alpha)][2]
bca3=c(a3,b3)
BCA.result=rbind(bca1,bca2,bca3)
row.names(BCA.result)=c("bca.mean","bca.median","bca.average.quantile")
colnames(BCA.result)=c("Lower Limit","Upper Limit")
BCA.result
## Lower Limit Upper Limit
## bca.mean -0.4676 0.3015
## bca.median -0.6242 0.5514
## bca.average.quantile -0.6848 0.5191