一、引言
二、算法手段:两类基于bootstrap思想的自编函数(有关均值分析)
三、基于第一类函数的覆盖率性质分析(有关均值分析)
四、基于第二类函数的覆盖率性质分析(有关均值分析)
五、基础函数(有关偏度分析)
六、性质验证:验证各个参数对覆盖概率的影响(有关偏度分析)
七、结语
人员分工:尹宏博PB17010389(均值分析:二、三、四),孙文明PB17010388(偏度分析:五、六)
Bootstrap方法作为非参数估计中的重要方法,在统计学中有着广泛的应用,bootstrap区间估计便是其中之一。在本研究中我们将采用bootstrap区间估计的思想方法进行有关CI覆盖率研究的自编算法实现,我们关心的统计量主要是均值和偏度,主要研究的影响因素为:样本(包括样本量与相关参数)、bootstrap次数、bootstrap区间类型。希望可以使用此算法对某些分布的覆盖率问题进行探索,希望以此揭示一些特定分布内蕴性质与意义。
针对我们感兴趣的置信区间覆盖率问题,自然地定义是使用bootstrap方法做出多组置信区间,以置信区间是否覆盖我们关心的标准统计量数值的概率作为区间覆盖率。这是目前广泛认可和使用的区间覆盖率的算法。另外我们还会很自然地想到一个对应的命题:我们只使用bootstrap做出一组置信区间,并对样本重复抽样计算出多组统计量,在区间内的统计量的比例这种算法从直观想象上和上一种相差不算太大。
我们或许可以猜测:这两种算法计算出来的覆盖率都会与置信水平较为接近。可是,针对不同的分布类型、不同的置信区间类型、不同的样本数据、不同的统计量,这其中的变量非常之多且繁杂,究竟这些变量在bootstrap方法下扮演了何种角色,又会对区间覆盖率产生怎样的影响?本着以上的目标,我们自编了两类函数,它们分别基于上面阐述的两种算法思想,其功能与实现结果侧重点也不尽相同(完全自主编写,不依赖BOOT包和BOOT函数)。
功能:(对样本性质进行分析的函数)对输入样本进行多次重抽样,对三种区间分别构造一个CI(normal,basic,percentile),并重复重抽样过程构造一组统计量,计算落入CI的概率,并进行三组绘图。
输入参数:x=样本;statistic=统计量,B=bootstrap次数,alpha=置信水平
返回:三组区间,三组区间的左右遗失率、总遗失率及其分别的三组覆盖图像
目标变量:统计量、bootstrap次数、样本信息
#自主编写的计算置信区间的函数BOOT
#BOOT的参数如下:x=样本;statistic=统计量,B=bootstrap次数,alpha=显著性水平
BOOT<-function(x=x0,statistic=mean,B=20,alpha=0.05){
statistic_hat<-statistic(x)
size0<-length(x)
#求statistic的一组有放回抽样值(我们采用的是bootstrap方法,这里我们bootstrap默认次数设为300)
statistic_hat_boot<-c(1:B)#对statistic_hat做bootstrap得到的一组statistic值
for(i in 1:B){
indices<-sample(1:size0,size0,replace = TRUE)#有放回抽样
x_bootstrap<-x[indices]
statistic_hat_boot[i]<-statistic(x_bootstrap)
}
# the standard normal bootstrap confidence interval
sd_hat_forstatistic<-sd(statistic_hat_boot)
standard_norm_CI<-c(statistic_hat+qnorm(alpha/2)*sd_hat_forstatistic,statistic_hat-qnorm(alpha/2)*sd_hat_forstatistic)
# the basic bootstrap confidence interval
basic_boot_CI<-c(2*statistic_hat-quantile(statistic_hat_boot,1-alpha/2),2*statistic_hat-quantile(statistic_hat_boot,alpha/2))
# the percentile confidence interval
percentile_CI<-c(quantile(statistic_hat_boot,alpha/2),quantile(statistic_hat_boot,1-alpha/2))
#制作返回列表
CI<-rbind(standard_norm_CI,basic_boot_CI,percentile_CI)
aa<-1-alpha/2
bb<-alpha/2
colnames(CI)<-c(bb,aa)
#计算the proportion of times that the confidence intervals miss on the left and right
#这里我们做m次重抽样(m=10000)并计算对应的statistic,以这个数组计算miss on the left and right的概率
CI<-data.frame(CI)
CI[,"left missing"]<-rep(0,nrow(CI))
CI[,"right missing"]<-rep(0,nrow(CI))
CI[,"total missing"]<-rep(0,nrow(CI))
m<-1000
sample_statistic0<-c(1:m)
for(i in 1:m){
indices<-sample(1:size0,size0,replace = TRUE)#有放回抽样
x_bootstrap<-x[indices]
statistic_sample<-statistic(x_bootstrap)
sample_statistic0[i]<-statistic_sample
for(j in 1:nrow(CI)){
if(statistic_sample<CI[j,1]){
CI[j,"left missing"]<-CI[j,"left missing"]+1
}else if(statistic_sample>CI[j,2]){
CI[j,"right missing"]<-CI[j,"right missing"]+1
}
}
}
CI$`left missing`<-CI$`left missing`/m
CI$`right missing`<-CI$`right missing`/m
CI$`total missing`<-CI$`left missing`+CI$`right missing`
par(mfrow = c(3,1))
plot_sample0<-plot(x=sample_statistic0,y=rep(0,m),xlab="bootstrap计算出的参数",ylab=" ",main = "bootstrap方法计算的参数分布(standard normal bootstrap confidence interval)")
plot_sample0<-rect(CI[1,1],-0.5,CI[1,2],0.5,col="gray")
plot_sample0<-points(x=sample_statistic0,y=rep(0,m))
plot_sample1<-plot(x=sample_statistic0,y=rep(0,m),xlab="bootstrap计算出的参数",ylab=" ",main = "bootstrap方法计算的参数分布(basic bootstrap confidence interval)")
plot_sample1<-rect(CI[2,1],-0.5,CI[2,2],0.5,col="gray")
plot_sample1<-points(x=sample_statistic0,y=rep(0,m))
plot_sample2<-plot(x=sample_statistic0,y=rep(0,m),xlab="bootstrap计算出的参数",ylab=" ",main = "bootstrap方法计算的参数分布(percentile confidence interval)")
plot_sample2<-rect(CI[3,1],-0.5,CI[3,2],0.5,col="gray")
plot_sample2<-points(x=sample_statistic0,y=rep(0,m))
return(CI)
}
例:第一类函数返回值示意:
x0<-rnorm(100,1,2)
BOOT(x=x0)
## X0.025 X0.975 left missing right missing
## standard_norm_CI 0.4632975 1.406034 0.008 0.009
## basic_boot_CI 0.4193788 1.273425 0.005 0.042
## percentile_CI 0.5959062 1.449953 0.046 0.005
## total missing
## standard_norm_CI 0.017
## basic_boot_CI 0.047
## percentile_CI 0.051
功能:(对具体分布类型进行分析的函数)对输入的对应分布的参数值构造多组样本,对三种区间分别构造一组CI(normal,basic,percentile),计算标准统计值落入这组CI的概率,并进行绘图。
函数分类:这里构造的分析分布性质的函数是不能混用的,本阶段着重研究二项与正态分布,故只做了分析这两种分布的两种函数
输入参数(以正态为例):size0=样本量,mean0=均值,sd0=标准差,statistic=统计量,standard_statistic=标准统计量值,B=bootstrap次数,alpha=显著性水平
返回:三组区间的覆盖率,及覆盖图像
正态分布的第二类函数:
#bootstrap求覆盖率函数(针对正态分布)
#BOOT_COVER_PROB_NORM的参数如下:size0=样本量,mean0=样本均值,sd0=样本标准差,statistic=统计量,standard_statistic=总体的标准参数值,B=bootstrap次数,alpha=显著性水平
BOOT_COVER_PROB_NORM<-function(size0=100,mean0=0,sd0=1,statistic=mean,standard_statistic=0,B=500,alpha=0.05){
sn_cp0<-c();bb_cp0<-c();per_cp0<-c()
sn_cp<-0;bb_cp<-0;per_cp<-0
sn_l<-c();sn_r<-c();bb_l<-c();bb_r<-c();per_l<-c();per_r<-c()
for(l in 1:500){
#求statistic的一组有放回抽样值(我们采用的是bootstrap方法,这里我们bootstrap默认次数设为500)
x<-rnorm(size0,mean0,sd0)
statistic_hat<-statistic(x)
statistic_hat_boot<-c(1:B)#对statistic_hat做bootstrap得到的一组statistic值
for(i in 1:B){
indices<-sample(1:size0,size0,replace = TRUE)#有放回抽样
x_bootstrap<-x[indices]
statistic_hat_boot[i]<-statistic(x_bootstrap)
}
# the standard normal bootstrap confidence interval
sd_hat_forstatistic<-sd(statistic_hat_boot)
sn_l[l]<-statistic_hat+qnorm(alpha/2)*sd_hat_forstatistic
sn_r[l]<-statistic_hat-qnorm(alpha/2)*sd_hat_forstatistic
# the basic bootstrap confidence interval
bb_l[l]<-2*statistic_hat-quantile(statistic_hat_boot,1-alpha/2)
bb_r[l]<-2*statistic_hat-quantile(statistic_hat_boot,alpha/2)
# the percentile confidence interval
per_l[l]<-quantile(statistic_hat_boot,alpha/2)
per_r[l]<-quantile(statistic_hat_boot,1-alpha/2)
#存储是否覆盖数组(由于总体的statistic数据未知,这里我们用原样本的statistic数据代替之)
sn_cp0[l]<-(sn_l[l]<standard_statistic)&(sn_r[l]>standard_statistic)
bb_cp0[l]<-(bb_l[l]<standard_statistic)&(bb_r[l]>standard_statistic)
per_cp0[l]<-(per_l[l]<standard_statistic)&(per_r[l]>standard_statistic)
}
#standard normal 绘图
max_lim<-max(sn_r-sn_l)
xx<-c(1:200)
yy<-c(1:200)
plot(xx,yy,type="n",ylim = c(standard_statistic-max_lim,standard_statistic+max_lim),xlab = "区间序号",ylab = "置信区间",main = "置信区间(standard normal)是否覆盖标准统计量值")
abline(h=standard_statistic,lwd=1,col=2)
for(i in 1:200){
segments(i,sn_l[i],i,sn_r[i],col=2-sn_cp0[i],lwd=2-sn_cp0[i])
}
#计算覆盖率
sn_cp<-mean(sn_cp0)
bb_cp<-mean(bb_cp0)
per_cp<-mean(per_cp0)
return(list(sn_cp=sn_cp,
bb_cp=bb_cp,
per_cp=per_cp))
}
例:第二类函数返回值示意:
BOOT_COVER_PROB_NORM(size0=100,mean0=1,sd0=2,statistic=mean,standard_statistic = 1)
## $sn_cp
## [1] 0.948
##
## $bb_cp
## [1] 0.938
##
## $per_cp
## [1] 0.948
二项分布的第二类函数
#bootstrap求覆盖率函数(针对二项分布)
#BOOT_COVER_PROB_BIN的参数如下:size0=样本量,p=二项分布参数,statistic=统计量,standard_statistic=标准统计量值,B=bootstrap次数,alpha=显著性水平
BOOT_COVER_PROB_BIN<-function(size0=100,p=0.5,statistic=mean,standard_statistic=0.5,B=500,alpha=0.05){
sn_cp0<-c();bb_cp0<-c();per_cp0<-c()
sn_cp<-0;bb_cp<-0;per_cp<-0
sn_l<-c();sn_r<-c();bb_l<-c();bb_r<-c();per_l<-c();per_r<-c()
for(l in 1:500){
x<-rbinom(size0,1,p)
statistic_hat<-statistic(x)
statistic_hat_boot<-c(1:B)#对statistic_hat做bootstrap得到的一组statistic值
for(i in 1:B){
indices<-sample(1:size0,size0,replace = TRUE)#有放回抽样
x_bootstrap<-x[indices]
statistic_hat_boot[i]<-statistic(x_bootstrap)
}
# the standard normal bootstrap confidence interval
sd_hat_forstatistic<-sd(statistic_hat_boot)
sn_l[l]<-statistic_hat+qnorm(alpha/2)*sd_hat_forstatistic
sn_r[l]<-statistic_hat-qnorm(alpha/2)*sd_hat_forstatistic
# the basic bootstrap confidence interval
bb_l[l]<-2*statistic_hat-quantile(statistic_hat_boot,1-alpha/2)
bb_r[l]<-2*statistic_hat-quantile(statistic_hat_boot,alpha/2)
# the percentile confidence interval
per_l[l]<-quantile(statistic_hat_boot,alpha/2)
per_r[l]<-quantile(statistic_hat_boot,1-alpha/2)
#存储是否覆盖数组(由于总体的statistic数据未知,这里我们用原样本的statistic数据代替之)
sn_cp0[l]<-(sn_l[l]<standard_statistic)&(sn_r[l]>standard_statistic)
bb_cp0[l]<-(bb_l[l]<standard_statistic)&(bb_r[l]>standard_statistic)
per_cp0[l]<-(per_l[l]<standard_statistic)&(per_r[l]>standard_statistic)
}
#standard normal 绘图
max_lim<-max(sn_r-sn_l)
xx<-c(1:200)
yy<-c(1:200)
plot(xx,yy,type="n",ylim = c(standard_statistic-max_lim,standard_statistic+max_lim),xlab = "区间序号",ylab = "置信区间",main = "置信区间(standard normal)是否覆盖标准统计量值")
abline(h=standard_statistic,lwd=1,col=2)
for(i in 1:200){
segments(i,sn_l[i],i,sn_r[i],col=2-sn_cp0[i],lwd=2-sn_cp0[i])
}
#计算覆盖率
sn_cp<-mean(sn_cp0)
bb_cp<-mean(bb_cp0)
per_cp<-mean(per_cp0)
return(list(sn_cp=sn_cp,
bb_cp=bb_cp,
per_cp=per_cp))
}
在本阶段我们侧重想要观察的变量是bootstrap次数、样本量、置信区间类型、样本分布的参数,故我们选择侧重研究二项分布,因其涉及参数较少,或许可以更直观的显现出某种性质,且二项分布在生产生活中应用较为广泛,或许可以有更积极地现实意义与价值。
在固定同一组样本的前提下,我们分别对三种置信区间使用了1-1000中较有代表性的一些数值作为bootstrap次数,来观察覆盖率与bootstrap次数的关系。
结论:1、以这种算法计算出的覆盖率会随着B的增大而趋向于置信水平(90%),基本在大于200次之后就已经较好的接近了置信水平,且我们可以发现这并不是严格的单调递增,而是有一定的不稳定性的;
2、在B较小时经验覆盖数值整体偏小,这种偏差不能忽略,且需要引起重视。Normal CI偏差最小,Percentile CI偏差最大,甚至已经低于0.6,偏差很大。这也提醒我们,在做这一类覆盖率研究时,bootstrap次数不能选的太小,选的太小会与置信水平有很大的差距,严重影响了结果的可信度。
我们简单地想象,或许样本量越大,数据越稳定,覆盖率越接近置信水平。但我们通过测试可以发现,这种定义下的覆盖率与样本量n呈现不出明显的相关性(normal basic percentile均是如此)。且覆盖率集中在置信水平左右(有相当一部分数据超出置信水平,这与后面将要讨论的第二种定义有很大的区别)。这里波动性比较大是因为这里的bootstrap次数较小。
下图分别为normal,basic和percentile
与和样本量关系的讨论相似地,这里也没有显现出与参数p的相关性,且数据分布在置信水平左右。值得注意地,从图像上可以看出,在p较小或较大时覆盖率集中在100%处,这是由于样本过于集中所带来的误差。
下图分别为normal,basic和percentile
综上所述,我们可以发现,在这种覆盖率的定义下(一组CI多组统计量),在bootstrap次数保证的前提下,覆盖率平凡地接近于置信水平,这是符合我们直观想象的。
但这种覆盖率展现出的性质匮乏,相关性很低,研究价值不高,没有展现出分布的内蕴性质。或许在数据模拟估计上这种定义是有实用性的,但在内在性质研究方面我们没有看到这种定义方法带来的好处与思考价值。
第三章的代码实现如下:
#二项分布,standard normal
#参数p
#n=50
p<-c()
N<-50
fugailv_sn<-c()
fugailv_bb<-c()
fugailv_per<-c()
for(i in 1:1000){
y0<-rbinom(N,1,i/1001)
A<-BOOT(x=y0,statistic = mean,B=100,alpha=0.05)
fugailv_sn[i]<-1-A[1,5]
fugailv_bb[i]<-1-A[2,5]
fugailv_per[i]<-1-A[3,5]
}
par(mfrow = c(1,1))
plot(c(1:1000)/1001,fugailv_sn,type ="p",xlab="二项分布中的p",ylab="经验覆盖比率",main="二项分布(alpha=0.05,n=50,p=0~1)")
abline(a=0.95,b=0,col="red")
plot(c(1:1000)/1001,fugailv_bb,type ="p",xlab="二项分布中的p",ylab="经验覆盖比率",main="二项分布(alpha=0.05,n=50,p=0~1)")
abline(a=0.95,b=0,col="red")
plot(c(1:1000)/1001,fugailv_per,type ="p",xlab="二项分布中的p",ylab="经验覆盖比率",main="二项分布(alpha=0.05,n=50,p=0~1)")
abline(a=0.95,b=0,col="red")
#二项分布,standard normal
#样本量n
#n=1~100,p=0.5
p<-0.5
N<-1000
fugailv_sn<-c()
fugailv_bb<-c()
fugailv_per<-c()
for(n in 1:N){
y0<-rbinom(n,1,p)
A<-BOOT(x=y0,statistic = mean,B=200,alpha=0.05)
fugailv_sn[n]<-1-A[1,5]
fugailv_bb[n]<-1-A[2,5]
fugailv_per[n]<-1-A[3,5]
}
par(mfrow = c(1,1))
plot(c(1:N),fugailv_sn,type ="p",xlab="二项分布中的n",ylab="经验覆盖比率",main="二项分布(alpha=0.05,p=0.5,n=1~1000)")
abline(a=0.95,b=0,col="red")
plot(c(1:N),fugailv_bb,type ="p",xlab="二项分布中的n",ylab="经验覆盖比率",main="二项分布(alpha=0.05,p=0.5,n=1~1000)")
abline(a=0.95,b=0,col="red")
plot(c(1:N),fugailv_per,type ="p",xlab="二项分布中的n",ylab="经验覆盖比率",main="二项分布(alpha=0.05,p=0.5,n=1~1000)")
abline(a=0.95,b=0,col="red")
#standard normal CI 经验覆盖与B的关系(二项分布 p=0.5 n=100)
y0<-rbinom(500,1,0.5)
bb<-c()
bb5<-c();bb10<-c();bb20<-c();bb50<-c();bb200<-c();bb1000<-c();bb5000<-c()
for(i in 1:10){
bb5[i]<-1-BOOT(x=y0,statistic = mean,B=5,alpha = 0.1)[1,5]
}
bb[1]<-mean(bb5)
for(i in 1:10){
bb10[i]<-1-BOOT(x=y0,statistic = mean,B=10,alpha = 0.1)[1,5]
}
bb[2]<-mean(bb10)
for(i in 1:10){
bb20[i]<-1-BOOT(x=y0,statistic = mean,B=20,alpha = 0.1)[1,5]
}
bb[3]<-mean(bb20)
for(i in 1:10){
bb50[i]<-1-BOOT(x=y0,statistic = mean,B=50,alpha = 0.1)[1,5]
}
bb[4]<-mean(bb50)
for(i in 1:10){
bb200[i]<-1-BOOT(x=y0,statistic = mean,B=200,alpha = 0.1)[1,5]
}
bb[5]<-mean(bb200)
for(i in 1:5){
bb1000[i]<-1-BOOT(x=y0,statistic = mean,B=1000,alpha = 0.1)[1,5]
}
bb[6]<-mean(bb1000)
par(mfrow = c(1,1))
ind<-c(5,10,20,50,200,1000)
plot(ind,bb,type = "b",ylim=c(0.7,1),xlab="bootstrap次数B",ylab="覆盖率",main="standard normal经验覆盖与bootstrap次数B的关系图(二项分布)")
#basic CI 经验覆盖与B的关系(二项分布 p=0.5 n=100)
y0<-rbinom(500,1,0.5)
bb<-c()
bb5<-c();bb10<-c();bb20<-c();bb50<-c();bb200<-c();bb1000<-c();bb5000<-c()
for(i in 1:10){
bb5[i]<-1-BOOT(x=y0,statistic = mean,B=5,alpha = 0.1)[2,5]
}
bb[1]<-mean(bb5)
for(i in 1:10){
bb10[i]<-1-BOOT(x=y0,statistic = mean,B=10,alpha = 0.1)[2,5]
}
bb[2]<-mean(bb10)
for(i in 1:10){
bb20[i]<-1-BOOT(x=y0,statistic = mean,B=20,alpha = 0.1)[2,5]
}
bb[3]<-mean(bb20)
for(i in 1:10){
bb50[i]<-1-BOOT(x=y0,statistic = mean,B=50,alpha = 0.1)[2,5]
}
bb[4]<-mean(bb50)
for(i in 1:10){
bb200[i]<-1-BOOT(x=y0,statistic = mean,B=200,alpha = 0.1)[2,5]
}
bb[5]<-mean(bb200)
for(i in 1:5){
bb1000[i]<-1-BOOT(x=y0,statistic = mean,B=1000,alpha = 0.1)[2,5]
}
bb[6]<-mean(bb1000)
par(mfrow = c(1,1))
ind<-c(5,10,20,50,200,1000)
plot(ind,bb,type = "b",ylim=c(0.7,1),xlab="bootstrap次数B",ylab="覆盖率",main="basic经验覆盖与bootstrap次数B的关系图(二项分布)")
#percentile CI 经验覆盖与B的关系(二项分布 p=0.5 n=100)
y0<-rbinom(500,1,0.5)
bb<-c()
bb5<-c();bb10<-c();bb20<-c();bb50<-c();bb200<-c();bb1000<-c();bb5000<-c()
for(i in 1:10){
bb5[i]<-1-BOOT(x=y0,statistic = mean,B=5,alpha = 0.1)[3,5]
}
bb[1]<-mean(bb5)
for(i in 1:10){
bb10[i]<-1-BOOT(x=y0,statistic = mean,B=10,alpha = 0.1)[3,5]
}
bb[2]<-mean(bb10)
for(i in 1:10){
bb20[i]<-1-BOOT(x=y0,statistic = mean,B=20,alpha = 0.1)[3,5]
}
bb[3]<-mean(bb20)
for(i in 1:10){
bb50[i]<-1-BOOT(x=y0,statistic = mean,B=50,alpha = 0.1)[3,5]
}
bb[4]<-mean(bb50)
for(i in 1:10){
bb200[i]<-1-BOOT(x=y0,statistic = mean,B=200,alpha = 0.1)[3,5]
}
bb[5]<-mean(bb200)
for(i in 1:5){
bb1000[i]<-1-BOOT(x=y0,statistic = mean,B=1000,alpha = 0.1)[3,5]
}
bb[6]<-mean(bb1000)
par(mfrow = c(1,1))
ind<-c(5,10,20,50,200,1000)
plot(ind,bb,type = "b",ylim=c(0.5,1),xlab="bootstrap次数B",ylab="覆盖率",main="percentile经验覆盖与bootstrap次数B的关系图(二项分布)")
在本阶段我们侧重想要观察的变量是bootstrap次数、样本量、置信区间类型、样本分布的参数,我们选择的分布是二项分布与正态分布,希望可以通过这种定义方法找到关于分布的内蕴性质。(这里我们主要讨论的统计量为均值)
a、二项分布(normal、basic、percentile)
这里的研究方法与上面第一类函数相似。
结论:1、以这种算法计算出的覆盖率会随着B的增大而趋向于置信水平(90%),基本在大于200次之后就已经较好的接近了置信水平,且我们可以发现这并不是严格的单调递增,而是有一定的不稳定性的;
2、在B较小时经验覆盖数值整体偏小,这种偏差不能忽略,且需要引起重视。Normal CI偏差最小,Percentile CI偏差最大。
b、正态分布(normal、basic、percentile)
这里得到的结果几乎与上面完全相同,故不再赘述。
综上,我们可以对bootstrap次数有一个比较广泛的结论:在bootstrap次数较小时(小于100),进行的bootstrap覆盖率计算都会偏离置信水平,次数越小偏差越大;其中以normal偏差最小,percentile偏差最大。这个结论在二项、正态样本均是成立的,在两种覆盖率定义下也均是成立的。这也提醒我们,在做bootstrap覆盖率研究时,次数不能选的太小。
在之前的第一类函数定义中,我们并没有观察到覆盖率与样本量之间的相关性,但在多组CI的覆盖率定义下,我们观察到了不平凡的现象。
a、二项分布(normal、basic、percentile)(p=0.005)
Normal CI:
我们从图中可以发现,,样本量在0-400左右时,覆盖率呈现出单增趋势。(这是我们在之前的定义中没有的现象)并且在非常接近置信水平后出现了明显下降趋势,紧接着出现连续上升下降的振荡效果,且随着样本量的增加,震荡幅度越来越小。可以想象,在更大的样本量情况下,振荡效果会变小,整体趋于置信水平。
另外,我们从图中发现没有点超出置信水平0.95,在这组测试中最大值在n=889出现,覆盖率为94.8%。
Percentile CI:
与normal CI类似地,我们也发现了不平凡的振荡效果。如果我们仔细观察振荡点的分布,这里与normal有一点点细微的差异:normal中点的下降似乎是平滑下降,而percentile中的下降更像是突然下降,这是一个值得思考的现象。
另外,我们从图中发现仅有一个点超出0.95,在这组测试中最大值在n=924出现,覆盖率为95.2%。
Basic CI:
比起前两张图的振荡效果,basic区间的振荡更加令人吃惊。这里并不像以上两种区间,在达到置信水平时才有下降,而是更早地与置信水平有一定距离时便出现了下跌,且下跌迅速,断层很大。整体上仍呈现出趋向于置信水平的效果。
另外,我们从图中发现没有点超出0.95,在这组测试中最大值在n=1169出现,覆盖率为93.6%。
b、二项分布(normal、basic、percentile)(p=0.5)
我们或许会注意到,上面的测试中二项分布的p取的很小,自然地我们思考是否因为p值过小而造成了振荡呢?这里我们又对p=0.5进行了测试,如我们所料,振荡效果消失了,覆盖率非常接近与置信水平。
a、样本量n=20
从上面的讨论我们发现关于样本量的振荡效果与p的选取有密切联系,自然地想到关于参数p的覆盖率是否会同样出现不平凡的效果。下图为p在0-1之间覆盖率的图像,选取n=20(仅以normal CI为例)。
由图像我们发现,振荡效果关于参数p再次出现。P越偏离0.5,覆盖率越偏离置信水平,振荡效果越明显。P很接近0.5时振荡效果不明显,覆盖率与置信水平较为接近。
另外,我们发现仅有一个点超出0.99,在这组测试中最大值在p=0.33出现,覆盖率为99.2%
b、样本量n=100
由上面B的讨论我们知道,n较小时对比较极端的p值覆盖率有偏差,这也是我们在上面n=20的情况下看到的,产生了偏差和震荡现象。自然地,我们将n设置为100,如下图,我们发现振荡效果变得非常不明显,几乎均很好地趋向于显著性水平。
二项分布中若固定p值:
1、p值较小或较大时,覆盖率关于样本量振荡明显(且以basic CI最为显著)
2、p值较接近0.5时,振荡效果几乎消失
二项分布中若固定n:
1、n较小时振荡效果明显(振荡出现在p较小或较大的区间)
2、n较大时振荡效果几乎消失
综上,我们可以观察到,振荡效果的出现范围基本满足两种条件:p值较小、较大或n值较小。我们可以合理的猜测,这是在参数较为偏离正常情况时出现的一种不稳定性质。
查询资料后得知:这种现象主要来源于二项分布的潜在离散性和潜在偏斜型。
注意:在上面的研究过程中,还发现了一个有趣的现象,在多数的实验中,均出现了覆盖率都是趋近于置信水平,但达不到置信水平的现象。这其中的具体原理我们恕本次研究时间有限,未能展开探讨,我们猜测这种现象或许是由区间类型造成的,这是否也提醒我们可以构造出更好的区间类型,是一个研究与思考的方向。
下面我们来探讨一下关于正态分布中覆盖率展现出来的性质
下图为正态分布均值与覆盖率的关系,我们可以发现覆盖率很好地趋向于置信水平,并未出现振荡效果。(basic percentile)
下图为正态分布样本量与覆盖率的关系,可以看出与之前相似地呈现出单增趋势,但并未体现出振荡效果。(basic percentile)
我们或许可以猜测,正态分布中并没有振荡现象。关于振荡现象是否为二项分布所独有这一命题,或许需要更严谨更深入的理论支持,在本次研究中无法展开。不过我们可以猜测,这是二项所独有的,因为振荡效果主要源于偏斜性与离散性,这是二项分布带来的。
以下为第四章的代码实现:
#以下为二项分布
#样本量n
#standard normal
#二项分布,p=0.005,alpha=0.05,B=200
size_bin<-c()
for(i in 1:200){
size_bin[i]<-BOOT_COVER_PROB_BIN(size0 =i*7,p=0.005,statistic = mean,standard_statistic=0.005,B=200,alpha = 0.05)[[1]]
}
par(mfrow = c(1,1))
plot(c(1:200)*7,c(1:200),ylim = c(0,1),type = "n",xlab = "样本量n",ylab="覆盖率",main="二项分布覆盖率与样本量的关系")
points(c(1:200)*7,size_bin)
abline(a=0.95,b=0,col="red")
#发现没有点超出0.95,在这组测试中最大值在n=889出现,覆盖率为94.8%
#以下:对上述的离散点进行多项式拟合(degree=25)
y<-size_bin
y0<-as.numeric(y)
x0<-c(1:200)*7
x0<-data.frame(x0)
data0<-cbind(x0,y0)
lm1<-lm(y0~poly(x0,25),data = data0)
x0<-c(1:1400)
y0hat= predict(lm1, data.frame(x=x0))
plot(c(1:1400),c(1:1400),ylim = c(0,1),type = "n",xlab = "样本量n",ylab="覆盖率",main="二项分布覆盖率与样本量的多项式拟合")
lines(x0, y0hat, col=3, lwd=2)
abline(a=0.95,b=0,col="red")
#震荡效果明显
#样本量n
#standard normal
#二项分布,p=0.5,alpha=0.05,B=200
size_bin<-c()
for(i in 1:200){
size_bin[i]<-BOOT_COVER_PROB_BIN(size0 =i*7,p=0.5,statistic = mean,standard_statistic=0.5,B=200,alpha = 0.05)[[1]]
}
par(mfrow = c(1,1))
plot(c(1:200)*7,c(1:200),ylim = c(0,1),type = "n",xlab = "样本量n",ylab="覆盖率",main="二项分布覆盖率与样本量的关系 p=0.5")
lines(c(1:200)*7,size_bin)
abline(a=0.95,b=0,col="red")
#这时p=0.5,对称性较好,没有出现明显震荡效果
#样本量n
#percentile!
#二项分布,p=0.005,alpha=0.05,B=200
size_bin<-c()
for(i in 1:200){
size_bin[i]<-BOOT_COVER_PROB_BIN(size0 =i*7,p=0.005,statistic = mean,standard_statistic=0.005,B=200,alpha = 0.05)[[3]]
}
par(mfrow = c(1,1))
plot(c(1:200)*7,c(1:200),ylim = c(0,1),type = "n",xlab = "样本量n",ylab="覆盖率",main="二项分布覆盖率与样本量的关系 percentile")
points(c(1:200)*7,size_bin)
abline(a=0.95,b=0,col="red")
#发现仅有一个点超出0.95,在这组测试中最大值在n=924出现,覆盖率为95.2%
#以下:对上述的离散点进行多项式拟合(degree=25)
y<-size_bin
y0<-as.numeric(y)
x0<-c(1:200)*7
x0<-data.frame(x0)
data0<-cbind(x0,y0)
lm1<-lm(y0~poly(x0,25),data = data0)
x0<-c(1:1400)
y0hat= predict(lm1, data.frame(x=x0))
plot(c(1:1400),c(1:1400),ylim = c(0,1),type = "n",xlab = "样本量n",ylab="覆盖率",main="二项分布覆盖率与样本量的多项式拟合 percentile")
lines(x0, y0hat, col=3, lwd=2)
abline(a=0.95,b=0,col="red")
#震荡效果明显
#样本量n
#basic!
#二项分布,p=0.005,alpha=0.05,B=200
size_bin<-c()
for(i in 1:200){
size_bin[i]<-BOOT_COVER_PROB_BIN(size0 =i*7,p=0.005,statistic = mean,standard_statistic=0.005,B=200,alpha = 0.05)[[2]]
}
par(mfrow = c(1,1))
plot(c(1:200)*7,c(1:200),ylim = c(0,1),type = "n",xlab = "样本量n",ylab="覆盖率",main="二项分布覆盖率与样本量的关系 basic")
points(c(1:200)*7,size_bin)
abline(a=0.95,b=0,col="red")
#发现没有点超出0.95,在这组测试中最大值在n=1169出现,覆盖率为93.6%;且从趋势来看样本量越大一些峰值点将越接近0.95
#震荡效果明显
#参数p
#standard normal,alpha=0.01,B=200
#二项分布,样本量100
p_bin<-c()
for(i in 1:200){
p_bin[i]<-BOOT_COVER_PROB_BIN(size0 =100,p=i/201,statistic = mean,standard_statistic=i/201,B=200,alpha = 0.01)[[1]]
}
par(mfrow = c(1,1))
plot(c(1:200)/201,c(1:200),ylim = c(0,1),type = "n",xlab = "二项分布参数p",ylab="覆盖率",main="二项分布覆盖率与参数p的关系")
points(c(1:200)/201,p_bin)
abline(a=0.99,b=0,col="red")
#以下:对上述的离散点进行多项式拟合(degree=24)
y<-p_bin
y0<-as.numeric(y)
x0<-c(1:200)/201
x0<-data.frame(x0)
data0<-cbind(x0,y0)
lm1<-lm(y0~poly(x0,24),data = data0)
x0<-seq(0,1,0.001)
y0hat= predict(lm1, data.frame(x=x0))
plot(seq(0,1,0.001),seq(0,1,0.001),ylim = c(0,1),type = "n",xlab = "二项分布参数p",ylab="覆盖率",main="二项分布覆盖率与参数p的多项式拟合")
lines(x0, y0hat, col=3, lwd=2)
abline(a=0.99,b=0,col="red")
#震荡效果不明显
#参数p
#standard normal
#二项分布,样本量20,alpha=0.01,B=200
p_bin<-c()
for(i in 1:200){
p_bin[i]<-BOOT_COVER_PROB_BIN(size0 =20,p=i/201,statistic = mean,standard_statistic=i/201,B=200,alpha = 0.01)[[1]]
}
par(mfrow = c(1,1))
plot(c(1:200)/201,c(1:200),ylim = c(0,1),type = "n",xlab = "二项分布参数p",ylab="覆盖率",main="二项分布覆盖率与参数p的关系")
points(c(1:200)/201,p_bin)
abline(a=0.99,b=0,col="red")
#发现仅有一个点超出0.99,在这组测试中最大值在p=0.33出现,覆盖率为99.2%(which(p_bin>0.99))
#以下:对上述的离散点进行多项式拟合(degree=24)
y<-p_bin
y0<-as.numeric(y)
x0<-c(1:200)/201
x0<-data.frame(x0)
data0<-cbind(x0,y0)
lm1<-lm(y0~poly(x0,24),data = data0)
x0<-seq(0,1,0.001)
y0hat= predict(lm1, data.frame(x=x0))
plot(seq(0,1,0.001),seq(0,1,0.001),ylim = c(0,1),type = "n",xlab = "二项分布参数p",ylab="覆盖率",main="二项分布覆盖率与参数p的多项式拟合")
lines(x0, y0hat, col=3, lwd=2)
abline(a=0.99,b=0,col="red")
#震荡效果明显
#以下为正态总体
#参数mean
#percentile
#正态分布,sd0=1,alpha=0.05,B=200
mean_norm<-c()
for(i in 1:100){
mean_norm[i]<-BOOT_COVER_PROB_NORM(size0 =20,mean0 = i/50,sd0=1,statistic = mean,standard_statistic=i/50,B=200,alpha = 0.05)[[3]]
}
par(mfrow = c(1,1))
plot(c(1:100)/50,c(1:100),ylim = c(0,1),type = "n",xlab = "均值",ylab="覆盖率",main="正态分布覆盖率与参数均值的关系 percentile")
points(c(1:100)/50,mean_norm)
abline(a=0.95,b=0,col="red")
#没有点超出0.95,最大值0.94在均值为0.08和0.3处达到
#无震荡
#参数mean
#basic
#正态分布,sd0=1,alpha=0.05,B=200
mean_norm<-c()
for(i in 1:100){
mean_norm[i]<-BOOT_COVER_PROB_NORM(size0 =20,mean0 = i/50,sd0=1,statistic = mean,standard_statistic=i/50,B=200,alpha = 0.05)[[2]]
}
par(mfrow = c(1,1))
plot(c(1:100)/50,c(1:100),ylim = c(0,1),type = "n",xlab = "均值",ylab="覆盖率",main="正态分布覆盖率与参数均值的关系 basic")
points(c(1:100)/50,mean_norm)
abline(a=0.95,b=0,col="red")
#仅有一个点恰达到0.95,均值为1.7
#无震荡
#样本量n
#basic
#正态分布,mean=0,sd=1,alpha=0.05,B=200
size_norm<-c()
for(i in 1:50){
size_norm[i]<-BOOT_COVER_PROB_NORM(size0 =i,mean0 = 0,sd0 = 1,statistic = mean,standard_statistic=0,B=200,alpha = 0.05)[[2]]
}
par(mfrow = c(1,1))
plot(c(1:50),c(1:50),ylim = c(0,1),type = "n",xlab = "样本量n",ylab="覆盖率",main="正态分布覆盖率与样本量的关系 basic")
points(c(1:50),size_norm)
abline(a=0.95,b=0,col="red")
#发现没有一个点超出0.95,在这组测试中最大值在n=37出现,覆盖率为94.8%
#以下:对上述的离散点进行多项式拟合(degree=25)
y<-size_norm
y0<-as.numeric(y)
x0<-c(1:50)
x0<-data.frame(x0)
data0<-cbind(x0,y0)
lm1<-lm(y0~poly(x0,25),data = data0)
x0<-c(1:50)
y0hat= predict(lm1, data.frame(x=x0))
plot(c(1:50),c(1:50),ylim = c(0,1),type = "n",xlab = "样本量n",ylab="覆盖率",main="正态分布覆盖率与样本量的多项式拟合 basic")
lines(x0, y0hat, col=3, lwd=2)
abline(a=0.95,b=0,col="red")
#无震荡效果
#样本量n
#percentile
#正态分布,mean=0,sd=1,alpha=0.05,B=200
size_norm<-c()
for(i in 1:50){
size_norm[i]<-BOOT_COVER_PROB_NORM(size0 =i,mean0 = 0,sd0 = 1,statistic = mean,standard_statistic=0,B=200,alpha = 0.05)[[3]]
}
par(mfrow = c(1,1))
plot(c(1:50),c(1:50),ylim = c(0,1),type = "n",xlab = "样本量n",ylab="覆盖率",main="正态分布覆盖率与样本量的关系 percentile")
points(c(1:50),size_norm)
abline(a=0.95,b=0,col="red")
#发现仅有一个点超出0.95,在这组测试中最大值在n=40出现,覆盖率为95.2%
#以下:对上述的离散点进行多项式拟合(degree=25)
y<-size_norm
y0<-as.numeric(y)
x0<-c(1:50)
x0<-data.frame(x0)
data0<-cbind(x0,y0)
lm1<-lm(y0~poly(x0,25),data = data0)
x0<-c(1:50)
y0hat= predict(lm1, data.frame(x=x0))
plot(c(1:50),c(1:50),ylim = c(0,1),type = "n",xlab = "样本量n",ylab="覆盖率",main="正态分布覆盖率与样本量的多项式拟合 percentile")
lines(x0, y0hat, col=3, lwd=2)
abline(a=0.95,b=0,col="red")
#无震荡效果
#basic bootstrap CI 覆盖率与B的关系(正态分布)
bb<-c()
bb5<-c();bb10<-c();bb20<-c();bb50<-c();bb200<-c();bb1000<-c();bb5000<-c()
for(i in 1:10){
bb5[i]<-BOOT_COVER_PROB_NORM(size0 = 1000,mean0=5,sd0=5,statistic = mean,standard_statistic = 5,B=5,alpha = 0.1)[[2]]
}
bb[1]<-mean(bb5)
for(i in 1:10){
bb10[i]<-BOOT_COVER_PROB_NORM(size0 = 1000,mean0=5,sd0=5,statistic = mean,standard_statistic = 5,B=10,alpha = 0.1)[[2]]
}
bb[2]<-mean(bb10)
for(i in 1:10){
bb20[i]<-BOOT_COVER_PROB_NORM(size0 = 1000,mean0=5,sd0=5,statistic = mean,standard_statistic = 5,B=20,alpha = 0.1)[[2]]
}
bb[3]<-mean(bb20)
for(i in 1:10){
bb50[i]<-BOOT_COVER_PROB_NORM(size0 = 1000,mean0=5,sd0=5,statistic = mean,standard_statistic = 5,B=50,alpha = 0.1)[[2]]
}
bb[4]<-mean(bb50)
for(i in 1:10){
bb200[i]<-BOOT_COVER_PROB_NORM(size0 = 1000,mean0=5,sd0=5,statistic = mean,standard_statistic = 5,B=200,alpha = 0.1)[[2]]
}
bb[5]<-mean(bb200)
for(i in 1:5){
bb1000[i]<-BOOT_COVER_PROB_NORM(size0 = 1000,mean0=5,sd0=5,statistic = mean,standard_statistic = 5,B=1000,alpha = 0.1)[[2]]
}
bb[6]<-mean(bb1000)
ind<-c(5,10,20,50,200,1000)
plot(ind,bb,type = "b",ylim=c(0.6,1),xlab="bootstrap次数B",ylab="覆盖率",main="basic区间覆盖率与bootstrap次数B的关系图(正态分布)")
#basic bootstrap CI 覆盖率与B的关系(二项分布)
bb<-c()
bb5<-c();bb10<-c();bb20<-c();bb50<-c();bb200<-c();bb1000<-c();bb5000<-c()
for(i in 1:10){
bb5[i]<-BOOT_COVER_PROB_BIN(size0 = 1000,p=0.5,statistic = mean,standard_statistic = 0.5,B=5,alpha = 0.1)[[2]]
}
bb[1]<-mean(bb5)
for(i in 1:10){
bb10[i]<-BOOT_COVER_PROB_BIN(size0 = 1000,p=0.5,statistic = mean,standard_statistic = 0.5,B=10,alpha = 0.1)[[2]]
}
bb[2]<-mean(bb10)
for(i in 1:10){
bb20[i]<-BOOT_COVER_PROB_BIN(size0 = 1000,p=0.5,statistic = mean,standard_statistic = 0.5,B=20,alpha = 0.1)[[2]]
}
bb[3]<-mean(bb20)
for(i in 1:10){
bb50[i]<-BOOT_COVER_PROB_BIN(size0 = 1000,p=0.5,statistic = mean,standard_statistic = 0.5,B=50,alpha = 0.1)[[2]]
}
bb[4]<-mean(bb50)
for(i in 1:10){
bb200[i]<-BOOT_COVER_PROB_BIN(size0 = 1000,p=0.5,statistic = mean,standard_statistic = 0.5,B=200,alpha = 0.1)[[2]]
}
bb[5]<-mean(bb200)
for(i in 1:5){
bb1000[i]<-BOOT_COVER_PROB_BIN(size0 = 1000,p=0.5,statistic = mean,standard_statistic = 0.5,B=1000,alpha = 0.1)[[2]]
}
bb[6]<-mean(bb1000)
ind<-c(5,10,20,50,200,1000)
plot(ind,bb,type = "b",ylim=c(0.6,1),xlab="bootstrap次数B",ylab="覆盖率",main="basic区间覆盖率与bootstrap次数B的关系图(二项分布)")
#percentile CI 覆盖率与B的关系(正态分布)
bb<-c()
bb5<-c();bb10<-c();bb20<-c();bb50<-c();bb200<-c();bb1000<-c();bb5000<-c()
for(i in 1:10){
bb5[i]<-BOOT_COVER_PROB_NORM(size0 = 1000,mean0=5,sd0=5,statistic = mean,standard_statistic = 5,B=5,alpha = 0.1)[[3]]
}
bb[1]<-mean(bb5)
for(i in 1:10){
bb10[i]<-BOOT_COVER_PROB_NORM(size0 = 1000,mean0=5,sd0=5,statistic = mean,standard_statistic = 5,B=10,alpha = 0.1)[[3]]
}
bb[2]<-mean(bb10)
for(i in 1:10){
bb20[i]<-BOOT_COVER_PROB_NORM(size0 = 1000,mean0=5,sd0=5,statistic = mean,standard_statistic = 5,B=20,alpha = 0.1)[[3]]
}
bb[3]<-mean(bb20)
for(i in 1:10){
bb50[i]<-BOOT_COVER_PROB_NORM(size0 = 1000,mean0=5,sd0=5,statistic = mean,standard_statistic = 5,B=50,alpha = 0.1)[[3]]
}
bb[4]<-mean(bb50)
for(i in 1:10){
bb200[i]<-BOOT_COVER_PROB_NORM(size0 = 1000,mean0=5,sd0=5,statistic = mean,standard_statistic = 5,B=200,alpha = 0.1)[[3]]
}
bb[5]<-mean(bb200)
for(i in 1:5){
bb1000[i]<-BOOT_COVER_PROB_NORM(size0 = 1000,mean0=5,sd0=5,statistic = mean,standard_statistic = 5,B=1000,alpha = 0.1)[[3]]
}
bb[6]<-mean(bb1000)
ind<-c(5,10,20,50,200,1000)
plot(ind,bb,type = "b",ylim=c(0.6,1),xlab="bootstrap次数B",ylab="覆盖率",main="percenti区间覆盖率与bootstrap次数B的关系图(正态分布)")
#percentile CI 覆盖率与B的关系(二项分布)
bb<-c()
bb5<-c();bb10<-c();bb20<-c();bb50<-c();bb200<-c();bb1000<-c();bb5000<-c()
for(i in 1:10){
bb5[i]<-BOOT_COVER_PROB_BIN(size0 = 1000,p=0.5,statistic = mean,standard_statistic = 0.5,B=5,alpha = 0.1)[[3]]
}
bb[1]<-mean(bb5)
for(i in 1:10){
bb10[i]<-BOOT_COVER_PROB_BIN(size0 = 1000,p=0.5,statistic = mean,standard_statistic = 0.5,B=10,alpha = 0.1)[[3]]
}
bb[2]<-mean(bb10)
for(i in 1:10){
bb20[i]<-BOOT_COVER_PROB_BIN(size0 = 1000,p=0.5,statistic = mean,standard_statistic = 0.5,B=20,alpha = 0.1)[[3]]
}
bb[3]<-mean(bb20)
for(i in 1:10){
bb50[i]<-BOOT_COVER_PROB_BIN(size0 = 1000,p=0.5,statistic = mean,standard_statistic = 0.5,B=50,alpha = 0.1)[[3]]
}
bb[4]<-mean(bb50)
for(i in 1:10){
bb200[i]<-BOOT_COVER_PROB_BIN(size0 = 1000,p=0.5,statistic = mean,standard_statistic = 0.5,B=200,alpha = 0.1)[[3]]
}
bb[5]<-mean(bb200)
for(i in 1:5){
bb1000[i]<-BOOT_COVER_PROB_BIN(size0 = 1000,p=0.5,statistic = mean,standard_statistic = 0.5,B=1000,alpha = 0.1)[[3]]
}
bb[6]<-mean(bb1000)
ind<-c(5,10,20,50,200,1000)
plot(ind,bb,type = "b",ylim=c(0.6,1),xlab="bootstrap次数B",ylab="覆盖率",main="percentile区间覆盖率与bootstrap次数B的关系图(二项分布)")
#standard normal CI 覆盖率与B的关系(二项分布)
bb<-c()
bb5<-c();bb10<-c();bb20<-c();bb50<-c();bb200<-c();bb1000<-c();bb5000<-c()
for(i in 1:10){
bb5[i]<-BOOT_COVER_PROB_BIN(size0 = 1000,p=0.5,statistic = mean,standard_statistic = 0.5,B=5,alpha = 0.1)[[1]]
}
bb[1]<-mean(bb5)
for(i in 1:10){
bb10[i]<-BOOT_COVER_PROB_BIN(size0 = 1000,p=0.5,statistic = mean,standard_statistic = 0.5,B=10,alpha = 0.1)[[1]]
}
bb[2]<-mean(bb10)
for(i in 1:10){
bb20[i]<-BOOT_COVER_PROB_BIN(size0 = 1000,p=0.5,statistic = mean,standard_statistic = 0.5,B=20,alpha = 0.1)[[1]]
}
bb[3]<-mean(bb20)
for(i in 1:10){
bb50[i]<-BOOT_COVER_PROB_BIN(size0 = 1000,p=0.5,statistic = mean,standard_statistic = 0.5,B=50,alpha = 0.1)[[1]]
}
bb[4]<-mean(bb50)
for(i in 1:10){
bb200[i]<-BOOT_COVER_PROB_BIN(size0 = 1000,p=0.5,statistic = mean,standard_statistic = 0.5,B=200,alpha = 0.1)[[1]]
}
bb[5]<-mean(bb200)
for(i in 1:5){
bb1000[i]<-BOOT_COVER_PROB_BIN(size0 = 1000,p=0.5,statistic = mean,standard_statistic = 0.5,B=1000,alpha = 0.1)[[1]]
}
bb[6]<-mean(bb1000)
ind<-c(5,10,20,50,200,1000)
plot(ind,bb,type = "b",ylim=c(0.6,1),xlab="bootstrap次数B",ylab="覆盖率",main="standard normal区间覆盖率与bootstrap次数B的关系图(二项分布)")
#standard normal CI 覆盖率与B的关系(正态分布)
bb<-c()
bb5<-c();bb10<-c();bb20<-c();bb50<-c();bb200<-c();bb1000<-c();bb5000<-c()
for(i in 1:10){
bb5[i]<-BOOT_COVER_PROB_NORM(size0 = 1000,mean0=5,sd0=5,statistic = mean,standard_statistic = 5,B=5,alpha = 0.1)[[1]]
}
bb[1]<-mean(bb5)
for(i in 1:10){
bb10[i]<-BOOT_COVER_PROB_NORM(size0 = 1000,mean0=5,sd0=5,statistic = mean,standard_statistic = 5,B=10,alpha = 0.1)[[1]]
}
bb[2]<-mean(bb10)
for(i in 1:10){
bb20[i]<-BOOT_COVER_PROB_NORM(size0 = 1000,mean0=5,sd0=5,statistic = mean,standard_statistic = 5,B=20,alpha = 0.1)[[1]]
}
bb[3]<-mean(bb20)
for(i in 1:10){
bb50[i]<-BOOT_COVER_PROB_NORM(size0 = 1000,mean0=5,sd0=5,statistic = mean,standard_statistic = 5,B=50,alpha = 0.1)[[1]]
}
bb[4]<-mean(bb50)
for(i in 1:10){
bb200[i]<-BOOT_COVER_PROB_NORM(size0 = 1000,mean0=5,sd0=5,statistic = mean,standard_statistic = 5,B=200,alpha = 0.1)[[1]]
}
bb[5]<-mean(bb200)
for(i in 1:5){
bb1000[i]<-BOOT_COVER_PROB_NORM(size0 = 1000,mean0=5,sd0=5,statistic = mean,standard_statistic = 5,B=1000,alpha = 0.1)[[1]]
}
bb[6]<-mean(bb1000)
ind<-c(5,10,20,50,200,1000)
plot(ind,bb,type = "b",ylim=c(0.6,1),xlab="bootstrap次数B",ylab="覆盖率",main="standard normal区间覆盖率与bootstrap次数B的关系图(正态分布)")
#首先定义偏度skewness函数
sample.skewness = function (sample) {
mu = mean(sample)
n = length(sample)
numerator<-sqrt(n)*sum((sample-mu)^3)
denominator<-(n-1)^(3/2)*sd(sample)^3
#num = 1/n * sum(sapply(sample, function (x) (x - mu)^3))
#denom = sd(sample)^3
return (numerator/denominator)
}
#计算置信区间standard normal,basic,and the percentile
conf.inters=function(xs,B=500,alpha=0.05){
#进行bootstrap抽样,xs为样本,B为抽样次数
theta.hat = sample.skewness(xs)
theta.hats.b = numeric(B)
for (b in 1:B) {
i = sample(1:n, size=n, replace=TRUE)
xs.b = xs[i]
theta.hats.b[b] = sample.skewness(xs.b)
}
#standard normal bootstrap confidence interval
se.theta.hat<-sd(theta.hats.b)
ci.stan.norm<-c(theta.hat+qnorm(alpha/2)*se.theta.hat,theta.hat-qnorm(alpha/2)*se.theta.hat)
#the basic bootstrap confidence interval
ci.basic <-c(2*theta.hat -quantile(theta.hats.b,1-alpha/2),2*theta.hat -quantile(theta.hats.b,alpha/2))
# the percentile confidence interval
ci.percentile = c(quantile(theta.hats.b,alpha/2),quantile(theta.hats.b,1-alpha/2))
#合并生成一个data.frame
conf.inter<-rbind(ci.stan.norm,ci.basic,ci.percentile)
conf.inter<-data.frame(conf.inter)
colnames(conf.inter)<-as.character(c(alpha/2,1-alpha/2))
#加入左偏和右偏
conf.inter$left.miss<-0
conf.inter$right.miss<-0
R<-1000
temp<-numeric(R)
for (r in 1:R) {
i = sample(1:n, n, replace = TRUE)
temp[r] = sample.skewness(xs[i])
for (k in 1:nrow(conf.inter)) {
lower = conf.inter[k,1]
upper = conf.inter[k,2]
conf.inter[k,'left.miss']<- length(temp[temp<lower])/R
conf.inter[k,'right.miss']<-length(temp[temp>upper])/R
}
}
return(conf.inter)
}
#覆盖概率 正态总体 偏度
cov.prob.norm<-function(mu=0,sigma=1,n=200,rep=300,plot=TRUE){
#mu是均值,sigma是方差,n为样本量,rep为bootstrap次数,plot为是否添加置信区间覆盖期望偏度图
stan.norm.cp<-numeric(rep)
basic.cp<-numeric(rep)
perc.cp<-numeric(rep)
low3<-numeric(rep)
high3<-numeric(rep)
#正态总体的偏度为0
skewness<-0
for(i in 1:rep){
xx<-rnorm(n,mu,sigma)
temp<-conf.inters(xx)[1:3,1:2]
low1<-temp[1,1]
high1<-temp[1,2]
stan.norm.cp[i]<-(low1<skewness)&(high1>skewness)
low2<-temp[2,1]
high2<-temp[1,2]
basic.cp[i]<-(low2<skewness)&(high2>skewness)
low3[i]<-temp[3,1]
high3[i]<-temp[3,2]
perc.cp[i]<-(low3[i]<skewness)&(high3[i]>skewness)
}
#对percentile画图
if(plot==TRUE){
max_lim<-max(high3-low3)
plot(c(1:rep),y=c(1:rep),type="n",ylim=c(skewness-max_lim,skewness+max_lim),xlab = "区间序号",ylab ="置信区间",main="percentile置信区间覆盖图")
abline(h=skewness,lwd=1,col=3)
for(i in 1:rep){
segments(i,low3[i],i,high3[i],col=2-perc.cp[i],lwd=2-perc.cp[i])
}}
#计算覆盖概率
mean1<-mean(stan.norm.cp)
mean2<-mean(basic.cp)
mean3<-mean(perc.cp)
cove.prob<-c(mean1,mean2,mean3)
names(cove.prob)<-c("stan.norm.cp","basic.cp","perc.cp")
return(cove.prob)
}
#卡方分布的覆盖概率,df是自由度,n是样本量,rep是抽样次数
cov.prob.chisq<-function(df=5,n=100,rep=200,plot=TRUE){
stan.norm.cp<-numeric(rep)
basic.cp<-numeric(rep)
perc.cp<-numeric(rep)
low3<-numeric(rep)
high3<-numeric(rep)
#卡方分布的偏度
skewness<-sqrt(8/df)
for(i in 1:rep){
xx<-rchisq(n,df)
temp<-conf.inters(xx)[1:3,1:2]
low1<-temp[1,1]
high1<-temp[1,2]
stan.norm.cp[i]<-(low1<skewness)&(high1>skewness)
low2<-temp[2,1]
high2<-temp[2,2]
basic.cp[i]<-(low2<skewness)&(high2>skewness)
low3[i]<-temp[3,1]
high3[i]<-temp[3,2]
perc.cp[i]<-(low3[i]<skewness)&(high3[i]>skewness)
}
#对percentile画图
if(plot==TRUE){
max_lim<-max(high3-low3)
plot(c(1:rep),y=c(1:rep),type="n",ylim=c(skewness-max_lim,skewness+max_lim),xlab = "区间序号",ylab ="置信区间",main="percentile置信区间覆盖图")
abline(h=skewness,lwd=1,col=3)
for(i in 1:rep){
segments(i,low3[i],i,high3[i],col=2-perc.cp[i],lwd=2-perc.cp[i])
}
}
mean1<-mean(stan.norm.cp)
mean2<-mean(basic.cp)
mean3<-mean(perc.cp)
cove.prob<-c(mean1,mean2,mean3)
names(cove.prob)<-c("stan.norm.cp","basic.cp","perc.cp")
return(cove.prob)
}
#二项分布
#n是样本量,size是重复了多少次0-1实验,prob是二项分布的概率,rep是bootstrap次数
cov.prob.bin<-function(n=100,size=100,prob=0.5,rep=200,plot=TRUE){
stan.norm.cp<-numeric(rep)
basic.cp<-numeric(rep)
perc.cp<-numeric(rep)
low3<-numeric(rep)
high3<-numeric(rep)
#卡方分布的偏度
skewness<-(1-2*prob)/sqrt(size*prob*(1-prob))
for(i in 1:rep){
xx<-rbinom(n,size,prob)
temp<-conf.inters(xx)[1:3,1:2]
low1<-temp[1,1]
high1<-temp[1,2]
stan.norm.cp[i]<-(low1<skewness)&(high1>skewness)
low2<-temp[2,1]
high2<-temp[2,2]
basic.cp[i]<-(low2<skewness)&(high2>skewness)
low3[i]<-temp[3,1]
high3[i]<-temp[3,2]
perc.cp[i]<-(low3[i]<skewness)&(high3[i]>skewness)
}
#对percentile画图
if(plot==TRUE){
max_lim<-max(high3-low3)
plot(c(1:rep),y=c(1:rep),type="n",ylim=c(skewness-max_lim,skewness+max_lim),xlab = "区间序号",ylab ="置信区间",main="percentile置信区间覆盖图")
abline(h=skewness,lwd=1,col=3)
for(i in 1:rep){
segments(i,low3[i],i,high3[i],col=2-perc.cp[i],lwd=2-perc.cp[i])
}
}
mean1<-mean(stan.norm.cp)
mean2<-mean(basic.cp)
mean3<-mean(perc.cp)
cove.prob<-c(mean1,mean2,mean3)
names(cove.prob)<-c("stan.norm.cp","basic.cp","perc.cp")
return(cove.prob)
}
首先比较正态分布的相关性质 #1 验证正态样本量和percentile覆盖概率关系
plot(c(100:2000),c(100:2000),type="n",ylim=c(0.6,1),main="正态样本量和percentile覆盖概率关系",xlab="样本量",ylab="prob")
abline(h=0.95,col=3)
l<-length(seq(100,2000,50))
for(i in seq(100,2000,50)){
points(i,cov.prob.norm(n=i,plot=FALSE)[3],cex=0.5)}
#2 验证正态bootstrap次数和percentile覆盖概率关系
plot(c(100:1000),c(100:1000),type="n",ylim=c(0.6,1),main="正态bootstrap次数和percentile覆盖概率关系",xlab="样本量",ylab="prob")
l<-length(seq(100,5000,50))
for(i in seq(100,1000,50)){
points(i,cov.prob.norm(rep=i,plot=FALSE)[3],cex=0.5,col=c(1:5))}
由图可知,正态样本量和正态bootstrap次数对percentile覆盖概率没有影响,当自变量变化时,覆盖概率都在理论值0.95附近波动,说明对正态分布而言,样本量与bootstrap次数与覆盖概率无关
第二步来验证卡方分布的相关性质
由正态分布的性质已知,样本量与bootstrap次数与覆盖概率无关,而对于卡方分布来说是否满足相同的性质?以下为编程验证
#3 验证卡方5分布样本量和percentile覆盖概率关系
plot(c(100:2000),c(100:2000),type="n",ylim=c(0.6,1),main="样本量和percentile覆盖概率关系")
l<-length(seq(100,2000,50))
for(i in seq(100,2000,50)){
points(i,cov.prob.chisq(df=5,n=i,plot=FALSE)[3],cex=0.5)}
#4 验证卡方100bootstrap次数和percentile覆盖概率关系
plot(c(10:1000),c(10:1000),type="n",ylim=c(0.6,1),main="卡方bootstrap次数和percentile覆盖概率关系",xlab="bootstrap次数",ylab="prob")
l<-length(seq(10,1000,50))
for(i in seq(10,1000,50)){
points(i,cov.prob.chisq(df=100,rep=i,plot=FALSE)[3],cex=0.5)}
由图可知,当卡方分布的自由度固定时,样本量和bootstrap次数都与覆盖概率无关,但是图片显示出卡方5的覆盖概率明显低于理论值0.95,而卡方100的覆盖概率则在0.95附近,于是猜想卡方分布的自由度对覆盖概率有影响
#5 验证卡方自由度与percentile覆盖概率关系
plot(c(1:300),c(1:300),type="n",ylim=c(0.6,1),main="自由度与percentile覆盖概率关系")
for(i in seq(1,300,5)){
points(i,cov.prob.chisq(df=i,plot="FALSE")[3],cex=0.5)}
当自由度从1开始递增到300时,在(0,100)可以看出覆盖概率有明显的递增趋势,当自由度大于100时,覆盖概率在0.95附近波动。由理论上的知识,卡方分布的偏度为sqrt(8/n),n为自由度,当自由度增加时,偏度趋于0,推测自由度的变化带来偏度的变化,并且偏度越大,覆盖概率越小,为了防止只抽取一种样本带来特殊情况,在选取二项分布验证推测
#6 验证二项prob和percentile覆盖概率关系
plot(seq(0,1,0.05),seq(0,1,0.05),type="n",ylim=c(0.6,1),main="二项prob和percentile覆盖概率关系",xlab="二项概率",ylab="prob")
abline(h=0.95,col=3)
for(i in seq(0.005,0.05,0.005)){
points(i,cov.prob.bin(prob=i,plot=FALSE)[3],cex=0.5)}
for(i in seq(0.95,0.995,0.005)){
points(i,cov.prob.bin(prob=i,plot=FALSE)[3],cex=0.5)}
for(i in seq(0.05,0.95,0.1)){
points(i,cov.prob.bin(prob=i,plot=FALSE)[3],cex=0.5)}
二项分布偏度 (1-2prob)/sqrt(sizeprob*(1-prob)),其中size为0-1实验次数,prob为成功率
出于对绘图时间的考虑,只在重点区间(即两侧)密集取点,在中间稀疏取点。
在prob从0.005到0.995变化时,偏度由正偏度到零到负偏度,偏度的绝对值由大到零再增大,由图表可知,在坐标轴的左侧,出现了明显的递增趋势,而在坐标轴的右侧,出现了明显的递减趋势,在坐标轴中间,覆盖概率在0.95附近波动。因此,总体偏度越大,percentile方法的覆盖概率越小,即偏度越大,越“测不准”
1.样本量,bootstrap次数与覆盖概率无关
2.决定总体偏度的参数与覆盖概率有关,且总体偏度越大覆盖概率越小,总体偏度接近0时,覆盖概率在理论值0.95附近波动
思考:
1.在编写boot函数的过程中,一度出现了覆盖概率一直为1的情况,后来发现是因为每次的boot抽样没有放入循环中,这是以前所忽视的
2.在验证各种性质时,出现的结果并不是一直与猜想的结果相同,有的甚至是意料之外的,但是这些意料之外的结果还是符合其中的数学原理的,体现了科学的严谨性
3.通过用bootstrap方法来计算各种分布的覆盖概率令我们更加深入的了解了bootstrap方法,同时也证明了覆盖概率的一些基本性质,收获颇丰
经过对均值与偏度的有关覆盖率分析,让我们对bootstrap方法与bootstrap置信区间有了较深的理解,并且深刻剖析了bootstrap方法中诸如boot次数、区间类型、样本等参数与覆盖率的内在联系与相关性。更在研究过程中详细讨论了二项分布的离散与偏斜性在覆盖率上的展现,有关覆盖率振荡等效果的体现与阐述。本文对覆盖率的参数相关性与内蕴性质的研究,不仅对相关方面的研究,更对现实中置信区间覆盖率的应用有所启示与参考。
人员分工详见目录
参考文献:Statistical Computing with R;bootstrap方法下的poisson分布置信区间估计; 二项分布参数置信区间的比较