R语言自己编程实现bootstrap函数和调用函数包实现

dataset=read.table(e:/CH01TA01.txt,T)#使用数据样本head(dataset)lm=lm(Y~X,dataset)summary(lm)boooootstrap-function(dataset,m,n)#自编函数实现bootstrap功能{ff=as.formula(paste(names(dataset)[m],~.,sep=))row=nrow(dataset);col=ncol(dataset)c=matrix(0,n,col)a=lm(ff,dataset)dataset_new=datasetfor(iin1:n){num=sample(1:row,row,replace=T)dataset_new[,m]=dataset[,m]+(a$resid)[num]a1=lm(ff,dataset_new)c[i,]=as.vector(a1$coef)}c=as.data.frame(c)b=list()for(iin1:col)b[[i]]=c[,i]names(b)=c(Intercept,colnames(dataset)[-m])return(b)}(b=boooootstrap(dataset,2,500))hist(b[[i]])#i为第i个参数###################函数boooootstrap的输入：1.dataset为数据集#2.m为数据集中因变量所在列#3.n为所需boooootstrap的次数;#输出：列表b，每一层为截距以及各自变量在boooootstrapn次所得的值CI-function(dataset,m,n,x){a=1-xb=boooootstrap(dataset,m,n)ff=as.formula(paste(names(dataset)[m],~.,sep=))fit=lm(ff,dataset)d=list()for(iin1:length(b)){d1=as.vector(2*fit$coef[i]-quantile(b[[i]],1-a/2))d2=as.vector(2*fit$coef[i]-quantile(b[[i]],a/2))d[[i]]=c(d1,d2)}names(d)=c(Intercept,colnames(dataset[-m]))return(d)}CI(dataset,2,100,0.95)#函数CI的输入：1.dataset为数据集#2.m为数据集中因变量所在列#3.n为所需bootstrap的次数#4.x为置信度(0~1)#输出：列表d，每一层为截距以及各自变量的置信区间#运行函数CI前应先加载函数boooootstrap#############################引用bootstrap包计算install.packages(bootstrap)library(bootstrap)theta.b-function(x,dataset,i){ff=as.formula(paste(names(dataset)[i],~.,sep=))out=lm(ff,dataset,subset=x)coe=out$coefficientsreturn(coe)}set.seed(333)n=nrow(dataset)coe=bootstrap(1:n,500,theta.b,dataset,i)#1:n为抽样空间，可s为样本空间的子集；#500为bootstrap次数#theta为构造的函数；#dataset为数据集对象#i为因变量所在列apply(coe$thetastar,1,mean)#对拟合系数按行求均值hist(coe$thetastar[n,])#b_(n-1)的频率分布图ff=as.formula(paste(names(dataset)[i],~.,sep=))#i需赋值为因变量所在列值fit=lm(ff,dataset)L=U=0for(iin1:length(dataset)){L[i]=2*fit$coef[i]-quantile(coe$thetastar[i,],1-0.05/2)U[i]=2*fit$coef[i]-quantile(coe$thetastar[i,],0.05/2)}(CI=data.frame(L,U,row.names=c(Intercept,colnames(dataset)[-i])))#b0至b_(n-1)的置信区间################################fit=fitted(dataset)e=residuals(dataset)x=model.matrix(dataset)################################alpha-function(dataset,para,i){ff=as.formula(paste(names(dataset)[i],~.,sep=))d=dataset[para,]out=lm(ff,d)return(coef(out)[2])}(result=boot(data=dataset,statistic=alpha,R=500,i=2))boot.ci(result)plot(result)