bootstrap.R-20170901

# #bootstrap验证
# #返回回归系数向量的自编函数
# bs <- function(formula, data, indices){
#  d <- data[indices,]
#  fit <- lm(formula,data=d)
#  return(coef(fit))
#  }
# #自助抽样
# set.seed(1234)
# a = lm(y~0+x1+x2+x17+as.factor(x3)+as.factor(x4)+as.factor(x5)+as.factor(x6)+as.factor(x7)+as.factor(x8)
#           +as.factor(x9)+as.factor(x10)+as.factor(x11)+as.factor(x12)+as.factor(x13)+as.factor(x14)
#           +as.factor(x15)+as.factor(x16),data=data.frame(train)) #anova模型
# bootobject <- boot(data=, statistic=, R=, ...)
# boot.ci(bootobject, conf=, type= )
# results <- bootstrap(data=mtcars, statistic = bs, R=1000, formula=a)
# #输出 boot 对象
# print(results)
# plot(results,index=2)
# #计算置信区间
# boot.ci(results,type = "bca", index=2)
# boot.ci(results,type = "bca", index=3)
# # index 参数用于指明需要分析的重抽样所得样本的列

data=as.matrix(oadata)
colnames(data)=c('y','x1','x2','x17',paste("x",3:16,sep=""))
#置换检验
a_all = lm(y~0+x1+x2+x17+as.factor(x3)+as.factor(x4)+as.factor(x5)+as.factor(x6)+as.factor(x7)+as.factor(x8)
          +as.factor(x9)+as.factor(x10)+as.factor(x11)+as.factor(x12)+as.factor(x13)+as.factor(x14)
          +as.factor(x15)+as.factor(x16),data=data.frame(data)) #全部样本构造anova模型
#summary(a_all)
#assign(paste("beta",k,sep=""),data.matrix(coef(a_all))) #系数betak
#ytest=predict(aov(a_all),data.frame(test[,2:18])) 
#assign(paste("ytest",k,sep=""),ytest) #test集预测结果ytestk
#resi=abs(ytest-test[,1])/test[,1]
#r=mean(resi1) #误差r1(k)
#print(mean(resi1))
#p=Anova(a_all,singular.ok = TRUE,type="III")$Pr #p-value
f=Anova(a_all,singular.ok = TRUE,type="III")$F #F-value

set.seed(1000)
R=5000 #有放回的抽样R次
pf=matrix(0,R,18)
for (i in 1:R){
  #产生抽样样本
  k=sample(1:length(data[,1]), size=1000, replace=TRUE)
  pdata=data[k,] #构造1000样本
  #colnames(pdata)=c('y','x1','x2','x17',paste("x",3:16,sep=""))
  pa = lm(y~0+x1+x2+x17+as.factor(x3)+as.factor(x4)+as.factor(x5)+as.factor(x6)+as.factor(x7)+as.factor(x8)
             +as.factor(x9)+as.factor(x10)+as.factor(x11)+as.factor(x12)+as.factor(x13)+as.factor(x14)
             +as.factor(x15)+as.factor(x16),data=data.frame(pdata)) 
  pf[i,]=Anova(pa,singular.ok = TRUE,type="III")$F #计算样本的f值
}
pfsort=matrix(0,nrow(pf),ncol(pf))
for (i in 1:ncol(pf)-1) {
  pfsort[,i]=sort(pf[,i])}

pp=matrix(0,1,ncol(pf)-1)
for (i in 1:ncol(pf)-1){
pp[i]=length(which(f[i]>pfsort[,i]))/nrow(pf)} #构造得出的pvalue
colnames(pp)=c('x1','x2','x17',paste("x",3:16,sep=""))
print(pp)