bootstrap怎么强制不换行_Bootstrap方法的软件实现

注：本文是叶楚韵对Bootstrap方法在R、MATLAB、Python以及SAS中实现的介绍

通过之前关于Bootstrap的介绍 (重抽样方法)，相信大家已经对Bootstrap的原理和用途有了一定的了解。本文将为大家介绍在R、MATLAB、Python以及SAS中实现Bootstrap的现成语句。注意Bootstrap的原理是重抽样，因此在所有语言中，均可以采用循环语句实现Bootstrap。这种做法的优点是形式灵活，可以实现更多不同方法的Bootstrap，缺点是可能对初学者不友好，编写代码时容易出纰漏。因此，下面的介绍将以实现Bootstrap的现成命令为主。

1 Bootstrap的R实现

目前比较常用于实现Bootstrap的R现成语句主要来自boot、bootstrap以及simpleboot包，均为第三方包，需运行install.packages('')下载并安装，以及每次使用前需运行library()加载。下面关于bootstrap的R实现的所有实例均采用R自带数据集cd4。

1.1 boot package

利用boot包中的函数boot()与boot.ci()函数，以及boot对象，我们可以实现一些简单的Bootstrap过程。下面将展示几个简单的例子并辅以说明：

• 单样本、成对两样本情况

mean.boot = function(data,indices){  d = data[indices,]  res = mean(d[,1])  return(res)}result1 = boot(data=cd4,statistic=mean.boot,R=500)boot.ci(result1,conf=0.95,type = "all")

boot()函数几个重要参数分别为data，statistic以及R，其中，data参数传递原数据，R传递Bootstrap重复次数B。statistic传递统计量形式，在上例，即为数据第1列的均值。因此statistic参数是一个自定义函数，在上例即为mean.boot()。在一般情况下，boot函数要求这个自定义函数必须含输入参数data与indices。通过上例可以看到，参数indices就是在boot函数重抽样过程中抽取数据的索引。

boot函数的输出参数为boot对象。函数boot.ci()可以直接作用于这一类型变量并输出区间估计。调整参数conf可以获得不同置信度的区间估计，调整参数type可以获得不同类型的区间估计，目前支持输出5种区间估计，感兴趣的读者可以查阅boot包说明文件。

• 配对Bootstrap回归

paired.regcoe.boot = function(formula,data,indices){  d = data[indices,]  fit = lm(formula,data=d)  return(fit$coef)}result2 = boot(data=cd4,statistic=paired.regcoe.boot,R=500,formula=oneyear~baseline)boot.ci(result2,conf=0.95,type = "all",index = 1)boot.ci(result2,conf=0.95,type = "all",index = 2)

注意到通过indices参数可以实现按抽取到的索引提取数据，因此，通过boot()，我们也可以实现配对Bootstrap方法估计回归系数。而在使用boot.ci()函数求置信区间时，可利用index参数区分多维估计量中的不同分量。

1.2 bootstrap package

bootstrap包中的函数bootstrap()与boott()函数是我们实现简单Bootstrap的主要工具，下面展示几个简单的例子：

• 单样本情况、成对两样本情况、配对Bootstrap回归

mean.bootstrap = function(x) {    mean(x)}result3 = bootstrap(x=cd4[,1],nboot=500,theta=mean.bootstrap)boott(x=cd4[,1],theta=mean.bootstrap,nboott=500,perc=c(0.01,0.025,0.1,0.9,0.975,0.99))boott(x=cd4[,1],theta=mean,nboott=500,perc=c(0.01,0.025,0.1,0.9,0.975,0.99))

通过bootstrap()函数，我们可以对数据进行简单Bootstrap，通过boott()函数，我们可以计算Bootstrap-t置信区间。其中，参数x传递数据(向量类型),theta传递统计量形式，nboott传递Bootstrap重复次数B，perc传递对应的概率。bootstrap()函数输出包含此次Bootstrap信息的列表，boott()函数输出置信区间端点。

与boot包不同，此处的参数theta虽然要求是函数形式，但在形式上不要求含indices输入参数，读者可以对比自定义函数mean.boot与mean.bootstrap的不同。同时，theta也支持R中原有的函数输入，读者可以对比上例第五行与第六行。

注意到参数x仅支持输入向量，那么如果我们讨论的数据类型比向量更复杂，应该怎么处理呢？我们看看计算样本相关系数的例子：

cor.bootstrap = function(x,data){  cor(data[x,1],data[x,2])}boott(x=1:nrow(dat),theta=theta,data=cd4,nboott=500,perc=c(0.01,0.025,0.1,0.9,0.975,0.99))

注意到此处统计量函数cor.bootstrap()有两个输入参数，其中x代表抽取数据的索引，是一个向量，满足boott()函数的要求。以相似的思路，我们利用boott()与bootstrap()函数，可以完成配对Bootstrap回归系数估计。

1.3 simpleboot package

simpleboot包包含了一些在boot包基础上按功能集成的Bootstrap代码，在实现简单Bootstrap这方面比直接使用boot包与bootstrap包更方便。下面将展示几个例子，注意在使用simpleboot包前除了需要加载simpleboot包本身，还需要加载boot包。

• 单样本情况

library(boot)mav.simpleboot = function(x) {    c(mean(x),var(x))}result4 = one.boot(data=cd4[,1],FUN=mav.simpleboot,R=500)boot.ci(result4,index=1)boot.ci(result4,index=2)

通过one.boot()函数，可以实现单样本Bootstrap，其中输入参数FUN传递统计量形式，注意这里代表统计量形式的自定义函数不需要含输入参数indices，且可以使用R自带函数，这一设定反而与bootstrap包更接近。one.boot()函数输出boot对象，我们同样可以用boot.ci()作用于boot对象求Bootstrap不同类型的置信区间。

• 成对两样本情况

simpleboot包还提供了实现成对两样本Bootstrap的专门函数pairs_boot()。

cor.simpleboot = function(x,y){  cor(x,y)}result5 = pairs_boot(x=cd4[,1],y=cd4[,2],FUN=cor.simpleboot,R=500)boot.ci(result5)

• 不成对两样本情况

利用boot包与bootstrap包都不能比较直接地实现不成对两样本的Bootstrap推断。simpleboot包提供了two.boot()函数，实现对不配对两样本的统计量之差的Bootstrap推断。

quantile.simpleboot = function(x,probs){  quantile(x,probs)}result6 = two.boot(sample1=cd4[,1],sample2=cd4[,2],FUN=quantile.simpleboot,R=500,probs=0.75)boot.ci(result6)

其中参数sample1与sample2分别代表两组样本。

• 配对Bootstrap回归、基于残差的Bootstrap回归

simpleboot包中的lm.boot()函数可以作用于lm对象，实现回归配对Bootstrap或回归基于残差Bootstrap。rows参数为TRUE则为配对Bootstrap，rows参数为FALSE则为基于残差Bootstrap。通过simpleboot包中的perc函数，可以求出对应的百分位置信区间。

# paired Bootstraplmodel1 = lm(formula=oneyear~baseline,data=cd4)result7 = lm.boot(lm.object=lmodel1,R=500,rows=T)summary(result7)perc(result7,p=c(0.025,0.05,0.95,0.975))# residual Bootstrapresult8 = lm.boot(lm.object=lmodel1,R=500,rows=F)summary(result8)perc(result8,p=c(0.025,0.05,0.95,0.975))

1.4 Wild-Bootstrap

注意到，在前面关于回归Bootstrap的介绍中，除了配对Bootstrap与基于残差Bootstrap，还有Wild Bootstrap。fANCOVA包提供了wild.boot()函数实现Wild Bootstrap，下面是一个简单的例子。

lmodel2 = lm(oneyear~baseline,data=cd4)res.boot = wild.boot(lmodel2$residuals,nboot=500)y.boot = matrix(rep(lmodel2$fitted.values,time=500),ncol=500) + res.bootres = matrix(rep(0,2*500),nrow=2)for (i in 1:500) {  res[,i] = lm(y.boot[,i]~cd4[,1])$coef}#loopingres=apply(X=y.boot,MARGIN=2,FUN=function(y){  lm(y~cd4[,1])$coef})#apply function

注意wild.boot()函数作用于残差序列，仅实现了对残差处理与重抽样，剩余的其他步骤可以通过循环完成，或者为了提升运算速度，可以通过apply类函数完成。

1.5 循环与向量化运算

根据Bootstrap的原理，使用循环语句可以实现Bootstrap。R的向量化运算函数(apply类函数)可以更高效地完成循环语句的工作。下面给出一个向量化运算函数replicate()实现成对Bootstrap回归的例子。

beta = function(formula,data,indices){  d = data[indices,]  fit = lm(formula,data=d)  return(fit$coef)}result9 = replicate(n=500,expr={  n = nrow(dat)  ind = sample(1:n,n,replace=T)  beta(formula=oneyear~baseline,data=cd4,indices=ind)})

函数replicate()可以实现对一个函数的重复操作，利用这一性质，我们可以实现一些简单的Bootstrap。

2 Bootstrap的MATLAB实现

在MATLAB中，实现Bootstrap的现成语句主要是bootstrp()与bootci()，它们分别输出重抽样统计量序列与Bootstrap区间估计。简单看一下这两个语句的常用语法：

bootstat = bootstrp(nboot,bootfun,d)bootstat = bootstrp(nboot,bootfun,d1,...,dN)ci = bootci(nboot,{bootfun,...},'alpha',alpha)ci = bootci(nboot,{bootfun,...},...,'type',type)

其中，参数nboot传递Bootstrap重复次数，参数bootfun传递统计量形式，要求是函数句柄。参数d,d1,...,dN传递数据，参数alpha传递置信度，参数type传递置信区间类型,percentile表示使用百分位置信区间，student表示使用Bootstrap T置信区间。

通过这两个语句，我们可以实现对单样本Bootstrap、成对两样本的Bootstrap以及基于残差的Bootstrap回归。下面将通过几个MathWorks帮助中心提供的例子进行说明。

• 单样本情况

y = exprnd(5,100,1);m = bootstrp(500,@mean,y);bootci(500,{@mean,y},'type','percentile','alpha',0.95)bootci(500,{@mean,y},'type','student','alpha',0.95)

注意到这里将MATLAB函数mean()定义成了函数句柄@mean作为统计量形式，类似的操作我们还会在下面看到。

当考虑多个统计量时，我们可以进行类似的操作：

y = exprnd(5,100,1);stats = bootstrp(500,@(x)[mean(x) std(x)],y);bootci(500,{@(x)[mean(x) std(x)],y},'type','percentile','alpha',0.95)bootci(500,{@(x)[mean(x) std(x)],y},'type','student','alpha',0.95)

• 成对两样本

此时用到的语法为上面列出bootstrp()函数语法的第二行，注意到这种语法不仅可以处理成对两样本，还可以处理成对多样本。

load lawdatasize(lsat)size(gpa)[bootstat,bootsam] = bootstrp(1000,@corr,lsat,gpa);bootci(500,{@corr,lsat,gpa},'type','percentile','alpha',0.95)bootci(500,{@corr,lsat,gpa},'type','student','alpha',0.95)

• 基于残差的Bootstrap回归

按下例的思路，可以简洁地实现基于残差的Bootstrap回归。

load haldx = [ones(size(heat)),ingredients];y = heat;b = regress(y,x);yfit = x*b;resid = y - yfit;res_boot = bootstrp(500,@(bootr) regress(yfit+bootr,x),resid);bootci(500,{@(bootr) regress(yfit+bootr,x),resid},'type','percentile','alpha',0.95)bootci(500,{@(bootr) regress(yfit+bootr,x),resid},'type','student','alpha',0.95)

3 Bootstrap的Python实现

Python中有很多第三方模块可以帮助我们完成不同的任务，包括Bootstrap，其中比较常用的是bootstrapped模块。注意在使用第三方模块前需要先下载和安装。在这里给出官方介绍中的使用范例，其他细节可以参见官方介绍^[1]

import numpy as npimport bootstrapped.bootstrap as bsimport bootstrapped.stats_functions as bs_statsmean = 100stdev = 10population = np.random.normal(loc=mean, scale=stdev, size=50000)# take 1k 'samples' from the larger populationsamples = population[:1000]print(bs.bootstrap(samples, stat_func=bs_stats.mean))>> 100.08  (99.46, 100.69)print(bs.bootstrap(samples, stat_func=bs_stats.std))>> 9.49  (9.92, 10.36)

可以看到bootstrapped模块中的bootstrap函数是基于模块中的stats_functions中的函数进行自助法，并返回点估计与区间估计。

4 Bootstrap的SAS实现

在SAS中，我们也可以利用循环的思路实现Bootstrap。因此，简单地说，我们可以使用DO语句(循环)、SURVEYSELECT过程(重抽样)以及SAS/IML模块(交互式矩阵语言)实现Bootstrap。关于如何利用SAS实现更多功能，读者可以参考SAS官网^[2]，其中的support模块与blog模块都对不同功能的实现有着详尽的解释与研究。下面我们简要介绍其他实现Bootstrap的方法，并给出网络上的相关介绍。

• bootstrap Statement

在SAS在一些PROC步中提供可用的bootstrap语句，以TTEST PROC为例——

proc ttest data=Sample;   class Type;   var MPG_City;   bootstrap / seed=123 nsamples=10000 bootci=percentile;  /* or BOOTCI=BC */run;

bootstrap语句的相关选项可以参见NLIN PROC下关于bootstrap Statement的说明^[3]

• 宏

SAS MACRO模块提供%boot()与%bootci()函数实现Bootstrap，具体可以参见SAS SUPPORT中的例子^[4]

如果觉得本文不错，请点赞关注！

参考资料

Python: bootstrapped 0.0.2: https://pypi.org/project/bootstrapped/

SAS homepage: https://www.sas.com/zh_cn/home.html

SAS: The NLIN Procedure: https://support.sas.com/documentation/cdl/en/statug/68162/HTML/default/viewer.htm#statug_nlin_syntax02.htm

SAS: Sample 24982: Jackknife and Bootstrap Analyses: https://support.sas.com/kb/24/982.html