2013-2014 Moutaouagad Saïda enregistrement Rplot1.pdf

R version 2.7.2 (2008-08-25)
Copyright (C) 2008 The R Foundation for Statistical Computing
ISBN 3-900051-07-0

# charger la librairie nnet
> library(nnet)
> library(MASS)

# x entre 0 et 10
x <- sort(10*runif(50))
# fonction sin(x) avec un bruit
 y<-sin(x)+0.2*rnorm(x)
 data<-data.frame(x,y)
# Linout=TRUE => données entre 0 et 1
nn<-nnet(x,y,size=4,maxit=100,linout=TRUE)  
# weights:  13
initial  value 60.263858 
iter  10 value 20.190962
iter  20 value 8.576464
iter  30 value 7.490977
iter  40 value 7.353302
iter  50 value 7.001934
iter  60 value 6.250603
iter  70 value 3.522466
iter  80 value 1.902509
iter  90 value 1.497377
iter 100 value 1.481838
final  value 1.481838 
stopped after 100 iterations

> plot(x,y,col="blue",pch=16)

> x<-seq(0,10,by=0.01)

> lines(x1,sin(x1),col="green")


> lines(x1,predict(nn,data.frame(x=x1)),col="red")


> 

Call:
lines(x1, predict(nn, data.frame(x = x1)))

Coefficients:
[1]   0.68504  -0.03743


> modelPoly <- lm(y ~ x + I(x^2) + I(x^3) + I(x^4))
Erreur dans model.frame.default(formula = y ~ x + I(x^2) + I(x^3) + I(x^4),  : 
  les longueurs des variables diffèrent (trouvé pour 'x')
Avec la nouvelle version de R cela fonctionne :
> modelPoly <- lm(y ~ x + I(x^2) + I(x^3) + I(x^4))

lines(x,predict(modelPoly), col="orange") 

 > legend(3,1.3,c("données simulées","fonct","PMC","modèle polynomial"), lty=c(0,1,1,1),pch=c(16,-1,-1-1),col=c("blue","green","red","orange"))

enregistrements  Rplot1.pdf et PMCnnet_1.pdf

====
Soulier 2012-2013

> appindex = sample(1:nrow(iris), round(2*nrow(iris)/3), replace = FALSE)
> app = iris[appindex,]
> val = iris[-appindex,]
> library(nnet)
> model.dis = nnet(Species ~ ., data = app, size = 2, decay = 0.001)
# weights:  19
initial  value 121.270720 
iter  10 value 50.414138
iter  20 value 12.508239
iter  30 value 10.669142
iter  40 value 10.530162
iter  50 value 9.310038
iter  60 value 8.379050
iter  70 value 7.797647
iter  80 value 7.703485
iter  90 value 6.812072
iter 100 value 5.872002
final  value 5.872002 
stopped after 100 iterations
> summary(model.dis)
a 4-2-3 network with 19 weights
options were - softmax modelling  decay=0.001
 b->h1 i1->h1 i2->h1 i3->h1 i4->h1 
  3.87   0.64   0.69  -1.08  -2.72 
 b->h2 i1->h2 i2->h2 i3->h2 i4->h2 
 -0.46  -0.38  -3.53   5.19   2.40 
 b->o1 h1->o1 h2->o1 
 -7.37  15.93  -9.88 
 b->o2 h1->o2 h2->o2 
 -1.46   2.51   4.90 
 b->o3 h1->o3 h2->o3 
  8.83 -18.44   4.98 

===

Remarquons que notre modèle a été ajusté sur la base d'apprentissage et que
pour avoir une idée de sa performance, on peut calculer le taux d'erreur de
classement lorsque la prédiction est réalisée sur le corpus de validation :

> pred = predict(model.dis, newdata = val[, -5], type = "class")
> mat = table(pred, val[,5])
> taux = sum(diag(mat))/sum(mat)
> mat
            
pred         setosa versicolor virginica
  setosa         19          0         0
  versicolor      0         15         0
  virginica       0          1        15
> taux
[1] 0.98

val[,-5] : toutes les colonnes de val, sauf la 5 ième, celle que l'on désire
prédire avec predict et notre réseau model.dis

<\session|r|default>
  <\input-math>
    <with|color|red|\<gtr\> >
  <|input-math>
    val[, -5]
  </input-math>
</session>

<math|val[,5]> : la dernière colonne\ 

<\session|r|default>
  <\unfolded-io-math>
    <with|color|red|\<gtr\> >
  <|unfolded-io-math>
    val[, 5]
  <|unfolded-io-math>
    \ [1] setosa \ \ \ \ setosa \ \ \ \ setosa \ \ \ \ setosa \ \ \ \ setosa
    \ \ \ \ setosa \ \ \ 

    \ [7] setosa \ \ \ \ setosa \ \ \ \ setosa \ \ \ \ setosa \ \ \ \ setosa
    \ \ \ \ setosa \ \ \ 

    [13] setosa \ \ \ \ setosa \ \ \ \ setosa \ \ \ \ versicolor versicolor
    versicolor

    [19] versicolor versicolor versicolor versicolor versicolor versicolor

    [25] versicolor versicolor versicolor versicolor versicolor versicolor

    [31] versicolor versicolor versicolor virginica \ virginica \ virginica\ 

    [37] virginica \ virginica \ virginica \ virginica \ virginica
    \ virginica\ 

    [43] virginica \ virginica \ virginica \ virginica \ virginica
    \ virginica\ 

    [49] virginica \ virginica\ 

    Levels: setosa versicolor virginica
  </unfolded-io-math>

  <\input-math>
    <with|color|red|\<gtr\> >
  <|input-math>
    \;
  </input-math>
</session>

<\session|r|default>
  <\input-math>
    <with|color|red|\<gtr\> >
  <|input-math>
    pred = predict(model.dis, newdata = val[, -5], type = "class")
  </input-math>
</session>

Voyons notre prédiction et comparons-là ensuite dans la matrice mat avec la
réalité

<\session|r|default>
  <\unfolded-io-math>
    <with|color|red|\<gtr\> >
  <|unfolded-io-math>
    pred
  <|unfolded-io-math>
    \ [1] "setosa" \ \ \ \ "setosa" \ \ \ \ "setosa" \ \ \ \ "setosa"
    \ \ \ \ "setosa" \ \ \ 

    \ [6] "setosa" \ \ \ \ "setosa" \ \ \ \ "setosa" \ \ \ \ "setosa"
    \ \ \ \ "setosa" \ \ \ 

    [11] "setosa" \ \ \ \ "setosa" \ \ \ \ "setosa" \ \ \ \ "setosa"
    \ \ \ \ "setosa" \ \ \ 

    [16] "versicolor" "versicolor" "versicolor" "versicolor" "versicolor"

    [21] "versicolor" "versicolor" "versicolor" "versicolor" "versicolor"

    [26] "versicolor" "versicolor" "versicolor" "versicolor" "versicolor"

    [31] "versicolor" "versicolor" "versicolor" "versicolor" "virginica"\ 

    [36] "virginica" \ "virginica" \ "virginica" \ "virginica" \ "virginica"\ 

    [41] "virginica" \ "virginica" \ "versicolor" "virginica" \ "virginica"\ 

    [46] "virginica" \ "virginica" \ "virginica" \ "virginica" \ "virginica"\ 
  </unfolded-io-math>

  <\input-math>
    <with|color|red|\<gtr\> >
  <|input-math>
    \;
  </input-math>
</session>

<\session|r|default>
  <\input-math>
    <with|color|red|\<gtr\> >
  <|input-math>
    mat = table(pred, val[,5])
  </input-math>

  <\input-math>
    <with|color|red|\<gtr\> >
  <|input-math>
    taux = sum(diag(mat))/sum(mat)
  </input-math>

  <\unfolded-io-math>
    <with|color|red|\<gtr\> >
  <|unfolded-io-math>
    mat
  <|unfolded-io-math>
    \ \ \ \ \ \ \ \ \ \ \ \ 

    pred \ \ \ \ \ \ \ \ setosa versicolor virginica

    \ \ setosa \ \ \ \ \ \ \ \ 15 \ \ \ \ \ \ \ \ \ 0 \ \ \ \ \ \ \ \ 0

    \ \ versicolor \ \ \ \ \ 0 \ \ \ \ \ \ \ \ 18 \ \ \ \ \ \ \ \ 2

    \ \ virginica \ \ \ \ \ \ 0 \ \ \ \ \ \ \ \ \ 0 \ \ \ \ \ \ \ 15
  </unfolded-io-math>
</session>

A ``l'oeil nu'', il y a deux erreurs : deux virginica estimées versicolor ;
ce que confirme la matrice mat. Notre taux d'erreur est :

<\session|r|default>
  <\unfolded-io-math>
    <with|color|red|\<gtr\> >
  <|unfolded-io-math>
    taux
  <|unfolded-io-math>
    [1] 0.96
  </unfolded-io-math>
</session>

Pour notre première tentative, nous avons pris des paramètres par défauts. A
présent, l'objectif va être de trouver le size et le decay optimal pour notre
réseau de neurones. En effet, toute méthode d'apprentissage a ses propres
paramètres que l'utilisateur doit ajuster pour construire un bon modèle à la
fois adapté aux données de la base d'apprentissage et performant du point de
vue de la prédiction sur les nouvelles observations. Les 2 principaux
paramètres de notre modèle sont :\ 

- <with|color|blue|size> : le nombre de neurones de la couche cachée\ 

- <with|color|blue|decay> : le pas d'apprentissage (paramètre de
décomposition)

En faisant varier ces deux paramètres, on s'aperçoit que plus le modèle est
complexe (plus <with|color|blue|size> est grand), mieux il apprend sur les
données et plus il est performant du point de vue de la prédiction. On peut
également remarquer qu'au-delà d'un certain seuil, la performance n'augmente
plus, voire diminue, c'est le sur-apprentissage. En faisant varier le
<with|color|blue|decay>, nous constatons que plus la décomposition de nos
neurones est importante, plus l'est la performance et réciproquement, s'il y
a trop de décompositions, le modèle redevient moins performant. La
problématique consiste à trouver le bon paramétrage du modèle.

A cette fin, nous utilisons la fonction tune du package e1071 pour essayer de
répondre à cette question : library(e1071) ; il vous faudra probablement
installer ce paquetage (voir onglet ``Paquetages & Données))

> library(e1071)
> tune.model = tune.nnet(Species  ~ ., data = app, size = c(1, 3, 5), decay = c(0.1, 0.001, 0.000001))
> plot(tune.model) SAUVE dans PlotTune*.pdf
> tune.model

Parameter tuning of ‘nnet’:

- sampling method: 10-fold cross validation 

- best parameters:
 size decay
    3   0.1

- best performance: 0.03 

Ajustons notre réseau de neurones avec les paramètres récemment trouvés, puis regardons la performance de prédiction grâce à la matrice de confusion :

model = nnet(Species ~ ., data = app, size = 1, decay = 0.1, maxit = 100)
summary(model)
pred = predict(model, newdata = val[, -5], type = "class")
mat = table(pred, val[,5])
taux = sum(diag(mat))/sum(mat)

> model = nnet(Species ~ ., data = app, size = 3, decay = 0.1, maxit = 100)
# weights:  27
initial  value 124.850169 
iter  10 value 35.065507
iter  20 value 23.019100
iter  30 value 21.517937
iter  40 value 21.077388
iter  50 value 19.543053
final  value 19.542882 
converged
> summary(model)
a 4-3-3 network with 27 weights
options were - softmax modelling  decay=0.1
 b->h1 i1->h1 i2->h1 i3->h1 i4->h1 
 -0.21  -0.33  -1.10   1.83   0.83 
 b->h2 i1->h2 i2->h2 i3->h2 i4->h2 
 -1.76  -1.62  -1.40   1.95   3.77 
 b->h3 i1->h3 i2->h3 i3->h3 i4->h3 
  0.21   0.34   1.06  -1.70  -0.79 
 b->o1 h1->o1 h2->o1 h3->o1 
  0.50  -2.51  -0.90   2.99 
 b->o2 h1->o2 h2->o2 h3->o2 
  0.57   2.62  -3.85  -1.91 
 b->o3 h1->o3 h2->o3 h3->o3 
 -1.07  -0.11   4.75  -1.09 
> pred = predict(model, newdata = val[, -5], type = "class")
> mat = table(pred, val[,5])
> taux = sum(diag(mat))/sum(mat)
> mat
            
pred         setosa versicolor virginica
  setosa         19          0         0
  versicolor      0         15         1
  virginica       0          1        14
> taux
[1] 0.96

A présent, on peut réaliser une discrimination sur le même jeu de données. Nous allons utiliser le package rpart. On commence par réaliser le même découpage avec une partie pour l’apprentissage et la validation. Puis on calcule notre modèle et on trace l’arbre de décision :

library(rpart)
model = rpart(Species ~ ., data = app, method = "class")
plot(model, uniform = TRUE, branch = 0.5, margin = 0.1,main="Arbre de discrimination")
text(model, all = FALSE, use.n = TRUE)

Ensuite, comme pour les réseaux de neurones, nous allons rechercher les paramètres optimaux pour la fonction rpart. Ici, nous avons à faire varier les paramètres suivants : - minsplit : correspond au nombre minimal d'observations d'un noeud. - cp : est le coefficient de complexité du modèle. - xval : il s’agit du nombre de validations croisées - maxdepth : est la profondeur maximale de l'arbre, sachant qu'à la racine, la profondeur vaut 0.

tune.model = tune.rpart(Species ~ ., data = app, minsplit = c(15, 20, 25),
cp = c(0.00001, 0.000001), maxcompete = 4, maxsurrogate = 5,
usesurrogate = 2, xval = c(10, 15), surrogatestyle = 0, maxdepth = c(25, 30))

> tune.model

Parameter tuning of ‘rpart.wrapper’:

- sampling method: 10-fold cross validation 

- best parameters:
 minsplit    cp maxcompete maxsurrogate usesurrogate xval surrogatestyle maxdepth
       15 1e-05          4            5            2   10              0       25

- best performance: 0.04 

On peut déjà remarquer qu?en utilisant le PCM, la performance est meilleure
(0.03) qu'à l'utilisation de l'arbre de décision (best performance : 0.09).
Ainsi nous pouvons utiliser la fonction avec les paramètres optimaux :

model = rpart(Species ~ ., data = app, method = "class",
control = rpart.control(cp = 0.00001, minsplit = 15, maxcompete = 4,
maxsurrogate = 5, usesurrogate = 2, xval = 10, surrogatestyle = 0, maxdepth = 25))

L'arbre qui nous intéresse est celui qui minimise l’erreur standard (xerror + xstd dans notre modèle). Si plusieurs arbres minimisent cette quantité, on privilègie l'arbre le plus petit. Ainsi il faut réaliser un élagage puis on peut afficher l’arbre de décision :
cpTab = as.data.frame(model$cptable)
ind = cpTab$xerror + cpTab$xstd
cpTab = cbind.data.frame(cpTab, ind)
monCp = cpTab[which.min(cpTab$ind),"CP"]
elag = prune(model, cp = monCp)
#Arbre
plot(elag, uniform = TRUE, branch = 0.5, margin = 0.1,main="Arbre de décision 2")
text(elag, all = FALSE, use.n = TRUE, cex = 0.6)

# Sauvé dans un fichier pdf comme le premier arbre ArbreDec_2.pdf

Enfin, nous allons utiliser notre arbre élagué pour la prédiction de notre base de validation et construire la matrice de confusion :
prev = predict(elag, newdata = val[,-5], type = "class")
mat = table(prev, val[,5])
taux = sum(diag(mat))/sum(mat)

> mat 
            
prev         setosa versicolor virginica
  setosa         19          0         0
  versicolor      0         15         2
  virginica       0          1        13
> taux
[1] 0.94

taux