dividedata.R-20170830

#data=cbind(ypart,xcovpart2)
#sample(2, nrow(data), replace=T, prob=c(0.7,0.3))
library(caret) 

data=as.matrix(oadata)
colnames(data)=c('y','x1','x2','x17',paste("x",3:16,sep=""))

c=5000 #验证c次
set.seed(5000)
r1=matrix(0,c,1)
r2=matrix(0,c,1)
r1_c=matrix(0,c,1)
r2_c=matrix(0,c,1)
p1=matrix(0,12,c)
p1_c=matrix(0,12,c)

for (k in 1:c) {
index <-createDataPartition(data[,1], time=1, p=0.8, list=F)
# if (table(test[,5])>0 && table(test[,6])>0 && table(test[,7])>0 && table(test[,8])>0 &&
#     table(test[,9])>0 && table(test[,10])>0 && table(test[,11])>0 && table(test[,12])>0 &&
#     table(test[,13])>0 && table(test[,14])>0 && table(test[,15])>0 && table(test[,16])>0 &&
#     table(test[,17])>0 && table(test[,18])>0 && table(train[,5])>0 && table(train[,6])>0 && 
#    table(train[,7])>0 && table(train[,8])>0 && table(train[,9])>0 && table(train[,10])>0 && 
#    table(train[,11])>0 && table(train[,12])>0 && table(train[,13])>0 && table(train[,14])>0 && 
#    table(train[,15])>0 && table(train[,16])>0 && table(train[,17])>0 && table(train[,18])>0)
#   { 
#     print(k)
  train=data[index, ]
  test=data[-index, ]
#print(nrow(train))
#print(nrow(test))
colnames(train)=c('y','x1','x2','x17',paste("x",3:16,sep=""))
colnames(test)=c('y','x1','x2','x17',paste("x",3:16,sep=""))
#assign(paste("train",k,sep=""),data[index, ])
#assign(paste("test",k,sep=""),data[-index, ])

a.lm = 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)) #train集构造anova模型,删去x6，x16
#assign(paste("a.lm",k,sep=""),a.lm) #记为ak
#summary(a.lm)
#assign(paste("beta",k,sep=""),data.matrix(coef(a.lm))) #系数betak
ytest=predict(aov(a.lm),data.frame(test[,2:18])) 
#assign(paste("ytest",k,sep=""),ytest) #test集预测结果ytestk
resi1=abs(ytest-test[,1])/test[,1]
r1[k]=mean(resi1) #误差r1(k)
#print(mean(resi1))
p1[,k]=Anova(a.lm,singular.ok = TRUE,type="III")$Pr #p-value
}

# #交叉验证
# a_c.lm = 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(test))
# #assign(paste("a_c.lm",k,sep=""),a_c.lm)
# #summary(a_c.lm)
# #assign(paste("beta2_",k,sep=""),data.matrix(coef(a_c.lm))) #系数beta2_k
# ytrain=predict(a_c.lm,data.frame(train[,2:18]))
# #assign(paste("ytrain",k,sep=""),ytrain)
# resi1_c=abs(ytrain-train[,1])/train[,1]
# r1_c[k]=mean(resi1_c) #误差
# #print(mean(r1_c))
# p1_c[,k]=Anova(a_c.lm,singular.ok = TRUE,type="III")$Pr #p-value
# }

m<-svm(train[,2:18],train[,1])
#assign(paste("m",k,sep=""),m)
#summary(m)
ytest2=predict(m,test[,2:18])
#assign(paste("ytest2_",k,sep=""),ytest2)
resi2=abs(ytest2-test[,1])/test[,1]
r2[k]=mean(resi2) #误差
#print(mean(resi2))

# #交叉验证
# m_c<-svm(test[,2:18],test[,1])
# #assign(paste("m_c",k,sep=""),m_c)
# #summary(m_c)
# ytrain2=predict(m_c,train[,2:18])
# #assign(paste("ytrain2_",k,sep=""),ytrain2)
# resi2_c=abs(ytrain2-train[,1])/train[,1]
# r2_c[k]=mean(resi2_c) #误差
# #print(mean(resi2_c))



mean(r1)
#mean(r1_c)
mean(r2)
#mean(r2_c)
#rowMeans(p1)
#rowMeans(p1_c)

# library(ISLR)
# library(caret)
# data(Wage)
# Wage<- subset(Wage,select=-c(logwage))
# #Split the data
# inBuild<- createDataPartition(y=Wage$wage,p=0.7,list=FALSE)
# validation<- Wage[inBuild,]
# buildData<- Wage[-inBuild,]
# inTrain<- createDataPartition(y=buildData$wage,p=0.7,list=FALSE)
# training<- buildData[inTrain,]
# testing<- buildData[-inTrain,]
# mod1<- train(wage~.,method="glm",data=training)
# mod2<- train(wage~.,method="rf",data=training,trcontrol=trainControl(method="CV"),number=3)
# pred1<- predict(mod1,testing)
# pred2<- predict(mod2,testing)
# qplot(pred1,pred2,col=wage,data=testing)
#随机森林
# mod2<- train(y~.,method="rf",data=train,trcontrol=trainControl(method="CV"),number=3)
# pred2<- predict(mod2,test)
# print(mean(abs(predict(mod2,data.frame(test[,2:18]))-test[,1])/test[,1]))