#交叉验证模板
#要确保数据没有缺失值
#确保因子化所有定性变量而且去掉诸如ID之类的多余变量
#先运行下面程序, 把它们放入内存:
####分类分类######################
Fold=function(Z=10,w,D,seed=7777){
n=nrow(w);d=1:n;dd=list()
e=levels(w[,D]);T=length(e)#因变量T类
set.seed(seed)
for(i in 1:T){
d0=d[w[,D]==e];j=length(d0)
ZT=rep(1:Z,ceiling(j/Z))[1:j]
id=cbind(sample(ZT,length(ZT)),d0);dd[]=id}
#上面每个dd[]是随机1:Z及i类的下标集组成的矩阵
mm=list()
for(i in 1:Z){u=NULL;
for(j in 1:T)u=c(u,dd[[j]][dd[[j]][,1]==i,2])
mm[]=u} #mm[]为第i个下标集i=1,...,Z
return(mm)}#输出Z个下标集
######回归回归回归回归###############
CV=function(n,Z=10,seed=888){
z=rep(1:Z,ceiling(n/Z))[1:n]
set.seed(seed);z=sample(z,n)
mm=list();for (i in 1:Z) mm[]=(1:n)[z==i]
return(mm)}#输出Z个下标集;mm[]为第i个下标集i=1到Z
####################
####################
library(MASS);library(rpart.plot);library(adabag)
library(randomForest);library(kernlab)
library(mda);library(kknn);detach(package:ipred)
library(Ecdat)#两个数据所在程序包
################################################################
#第一个数据Bids (分类)
for(i in 8:12)Bids[,i]=factor(Bids[,i])#几个定性变量因子化
w=Bids[,-c(1,4)]#删除两个变量
D=10;ff=(whtknght~.)#因变量在w是第10个, 名字叫whtknght, 这里ff是分类公式
(n=nrow(w))#n=126
Z=10#10折
mm=Fold(Z,w,D)
############
ERROR=matrix(0,Z,10)
for(i in 1:Z)
{m=mm[]
a=mda(ff,w[-m,])#灵活线性判别
ERROR[i,1]=sum(w[m,D]!=predict(a,w[m,]))/length(m)}
for(i in 1:Z)
{m=mm[]
a2=rpart(ff,w[-m,]) #决策树分类
ERROR[i,2]=sum(w[m,D]!=predict(a2,w[m,],type="class"))/length(m)}
set.seed(1010);for(i in 1:Z)
{m=mm[]
a3=boosting(ff,w[-m,]) #boosting分类
ERROR[i,3]=sum(w[m,D]!=predict(a3,w[m,])$class)/length(m)}
set.seed(1010);for(i in 1:Z)
{m=mm[]
a4=bagging(ff,data =w[-m,])
ERROR[i,4]=sum(w[m,D]!=predict(a4,w[m,])$class)/length(m)}
set.seed(1010);for(i in 1:Z)
{m=mm[]
a5=randomForest(ff,data=w[-m,])
ERROR[i,5]=sum(w[m,D]!=predict(a5,w[m,]))/length(m)}
for(i in 1:Z)
{m=mm[]
a6=kknn(ff, train= w[-m,],test=w[m,])#KKNN
ERROR[i,6]=sum(w[m,D]!=a6$fit)/length(m)}
for(i in 1:Z)
{m=m=mm[]
a7=ksvm(ff,w[-m,])#
ERROR[i,7]=sum(w[m,D]!=predict(a7,w[m,]))/length(m)}
for(i in 1:Z)
{m=m=mm[]
a8=fda(ff,w[-m,])#
ERROR[i,8]=sum(w[m,D]!=predict(a8,w[m,]))/length(m)}
for(i in 1:Z)
{m=m=mm[]
a9=lda(ff,w[-m,])#
ERROR[i,9]=sum(w[m,D]!=predict(a9,w[m,])$class)/length(m)}
for(i in 1:Z)
{m=mm[]
a10=glm(ff,w[-m,],family="binomial")#
pp=predict(a10,w[m,],type="response")
bb=levels(w[,D])
ss=rep(bb[1],length(m));ss[pp>.5]=bb[2]
p=sum(w[m,D]!=ss)/length(m)
ERROR[i,10]=p}
ERROR=data.frame(ERROR)
names(ERROR)=c("mda","tree","boost","bagging","random forest",
"kknn","svm","fda","lda","logistic")
(ME=apply(ERROR,2,mean));ERROR
###################################################################
对上面数据单独做Logistic回归
a=glm(ff,w,family="binomial")
pp=predict(a,w,type="response")
bb=levels(w[,D])
ss=rep(bb[1],n);ss[pp>.5]=bb[2]
p=sum(w[,D]!=ss)/nrow(w)
(Error=p)
###################################################################
#第二个数据Clothing (回归)
dim(Clothing)#400 13
(n=nrow(Clothing))#n=400
#统一用w表示数据集:
w=Clothing
D=1;gg=(tsales~.)#因变量在Clothing是第1个, 名字叫tsales, 这里gg是分类公式
library(rpart.plot);library(ipred);library(mboost);
library(randomForest);library(kernlab)
library(nnet)#;library(rminer);detach(package:adabag)
zy=(1:ncol(w))[-D]
gg1=paste(names(w)[D],"~btree(",names(w)[zy[1]],")",sep="")
for(i in (1:ncol(w))[-D][-1])gg1=paste(gg1,"+btree(",names(w),")",sep="")
gg1=as.formula(gg1)
Z=10;mm=CV(nrow(w),Z,1010)
MSE=matrix(0,Z,7)
for(i in 1:Z)
{m=mm[];M=mean((w[m,D]-mean(w[m,D]))^2)
a1=lm(gg,w[-m,])#线性回归
MSE[i,1]=mean((w[m,D]-predict(a1,w[m,]))^2)/M}
for(i in 1:Z)
{m=mm[];M=mean((w[m,D]-mean(w[m,D]))^2)
a2=rpart(gg,w[-m,]) #决策树回归
MSE[i,2]=mean((w[m,D]-predict(a2,w[m,]))^2)/M}
set.seed(1010);for(i in 1:Z)
{m=mm[];M=mean((w[m,D]-mean(w[m,D]))^2)
a3=mboost(gg1,data =w[-m,])#####注意这里用gg1
MSE[i,3]=mean((w[m,D]-predict(a3,w[m,]))^2)/M}
set.seed(1010);for(i in 1:Z)
{m=mm[];M=mean((w[m,D]-mean(w[m,D]))^2)
a4=bagging(gg,data =w[-m,])
MSE[i,4]=mean((w[m,D]-predict(a4,w[m,]))^2)/M}
set.seed(1010);for(i in 1:Z)
{m=mm[];M=mean((w[m,D]-mean(w[m,D]))^2)
a5=randomForest(gg,data=w[-m,])
MSE[i,5]=mean((w[m,D]-predict(a5,w[m,]))^2)/M }
for(i in 1:Z)
{m=mm[];M=mean((w[m,D]-mean(w[m,D]))^2)
a6=nnet(w[-m,-D],w[-m,D],size=5,rang=0.1,decay=5e-4,maxit = 200)
MSE[i,6]=mean((w[m,D]-predict(a6,w[m,]))^2)/M }
for(i in 1:Z)
{m=mm[];M=mean((w[m,D]-mean(w[m,D]))^2)
a7=ksvm(gg,w[-m,],model="svm")
MSE[i,7]=mean((w[m,D]-predict(a7,w[m,]))^2)/M }
MSE=data.frame(MSE);names(MSE)=c("lm","tree","boost","bagging","random forest","nnet","svm")
(NMSE=apply(MSE,2,mean));MSE
#######
扫码加好友，拉您进群
全部版块
我的主页
收藏
