# Exercice 1
# 1
library(rpart)
?rpart
load("tennis.RData")
tennis
arbre = rpart(Jouer~Ciel+Temperature+Humidite+Vent,
      data=tennis,method="class")
arbre 
# on remarque que par defaut l'arbre n'est pas du tout developpe 
# Les regles d'arret sont atteintes des le debut : principalement,
# l'option minsplit est par defaut fixee a 20 : une branche n'est plus
# developpee des qu'elle contient 20 individus ou moins.
# On fait un nouvel essai avec minsplit=4
options = rpart.control(minsplit=4)
arbre = rpart(Jouer~Ciel+Temperature+Humidite+Vent,
              data=tennis,method="class",control=options)

# 2
arbre
plot(arbre)
text(arbre)

plot(arbre,branch=.3,uniform=TRUE,compress=TRUE,margin=.1)
text(arbre,all=TRUE,use.n=TRUE,fancy=TRUE,minlength = 10)

# 3
predict(arbre,tennis)
cl1 = predict(arbre,tennis,type="class")
sum(cl1!=tennis$Jouer)/14 # taux d'erreur d'apprentissage

# 4
data(iris)
# CART sur les donnees iris
arbre2 = rpart(Species~Sepal.Length+Sepal.Width+Petal.Length+Petal.Width,
               data=iris,method="class",control=options)
# idem mais en ecrivant Species~. pour simplifier la formule (c'est equivalent):
arbre2 = rpart(Species~.,data=iris,method="class",control=options)
# affichage de l'arbre
plot(arbre2,branch=.3,uniform=TRUE,compress=TRUE,margin=.1)
text(arbre2,all=TRUE,use.n=TRUE,fancy=TRUE,minlength = 10)
# utilisation de la fonction prp pour un affichage de l'arbre alternatif :
install.packages("rpart.plot")
library(rpart.plot)
prp(arbre2)
# prediction sur les donnees d'apprentissage :
cl2 = predict(arbre2,iris,type="class")
table(cl2,iris$Species)

# 5
printcp(arbre2)
plotcp(arbre2) 
# on remarque que la valeur du "cp" (complexity parameter) optimale
# est 0.094 (premiere valeur passant sous le seuil en pointilles)
# elagage de l'arbre pour cette valeur :
arbre3 = prune(arbre2,cp=0.094)
# affichage de l'arbre elague :
plot(arbre3,branch=.3,uniform=TRUE,compress=TRUE,margin=.1)
text(arbre3,all=TRUE,use.n=TRUE,fancy=TRUE,minlength = 10)

# 6 (hors programme)
# Utilisation de CART pour la regression
# Les donnees airquality sont quantitatives. On va chercher a 
# predire la variable Ozone en fonction des autres. Attention :
# ce n'est plus un probleme de classification.
data(airquality)
airquality
# On remarque que certaines donnees sont manquantes (NA). On va
# simplement supprimer les individus contenant des donnees 
# manquantes :
donnee = na.omit(airquality)
# construction de l'arbre
options = rpart.control(minsplit=2)
arbre4 = rpart(Ozone~.,data=donnee,method="anova",control=options)
# affichage
prp(arbre4)
# Chaque feuille de l'arbre est labellisee avec la valeur Ozone
# predite par la methode CART.
# remarque 1 : on peut remarquer que les variables Month et Day ne sont
# pas utilisees. C'est tout a fait conforme a l'intuition : ces donnees
# correspondent a des marqueurs temporels et n'ont aucune raison d'etre
# correlees avec la quantite d'ozone.
# remarque 2 : en fait rpart sait egalement gerer les donnees manquantes. 
# On peut donc egalement utiliser directement rpart sur airquality :
arbre5 = rpart(Ozone~.,data=airquality,method="anova")
prp(arbre5)
# Cependant ceci suppose de bien comprendre la façon dont sont gerees
# ces donnees manquantes...
# prediction avec l'arbre construit :
estim_Ozone = predict(arbre5,donnee)
# On  peut evaluer une erreur moyenne relative :
mean(abs(estim_Ozone-donnee$Ozone)/abs(donnee$Ozone))

# Exercice 2
# 1
indtrain = sample(1:150,100)
train = iris[indtrain,]
test = iris[-indtrain,]
# 2
# construction et affichage de l'arbre sur les donnees train :
arbre1 = rpart(Species~.,train,method="class")
plot(arbre1,branch=.3,uniform=TRUE,compress=TRUE,margin=.1)
text(arbre1,all=TRUE,use.n=TRUE,fancy=TRUE,minlength = 10)
# prediction sur les donnees test :
cl1 = predict(arbre1,test,type="class")
# on compare avec les vraies classes :
table(cl1,test$Species)
mean(cl1!=test$Species)
# on trouve 10% d'erreur

# 3
library(randomForest)
?randomForest
# utilisation de randomForest sur les donnees train :
rf = randomForest(Species~.,train)
plot(rf) # affiche des taux d'erreurs pour chacun des 500 arbres construits
# prediction sur les donnees test :
cl2 = predict(rf,test)
table(cl2,test$Species)
mean(cl2!=test$Species)
# on trouve aussi 10% d'erreur. Il faudrait sans doute evaluer mieux l'erreur
# en repetant l'evaluation un grand nombre de fois (ce qui revient plus ou moins
# a faire une validation croisee). On le fera a l'exercice 3, question 5.

# on compare avec CART developpe jusqu'au bout
options = rpart.control(minsplit=1)
arbre3 = rpart(Species~.,train,method="class",control=options)
cl3 = predict(arbre3,test,type="class")
table(cl3,test$Species)
mean(cl3!=test$Species)
# On trouve aussi autour de 10% d'erreur. Les donnees iris sont aussi sans doute 
# trop simples pour observer des differences entre les methodes.

# Exercice 3
# 1
load("fgl.RData")
?fgl
fgl

# 2
indtrain = sample(1:214,200)
train = fgl[indtrain,]
test = fgl[-indtrain,]

# 3
arbre = rpart(type~.,train,method="class")
prp(arbre)

# 4
cl = predict(arbre,test,type="class")
table(cl,test$type)
mean(cl!=test$type)
# environ 29% d'erreur

# 5 (avec 1000 essais plutôt que 100)
E1 = rep(0,1000)
for(i in 1:1000)
{
  indtrain = sample(1:214,200)
  train = fgl[indtrain,]
  test = fgl[-indtrain,]
  arbre = rpart(type~.,train,method="class")
  cl = predict(arbre,test,type="class")
  E1[i] = mean(cl!=test$type)
}
mean(E1) # on trouve environ 30% d'erreur en moyenne

# 6
# On refait l'evaluation en developpant les arbres jusqu'au bout
options = rpart.control(minsplit=1)
E2 = rep(0,1000)
for(i in 1:1000)
{
  indtrain = sample(1:214,200)
  train = fgl[indtrain,]
  test = fgl[-indtrain,]
  arbre = rpart(type~.,train,method="class",control=options)
  cl = predict(arbre,test,type="class")
  E2[i] = mean(cl!=test$type)
}
mean(E2) # on trouve environ 29% d'erreur en moyenne
# Le fait de developper jusqu'au bout donne un tres legere
# amelioration, mais ce n'est pas significatif. Le mode
# par defaut de rpart (utilise a la question 5) donnait
# une erreur presque identique tout en etant plus rapide
# puisque le developpement des arbres n'est pas fait 
# jusqu'au bout. Pour ces donnees qui sont de petite taille
# la question de la rapidite est sans importance, mais pour
# de vrais problemes de classification ça pourrait etre 
# crucial.
# D'autre part on sait qu'en general le developpement de
# l'arbre jusqu'au bout peut impliquer un sur-apprentissage.
# Ici ça ne semble pas vraiment le cas mais il faudrait comparer
# avec une procedure d'elagage de l'arbre.

# 7
bagging = function(train,test,N)
{
  m = nrow(train)
  C = matrix(0,N,nrow(test))
  for(k in 1:N)
  {
    traink = train[sample(1:m,rep=TRUE),]
    arbrek = rpart(type~., data=traink)
    clk = predict(arbrek,test,type="class")
    C[k,] = clk
  }
  apply(C,2,function(x) names(which.max(table(x))))
}

# 8
N = 25
E3 = rep(0,100)
for(i in 1:100)
{
  indtrain = sample(1:214,200)
  train = fgl[indtrain,]
  test = fgl[-indtrain,]
  clf = bagging(train,test,N)
  E3[i] = mean(clf!=as.numeric(test$type))
}
mean(E3) # on trouve environ 27% d'erreur en moyenne
# On observe une legere amelioration des performances en utilisant
# la technique du bagging, mais on verra qu'on peut faire mieux
# avec les forets aleatoires. En fait les forets aleatoires utilisent
# la technique du bagging mais egalement font une selection aleatoire
# des variables utilisees a chaque etape du developpement des
# arbres, ce qui n'est pas fait ici. Ce sont ces deux idees combinees
# qui permettent d'ameliorer les performances de la classification.

# 9
Erf = rep(0,100)
for(i in 1:100)
{
  indtrain = sample(1:214,200)
  train = fgl[indtrain,]
  test = fgl[-indtrain,]
  rf = randomForest(type~.,train)
  cl = predict(rf,test,type="class")
  Erf[i] = mean(cl!=test$type)
}
mean(Erf) 
# on trouve environ 20% d'erreur en moyenne, donc on observe bien
# cette fois une amelioration de la classification avec les forets
# aleatoires pour ces donnees.