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

pluie <- read.table("pluie.txt", header=TRUE, dec=",",sep="\t")
with(data=pluie, stripchart(hauteur~traitement, xlab="hauteur de pluie (mm)", ylab="traitement"))
#si on ne connait pas les stripchart, on peut faire à la main :
plot(pluie$hauteur,pluie$traitement,xlab="hauteur de pluie (mm)", ylab="traitement")

#premiers indicateurs
hauteur.nt<-pluie$hauteur[pluie$traitement==0]
hauteur.t<-pluie$hauteur[pluie$traitement==1]
summary(hauteur.nt); sd(hauteur.nt); quantile(hauteur.nt,c(0.05,0.95));quantile(hauteur.nt,0.75)-quantile(hauteur.nt,0.25)
summary(hauteur.t); sd(hauteur.t); quantile(hauteur.t,c(0.05,0.95));quantile(hauteur.t,0.75)-quantile(hauteur.t,0.25)

#on passe aux log
with(data=pluie, stripchart(log(hauteur)~traitement, xlab="logarithme de la pluviométrie", ylab="traitement"))
lhauteur.nt<-log(hauteur.nt)
lhauteur.t<-log(hauteur.t)
summary(lhauteur.nt); sd(lhauteur.nt); quantile(lhauteur.nt,c(0.05,0.95));quantile(lhauteur.nt,0.75)-quantile(lhauteur.nt,0.25)
summary(lhauteur.t); sd(lhauteur.t); quantile(lhauteur.t,c(0.05,0.95));quantile(lhauteur.t,0.75)-quantile(lhauteur.t,0.25)

#boxplots
with(data=pluie, boxplot(log(hauteur)~traitement, horizontal=TRUE, range=1,                         xlab="logarithme de la pluviométrie", ylab="traitement"))

#fonction de répartition empirique
FRht <- ecdf(lhauteur.t); FRhnt <- ecdf(lhauteur.nt)
abs <- seq(-3,6,0.1)
plot(abs, FRht(abs), type="s", col="blue", ylim=c(0,1), main="fonctions de répartition empiriques", xlab="logarithme de la pluviométrie", ylab="probabilité")
lines(abs, FRhnt(abs), type="s", col="red")
legend("topleft",col=c("red","blue"),lty=1,legend=c("non traités","traités"))

#histogramme (attention à l'argument "freq" pour obtenir un vrai histogramme)
hist(lhauteur.t, freq=FALSE, border="blue",xlab="log pluviométrie",main="nuages traités") 
#ajout des donnees sur le graphique
points(lhauteur.t, rep(0,length(lhauteur.t)), pch="I", cex=1.5,col="green", xlim=c(-2,6), ylim=c(0,0.4))

#les deux jeux de données
op <- par(mfrow=c(2, 1))
hist(lhauteur.t, freq=FALSE, border="blue",xlab="log pluviométrie",main="nuages traités") 
points(lhauteur.t, rep(0,length(lhauteur.t)), pch="I", cex=1.5,col="green", xlim=c(-2,6), ylim=c(0,0.4))
hist(lhauteur.nt, freq=FALSE, border="red",xlab="log pluviométrie",main="nuages non traités") 
points(lhauteur.nt, rep(0,length(lhauteur.nt)), pch="I", cex=1.5,col="green", xlim=c(-2,6), ylim=c(0,0.4))
par(op)

#influence du nombre de classes
op <- par(mfrow=c(1, 2))
hist(lhauteur.t,4, freq=FALSE, border="blue",xlab="log pluviométrie",main="nuages traités - 4 classes") 
points(lhauteur.t, rep(0,length(lhauteur.t)), pch="I", cex=1.5,col="green", xlim=c(-2,6), ylim=c(0,0.4))
hist(lhauteur.t,12, freq=FALSE, border="blue",xlab="log pluviométrie",main="nuages traités - 12 classes") 
points(lhauteur.t, rep(0,length(lhauteur.t)), pch="I", cex=1.5,col="green", xlim=c(-2,6), ylim=c(0,0.4))
par(op)

op <- par(mfrow=c(1, 2))
hist(lhauteur.nt,4, freq=FALSE, border="red",xlab="log pluviométrie",main="nuages non traités - 4 classes") 
points(lhauteur.nt, rep(0,length(lhauteur.nt)), pch="I", cex=1.5,col="green", xlim=c(-2,6), ylim=c(0,0.4))
hist(lhauteur.nt,12, freq=FALSE, border="red",xlab="log pluviométrie",main="nuages non traités - 12 classes") 
points(lhauteur.nt, rep(0,length(lhauteur.nt)), pch="I", cex=1.5,col="green", xlim=c(-2,6), ylim=c(0,0.4))
par(op)

#estimations de densité
plot(density(lhauteur.t), col="blue",main="estimations de densité",xlab="log pluviométrie")
lines(density(lhauteur.nt), col="red")
legend("topleft",col=c("red","blue"),lty=1,legend=c("non traités","traités"))

op <- par(mfrow=c(1, 2))
plot(density(lhauteur.t,bw=0.2), col="blue",main="estimations de densité pour h=0.2",xlab="log pluviométrie")
lines(density(hauteur.nt,bw=0.2), col="red")
legend("topleft",col=c("red","blue"),lty=1,legend=c("non traités","traités"))
plot(density(lhauteur.t,bw=0.6), col="blue",main="estimations de densité pour h=0.6",xlab="log pluviométrie")
lines(density(hauteur.nt,bw=0.6), col="red")
legend("topleft",col=c("red","blue"),lty=1,legend=c("non traités","traités"))
par(op)

op <- par(mfrow=c(1, 2))
plot(density(lhauteur.t,bw=1), col="blue",main="estimations de densité pour h=1",xlab="log pluviométrie")
lines(density(hauteur.nt,bw=1), col="red")
legend("topleft",col=c("red","blue"),lty=1,legend=c("non traités","traités"))
plot(density(lhauteur.t,bw=4), col="blue",main="estimations de densité pour h=4",xlab="log pluviométrie")
lines(density(hauteur.nt,bw=4), col="red")
legend("topleft",col=c("red","blue"),lty=1,legend=c("non traités","traités"))
par(op)