########################################
#TELECOM Saint-Etienne - année 2010/2011
#Cours n°5
#L. Carraro - O. Roustant
########################################

data(iris)
head(iris)
tail(iris)
plot(iris[1:4], pch=c("s", "e", "i")[as.numeric(iris$Species)],col=c("blue","red","green")[as.numeric(iris$Species)])

#résumés numériques
attach(iris) #permet d'utiliser les noms des variables directement
mL<-c(mean(Petal.Length[Species=="setosa"]),mean(Petal.Length[Species=="versicolor"]),mean(Petal.Length[Species=="virginica"]))
medL<-c(median(Petal.Length[Species=="setosa"]),median(Petal.Length[Species=="versicolor"]),median(Petal.Length[Species=="virginica"]))
sdL<-c(sd(Petal.Length[Species=="setosa"]),sd(Petal.Length[Species=="versicolor"]),sd(Petal.Length[Species=="virginica"]))
qL<-cbind(quantile(Petal.Length[Species=="setosa"],probs=c(0.05,0.25,0.75,0.95)),quantile(Petal.Length[Species=="versicolor"],probs=c(0.05,0.25,0.75,0.95)),quantile(Petal.Length[Species=="virginica"],probs=c(0.05,0.25,0.75,0.95)))

#boxplots
par(mfrow=c(1,2))
boxplot(iris[,3]~iris$Species,xlab="Longueur des pétales")
boxplot(iris[,4]~iris$Species,xlab="Largeur des pétales")
par(mfrow=c(1,1))
#on regarde de près une seule boîte
boxplot(iris[,3][iris$Species=="setosa"],main="variété setosa",xlab="Largeur des pétales")

#fonction de répartition empirique
Fr.s<-ecdf(Petal.Length[Species=="setosa"])
Fr.e<-ecdf(Petal.Length[Species=="versicolor"])
Fr.i<-ecdf(Petal.Length[Species=="virginica"])

plot(Fr.s,xlim=c(0.9,7),col="blue",main="fonctions de répartition empiriques",xlab="Longueur des pétales en mm")
lines(Fr.e,col="red")
lines(Fr.i,col="green")
legend("bottomright",col=c("blue","red","green"),lty=1,legend=c("setosa","versicolor","virginica"))

#histogramme (attention à l'argument "freq" pour obtenir un vrai histogramme)
hist.s<-hist(Petal.Length[Species=="setosa"], freq=FALSE, border="blue",xlab="Longueur des pétales en mm",main="Histogramme - variété setosa")
hist.e<-hist(Petal.Length[Species=="versicolor"], freq=FALSE, border="blue",xlab="Longueur des pétales en mm",main="Histogramme - variété versicolor")
hist.i<-hist(Petal.Length[Species=="virginica"], freq=FALSE, border="blue",xlab="Longueur des pétales en mm",main="Histogramme - variété virginica")
#les 3 histogrammes les uns sous les autres
par(mfrow=c(3, 1))
plot(hist.s,main="variété setosa",freq=FALSE, border="blue",xlab="Longueur des pétales en mm",xlim=c(0.9,7))
plot(hist.e,freq=FALSE, border="red",xlab="Longueur des pétales en mm",main="variété versicolor",xlim=c(0.9,7))
plot(hist.i,freq=FALSE, border="green",xlab="Longueur des pétales en mm",main="variété virginica",xlim=c(0.9,7))
#les 3 histogrammes sur la même figure
par(mfrow=c(1,1))
plot(hist.s,main="les 3 histogrammes",freq=FALSE, border="blue",xlab="Longueur des pétales en mm",xlim=c(0.9,7))
lines(hist.e,freq=FALSE, border="red")
lines(hist.i,freq=FALSE, border="green")
legend("topright",col=c("blue","red","green"),lty=1,legend=c("setosa","versicolor","virginica"))
#attention si on représente les effectifs au lieu des fréquences
par(mfrow=c(1,1))
plot(hist.s,main="les 3 histogrammes avec les effectifs",freq=TRUE, border="blue",xlab="Longueur des pétales en mm",xlim=c(0.9,7),ylim=c(0,20))
lines(hist.e,freq=TRUE, border="red")
lines(hist.i,freq=TRUE, border="green")
legend("topright",col=c("blue","red","green"),lty=1,legend=c("setosa","versicolor","virginica"))

#influence du nombre de classes
par(mfrow=c(1, 2))
hist(Petal.Length[Species=="virginica"],freq=FALSE, border="green",xlab="Longueur des pétales en mm",main="variété virginica - 5 classes")
hist(Petal.Length[Species=="virginica"],freq=FALSE, border="green",xlab="Longueur des pétales en mm",main="variété virginica - 12 classes",breaks=seq(from=4,to=7,by=0.25))

#estimations de densité
plot(density(Petal.Length[Species=="setosa"]),col="blue",main="estimations de densité",xlab="Longueur des pétales en mm",xlim=c(0.5,8))
lines(density(Petal.Length[Species=="versicolor"]),col="red")
lines(density(Petal.Length[Species=="virginica"]),col="green")
legend("topright",col=c("blue","red","green"),lty=1,legend=c("setosa","versicolor","virginica"))

#bandwidth=0.1 et 0.05
par(mfrow=c(1,2))
bandwidth=0.1
plot(density(Petal.Length[Species=="setosa"],bw=bandwidth),col="blue",xlab="Longueur des pétales en mm",xlim=c(0.5,8),main=paste("h = ",eval(bandwidth)))
lines(density(Petal.Length[Species=="versicolor"],bw=bandwidth),col="red")
lines(density(Petal.Length[Species=="virginica"],bw=bandwidth),col="green")
legend("topright",col=c("blue","red","green"),lty=1,legend=c("setosa","versicolor","virginica"))
bandwidth=0.05
plot(density(Petal.Length[Species=="setosa"],bw=bandwidth),col="blue",xlab="Longueur des pétales en mm",xlim=c(0.5,8),main=paste("h = ",eval(bandwidth)))
lines(density(Petal.Length[Species=="versicolor"],bw=bandwidth),col="red")
lines(density(Petal.Length[Species=="virginica"],bw=bandwidth),col="green")
legend("topright",col=c("blue","red","green"),lty=1,legend=c("setosa","versicolor","virginica"))

#bandwidth=0.3 et 0.6
par(mfrow=c(1,2))
bandwidth=0.3
plot(density(Petal.Length[Species=="setosa"],bw=bandwidth),col="blue",xlab="Longueur des pétales en mm",xlim=c(0.5,8),main=paste("h = ",eval(bandwidth)))
lines(density(Petal.Length[Species=="versicolor"],bw=bandwidth),col="red")
lines(density(Petal.Length[Species=="virginica"],bw=bandwidth),col="green")
legend("topright",col=c("blue","red","green"),lty=1,legend=c("setosa","versicolor","virginica"))
bandwidth=0.6
plot(density(Petal.Length[Species=="setosa"],bw=bandwidth),col="blue",xlab="Longueur des pétales en mm",xlim=c(0.5,8),main=paste("h = ",eval(bandwidth)))
lines(density(Petal.Length[Species=="versicolor"],bw=bandwidth),col="red")
lines(density(Petal.Length[Species=="virginica"],bw=bandwidth),col="green")
legend("topright",col=c("blue","red","green"),lty=1,legend=c("setosa","versicolor","virginica"))