############################################################################ ##### ##### Correction du TP2 - Statistiques Descriptives ##### ############################################################################ rm(list=ls()) # Ne pas mettre correctement les commentaires de telle sorte que ca ne tourne pas quand on fait un copie colle faisait perdre des points. # En particulier >B=as.matrix(B) # ne marche pas a cause du > # Attention aux accents et aux ' qui veulent dire quelquechose en R, ca ne tourne pas parfois a cause de ca. # Par principe, ne PAS METTRE D'ACCENTS ni toutes autres ponctuations si ce n'est pas en commentaire (et meme). # ----------------------------------------------------------------------------------- # 1. x= seq(0,10,0.1) # Prendre un pas suffisament petit pour obtenir un trace assez fin. plot(x, dexp(x,1), ylim=c(0,5), type='l', col='red', ylab='') par(new=TRUE) plot(x, dexp(x,2), ylim=c(0,5), type='l', col='blue', ylab='') par(new=TRUE) plot(x, dexp(x,4), ylim=c(0,5), type='l', col='green', ylab='') par(new=TRUE) plot(x, dexp(x,5), ylim=c(0,5), type='l', col='magenta', ylab='f(x)', main='Densite exponentielle') legend(6,4,c('lambda=1','lambda=2','lambda=4','lambda=5'),col=c('red','blue','green','magenta'),lty=c(1,1,1,1)) # Ne pas negliger de mettre les noms aux axes et les titres ou la legende quand necessaire. # 2. help(hist) # hist cree un histogramme. help(barplot) # barplot cree un diagramme en baton. help(table) # Calcule la table de contingence. # 3. # 4. setwd('...') # A vous de mettre le bon chemin vers votre repertoire. A= read.table('exo1b.txt') B= read.table('exo1b.txt', header=TRUE) C= read.table('exo2b.txt') D= read.table('exo2b.txt', header=TRUE) # read.table lit les donnees. # Comme il y a un entete sur exo1b.txt, il faut mettre la commande header=TRUE pour recuperer uniquement les donnees, 'attente' (le nom des donnees) etant alors mis comme nom de colonnes. # exo2b.txt n'a pas d'entete, si on prenait D on perdrait une donnee. # Les bonnes acquisitions de donnees sont donc B et C. # 5. is.data.frame(A) # A est un data frame, c'est toujours cette classe qui est renvoyee par read.table. # class(A) permet de savoir directement ce que c'est. class(A) # NB : un data.frame est un cas particulier de liste donc is.list(A) # renvoie TRUE egalement. # 6. B= as.matrix(B) # 7. Les donnees de B etant de nature quatitative continue, la representation adaptee est l'histogramme. # Remarquez que si vous n'aviez pas transforme B en matrice, vous ne pourriez pas lui appliquer hist, qui ne prend pas en entree des data.frame. hist(B) # Ne donne pas la bonne representation graphique car il n'y a que les effectifs en ordonnee et pas les hauteurs renormalisees de l'histogramme correct. hist(B,freq=FALSE) # Ici on represente bien les bonnes hauteurs. # NB: par defaut, hist utilise la methode de Sturges (vue en cours) pour determiner le nombre de classes. # Attention, ce graphe n'est pas suffisant, il faut les bons titres. hist(B, freq=FALSE, main='histogramme des donnees exo1b.txt', xlab='attente', ylab='densite') # 8. hist(B, main='histogramme des donnees exo1b.txt', xlab='attente', ylab='densite', breaks=c(0,1,2,3,5,7,10,15,20)) # Si on avait mis breaks=k, avec k un nombre entier, cela fixerait le nombre de classes qui seraient alors de meme taille. # breaks=v, avec v est un vecteur, donne les limites entre les classes. # NB : si on precise breaks, on n'a plus besoin de preciser freq=FALSE, c'est fait par defaut. # 9. # Une modelisation par variable exponentielle semble adaptee. # Il faut estimer le bon parametre lambda de la loi. # On sait qu'une variable exponentielle a comme esperance 1/lambda, donc l'inverse de la moyenne empirique est un bon estimateur de lambda. n= nrow(B) moy= sum(B)/n # Ou plus directement : moy=mean(B). lambdahat= 1/moy # Ou plus directement : lambdahat=1/mean(B). # 10. x= seq(0,20,0.1) y= dexp(x, lambdahat) # On calcule la densite exponentielle associee au lambda estime. par(new=TRUE) # On garde l'histogramme precedent. plot(x, y, type='l', col='red', ylim=c(0,0.25), xlim=c(0,20), xlab='', ylab='') # On regarde le graphe precedent pour avoir les limites en x et y et on les remet pour avoir un graphe qui ne change pas d'echelle, on enleve les titres des axes pour ne pas polluer le graphe precedent. # Une solution beaucoup plus propre (mais un peu plus gourmande en temps) est la suivante : # Pour pouvoir connaitre a l'avance les valeurs a mettre dans xlim ylim, on peut faire la chose suivante. H= hist(B,plot=FALSE,breaks=c(0,1,2,3,5,7,10,15,20)) # On connait les limites de l'axe des abscisses car on impose les classes. # Par contre, on ne connait pas la hauteur maximale de l'histogramme. # Pour recuperer cette information, on donne un nom a la commande hist et on ajoute plot=FALSE pour que l'histogramme ne soit pas dessine. ym1= max(H$density) # Cela permet de recuperer la hauteur maximale de l'histogramme car H$density est le vecteur contenant les hauteurs de chaque classe. # Maintenant on prend le maximum entre la hauteur maximale de l'histogramme et les ordonnees de la fonction de densite. ymax= max(y,ym1)*1.01 # La multiplication par 1.01 est presente simplement pour avoir une petite marge au dessus du point le plus haut sur les dessins. # On trace maintenant l'histogramme avec la densite en ordonnees, et avec la definition des axes en terme de limites. hist(B, freq=FALSE, breaks=c(0,1,2,3,5,7,10,15,20), xlim=c(0,20), ylim=c(0,ymax), xlab='attente', ylab='densite', main='histogramme et theorie') par(new=TRUE) plot(x, y, type='l', col='red', xlim=c(0,20), ylim=c(0,ymax), xlab='', ylab='') # On rajoute une legende. legend(10, 0.2, c('observe','theorique'), lty=c(1,1), col=c('black','red')) # Il y a aussi des solutions "economiques" avec lines et curve. # Pour lines voir ci-dessous. # Pour curve, si vous voulez l'utiliser, a vous de bien lire la doc. # 11. tc= table(C) # Construit la table de contingence. tc= tc/nrow(C) # Passage des effectifs aux frequences. # Construction du diagramme en baton avec les frequences : barplot(tc, main='Diagramme en baton des donnees exo2b.txt', ylab='frequence', xlab='valeurs') # 12. # Il n'y a que deux valeurs possibles, 0 et 1. C'est une variable de Bernoulli. # Pour estimer son parametre qui est la proba d'avoir un 1, on compte le nombre de 1 dans l'echantillon et on divise par la taille de l'echantillon, c'est en fait deja calcule. phat= mean(as.matrix(C)) # Attention C est un data.frame il faut le transformer en matrice pour pouvoir appliquer ce type de fonction. # On peut eventuellement dire que c'est une binomiale de parametres (1,p). # Mais pas une binomiale (n,p) (ca ce serait plutot la somme des valeurs de C). # 13. h= hist(B, breaks=c(0,1,2,3,5,7,10,15,20), plot=FALSE) h # Cela renvoie une liste avec les coupures, les effectifs par classe, les hauteurs de classe, le milieu des classes, le nom de l'axe (cad le nom des donnees par defaut), et si les classes sont de meme taille. # class(h) dit que c'est un histogramme mais is.list(h) dit que c'est une liste. C'est parce que histogramme est un cas particulier de list, comme data.frame. eff= h$counts # On recupere les effectifs. # 14. cumsum(eff) # Fait la somme cumulee des effectifs, ce qui va permettre de tracer la courbe des frequences cumulees. # 15. # Ce sont des variables quantitatives donc la courbe des frequence cumulees doit etre une ligne brisee. x= c(0,1,2,3,5,7,10,15,20) # Les abscisses ou l'on va calculer la ligne. y= c(0,cumsum(eff)/nrow(B)) # Les ordonnees correspondantes, cumsum fait la somme cumulee. # NB : on rajoute un zero au debut de y afin que x et y soit de meme taille. plot(x, y, type='l', main='Courbe des frequences cumulees', xlab='attente', ylab='Frequences cumulees', xaxt='n') axis(1,x) # xaxt='n' pour enlever les graduations de l'axe des abscisses, axis(1,x) pour mettre le vecteur x sur la coordonnee 1 (ie les abscisses) en tant que graduation de l'axe. # barplot est une mauvaise representation !!! # 16. xexp= seq(0,20,0.1) yexp= pexp(xexp, lambdahat) # pexp donne la fonction de repartition de la variable exponentielle. lines(xexp,yexp,col='red') # lines permet de superposer une ligne sans faire par(new=TRUE). # Ne modifie pas les graduations des axes et conserve les xlim et ylim du premier graphe. legend(10, 0.6, c('observe','theorique'), lty=c(1,1), col=c('black','red')) # 17. par(mfrow=c(1,2)) xexp= seq(0,20,0.1) yexp_densite= dexp(xexp, lambdahat) hist(B, main='histogramme des donnees exo1b.txt', xlab='attente', ylab='Densite', breaks=c(0,1,2,3,5,7,10,15,20)) lines(xexp, yexp_densite, col='red') legend(10, 0.2, c('observe','theorique'), lty=c(1,1), col=c('black','red')) plot(x, y,type='l', main='Courbe des frequences cumulees', xlab='attente', ylab='Frequences cumulees', xaxt='n') axis(1, x) lines(xexp, yexp, col='red') legend(10, 0.6, c('observe','theorique'), lty=c(1,1), col=c('black','red')) # 18. x= rnorm(1000) # Boxplot permet de tracer la boite a moustache d'une serie de donnees. boxplot(x, main='Boite a moustache pour x') # Le trait gras est la mediane de x, les deux traits adjacents sont les 1er et 3eme quartiles. # Les traits exterieurs sont les quantiles d'ordre 0.005 et 0.995. # Les cercles sont les valeurs extremes. # 19. y= rnorm(1000, 5, 2) boxplot(x, y, main='Comparaison des distributions de x et y', names=c('x', 'y')) # Utiliser names pour preciser le nom de chaque boite a moustache. # Ce graphique permet d'illustrer la difference de distribution entre x et y. # Celle de y est decalee vers le haut (plus grande moyenne) et est plus etalee (plus grande variance).