##################################### #regression - option maths applis #EMSE #L. Carraro (Telecom Saint-Etienne) #decembre 09 ##################################### ####################################### #boostrap et bagging ####################################### #exemple de la loi normale N(m,1) pour le bootstrap ################################################### n<-100 alpha<-0.05 x<-rnorm(n,mean=5,sd=1) #IC théorique a<-mean(x)+(1/sqrt(n))*qnorm(alpha/2) b<-mean(x)+(1/sqrt(n))*qnorm(1-(alpha/2)) #IC par bootstrap #NB : le package boot est beaucoup plus sophistiqué B<-5000 ech.boot<-vector(length=B) for (i in 1:B) { temp<-sample(x,size=n,replace=T) ech.boot[i]<-mean(temp) } a.boot<-as.numeric(quantile(ech.boot,probs=alpha/2)) b.boot<-as.numeric(quantile(ech.boot,probs=1-alpha/2)) c(a,a.boot) c(b,b.boot) #bagging pour kNN ################# #on montre que c'est inutile car le modèle est linéaire en les réponses B<-50 load(file="classification.Rdata") library(class) #plan de test pour graphique n<-200 x<-seq(-3,4,length.out=n) y<-seq(-2,3,length.out=n) plan<-expand.grid(x=x,y=y) class.b1<-matrix(0,nrow=n,ncol=n) for (i in 1:B) { order<-sample(1:n,size=n,replace=T) temp<-data.class$ech[order,] class.knn<-as.numeric(knn(train=temp[,1:2],test=plan,cl=temp[,3],k=1)) class.b1<-matrix(class.knn,nrow=n,ncol=n)+class.b1 } indic.knn.b1<-(class.b1>3*B/2) ech<-data.class$ech plot(ech[ech[,3]==1,],col="red",xlim=c(min(ech[,1]),max(ech[,1])),ylim=c(min(ech[,2]),max(ech[,2])),xlab="x",ylab="y") title(main=c("séparation des deux classes","kNN avec k=1")) points(ech[ech[,3]==0,],col="blue") points(t(data.class$m1),pch=23,col="red",bg="red") points(t(data.class$m2),pch=23,col="blue",bg="blue") contour(x,y,indic.knn.b1,add=T,lwd=2,nlevels=2) class.knn1<-knn(train=data.class$ech[,1:2],test=plan,cl=data.class$ech[,3],k=1) indic.knn1<-matrix(class.knn1,nrow=n) contour(x,y,indic.knn1,add=T,lty=3,lwd=1,col="orange",nlevels=2) legend("topright",cex=0.8,legend=c("frontière kNN","frontière après bagging : B=50"),lty=c(1,3),lwd=c(2,1),col=c("black","orange")) ############################################# #exemple synthetique de predicteurs correles ############################################# #plan d'exp a 30 points #modele y=2(x1+x2+x3)+N(0,0.1) #x1 N(0,1), x2/x1 N(3x1,0.1),x3/x1,x2 N(x1-x2,0.1) n<-30 x1<-rnorm(n,0,1) x2<-rnorm(n,3*x1,0.1) x3<-rnorm(n,x1-x2,0.1) y<-rnorm(n,2*(x1+x2+x3),0.5) pairs(cbind(x1,x2,x3,y)) plot(x1+x2+x3,y) kappa(cbind(x1,x2,x3)) mod1<-lm(y~x1+x2+x3) summary(mod1) #le choix de modele qui suit depend evidemment des simulations mod2<-lm(y~x2+x3-1) anova(mod2,mod1) summary(mod2) #plan test plan.test<-expand.grid(x1=seq(from=-2,to=2,length=11),x2=seq(from=-6,to=6,length=11),x3=seq(from=-4,to=4,length=11)) rep.test<-rnorm(2*(plan.test[,1]+plan.test[,2]+plan.test[,3]),0.5) pred2<-predict(mod2,as.data.frame(plan.test),interval="prediction") #ecart type est ecarts ecart2<-sqrt(sum((pred2[,1]-rep.test)^2)/1331) #proportion de points dans l'intervalle de prevision prop.test2<-mean((rep.test>=pred2[,2])&(rep.test<=pred2[,3])) ecart2 prop.test2 #idem pour le modele complet pred1<-predict(mod1,as.data.frame(plan.test),interval="prediction") #ecart type est ecarts ecart1<-sqrt(sum((pred1[,1]-rep.test)^2)/1331) #proportion de points dans l'intervalle de prevision prop.test1<-mean((rep.test>=pred1[,2])&(rep.test<=pred1[,3])) ecart1 prop.test1 save(x1,x2,x3,y,rep.test,file="simu_pred_correles.RData") save(x1,x2,x3,y,rep.test,file="simu_pred_correles.txt",ascii=TRUE) permut<-sample(1331,1331) plot(pred1[permut,1],col="blue",pch=19,ylim=c(-70,70),main="prévisions et intervalles de confiance") points(pred1[permut,3],pch=4,cex=0.8,col="green") points(pred1[permut,2],pch=3,cex=0.8,col="green") legend("topleft",cex=0.8,legend=c("prévision","borne inf prévision","borne sup prévision"),pch=c(19,3,4),col=c("blue","green","green")) #################################################### #ridge et lasso #librement inspire des exs du packages ElemStatLearn ##################################################### library(ElemStatLearn) library(leaps) str( prostate ) cor( prostate[,1:8] ) #on separe les donnes d'apprentissage de celles de test train <- subset( prostate, train==TRUE )[,1:9] test <- subset( prostate, train=FALSE )[,1:9] pairs( train[,1:9], col="violet" ) #on fait les études avec les donnees d'apprentissage prostate.leaps <- regsubsets( lpsa ~ . , data=train, nbest=70, really.big=TRUE )#on teste tous les sous-modeles prostate.leaps.sum <- summary( prostate.leaps ) prostate.models <- prostate.leaps.sum$which prostate.models.size <- as.numeric(dimnames(prostate.models)[[1]]) barplot( table(prostate.models.size),xlab="taille modèle",ylab="effectifs",main="histogramme des tailles de modèles" ) prostate.models.rss <- prostate.leaps.sum$rss prostate.models.best.rss <- tapply( prostate.models.rss, prostate.models.size, min ) prostate.models.best.rss #on ajoute le modèle sans prédicteur, qui donne la variance des donnees prostate.0 <- lm( lpsa ~ 1, data=train ) prostate.models.best.rss <- c( sum(resid(prostate.0)^2), prostate.models.best.rss) # Making a plot: plot( 0:8, prostate.models.best.rss, ylim=c(20, 100), type="b", xlab="taille modèle", ylab="Residual Sum Square",main="régression sur tous les sous-modèles", col="red" ) points( prostate.models.size, prostate.models.rss, pch=4, col="blue",cex=0.5 ) points(0,sum(resid(prostate.0)^2),col="blue",pch=4,cex=0.5) # ridge regression: prostate.ridge <- simple.ridge( train[,1:8], train[,9], df=seq(from=0.5,by=0.5,to=8) ) # coefficients: matplot( prostate.ridge$df, t(prostate.ridge$beta), type="b", col="blue", pch=17, cex=0.8,xlab="degrés de liberté",ylab="coefficients" ,main="régression ridge") # Calculations for the lasso: # library(lasso2) prostate.lasso <- l1ce( lpsa ~ ., data=train, trace=TRUE, sweep.out=~1, bound=seq(0,1,by=0.1) ) prostate.lasso.coef <- sapply(prostate.lasso, function(x) x$coef) colnames(prostate.lasso.coef) <- seq( 0,1,by=0.1 ) matplot( seq(0,1,by=0.1), t(prostate.lasso.coef[-1,]), type="b", xlab="borne L1 sur les coefficients", ylab="coefficients", main="régression LASSO", xlim=c(0, 1), col="blue", pch=17 )