mon master2 ISIFAR

ISIFAR
 
AccueilFAQRechercherS'enregistrerMembresGroupesConnexion

Partagez | 
 

 TP

Voir le sujet précédent Voir le sujet suivant Aller en bas 
AuteurMessage
Admin
Admin


Nombre de messages : 418
Date d'inscription : 27/09/2005

MessageSujet: TP   Ven 3 Fév à 0:44

Code:
rm(list=ls())
#dev.off()

#################################################################
#                                                              #
#                          INTROSPECTION                   #
#                                                              #
#################################################################

Xd=read.table('dataoneyear.txt' ,header=TRUE, sep=";");
X=Xd[,-1];

aaa=read.table('datatest05.txt' ,header=TRUE, sep=";");
bbb=read.table('datatest10.txt' ,header=TRUE, sep=";");
test1=aaa[,-1];
test2=bbb[,-1];


#visualisation des series temporelles
#donnees synchrones, agregees a l'heure
#nom=names(X)
#par(mfrow=c(2,3))
#for (i in 1:ncol(X))
#{ plot(X[,i],xlab="Temps",ylab=nom[i]) }

#ouverture d'une nouvelle fenetre et visualisation des  histogrammes
#x11()
#par(mfrow=c(2,3))
#for (i in 1:ncol(X))
#{ hist(X[,i],xlab=nom[i], main=paste("histogramme de",nom[i])) }

#for (i in 1:ncol(X))
# { summary(X[,i]) }
summary(X[,1]) #=> Y
summary(X[,2]) #=> Q
summary(X[,3]) #=> P1
summary(X[,4]) #=> P2
summary(X[,5]) #=> T1
summary(X[,6]) #=> T2



#################################################################
#                                                              #
#                          PRE-TRAITEMENT                   #
#                                                              #
#################################################################

########################
#je parcours chaque element de ma matrice et je supprime la ligne  fausse
#cree la matrice de taille (nbLigX x nbColX) remplie de 0
########################
U=matrix(rep(0,ncol(X)*nrow(X)),nrow(X),ncol(X));
#en parcourant chaque vecteur Colonne
for (i in 1:ncol(X))
{
  U[,i]=!is.finite(X[,i])
}
ind=(apply(U,1,sum)<1)
Xf=X[ind,];

#########################
#Pourcentage de donnees restantes apres enlevement des lignes de
#donnees manquantes
#########################
dim(X)
dim(Xf)
dim(Xf)[1]/dim(X)[1]
#=>95%


########################
#ATTENTION, je n'ai pas supprime les points extremes
########################
#On supprime les points extremes (le premier centile : 1%)
T=matrix(rep(0,ncol(Xf)*nrow(Xf)),nrow(Xf),ncol(Xf));
for (i in 1:ncol(Xf))
{
  T[,i]=!(Xf[,i]>quantile(Xf[,i], probs=0.01))
}
indx=(apply(T,1,sum)<1)
Xn=Xf[indx,];

#illustration
#x11()
#par(mfrow=c(2,3))
#for (i in 1:ncol(Xf))
#{ plot(Xn[,i],xlab="Temps",ylab=nom[i]) }


#########################
#Pourcentage de donnees restantes apres pre-traitement des  donnees
#########################
dim(X)
dim(Xn)
dim(Xn)[1]/dim(X)[1]
#=>91%




########################
#Correlations :
########################
########1ere idee(fausse)########
cor(Xf)
correlation1=round(cor(Xf),2)
#Y est correle avec Q et P2
#T1 correle avec T2

########2eme idee########
cor(Xn)
correlation2=round(cor(Xn),2)
#On constate que Y est seulement correle avec Q
#T1 correle avec T2


########################
#On verifie les correlations de toutes les variables avec Y en  plotant
#les donnees avant et après tous les pre-traitements
########################
#x11()
#par(mfrow=c(2,3))
#for (i in 1:(ncol(Xf)-1))
#{ plot(Xf[,i+1],Xf[,1],xlab=nom[i+1],ylab=nom[1],
#  main=(paste("correlation  :",correlation1[1,i+1])),col.main="red") }

#x11()
#par(mfrow=c(2,3))
#for (i in 1:(ncol(Xf)-1))
#{ plot(Xn[,i+1],Xn[,1],xlab=nom[i+1],ylab=nom[1],
#  main=(paste("correlation  :",correlation2[1,i+1])),col.main="red") }

#ouverture d'une nouvelle fenetre et visualisation des  histogrammes
#sur les donnees pre-traitees
#x11()
#par(mfrow=c(2,3))
#for (i in 1:ncol(X))
#{ hist(Xn[,i],xlab=nom[i], main=paste("histo. modif.  de",nom[i])) }


summary(Xn[,1]) #=> Y
summary(Xn[,2]) #=> Q
summary(Xn[,3]) #=> P1
summary(Xn[,4]) #=> P2
summary(Xn[,5]) #=> T1
summary(Xn[,6]) #=> T2


#################################################################
#                                                              #
#            "SEPARATION"  :  APPRENTISSAGE/TEST              #
#                                                              #
#################################################################
X_App=Xn

########################
# PRE-TRAITEMENTS fichiers de test, supprime "lignes NaN"
# pour garder un certain synchronisme
########################
U=matrix(rep(0,ncol(test1)*nrow(test1)),nrow(test1),ncol(test1));
#en parcourant chaque vecteur Colonne
for (i in 1:ncol(test1))
{
  U[,i]=!is.finite(test1[,i])
}
ind1=(apply(U,1,sum)<1)


V=matrix(rep(0,ncol(test2)*nrow(test2)),nrow(test2),ncol(test2));
#en parcourant chaque vecteur Colonne
for (j in 1:ncol(test2))
{
  V[,j]=!is.finite(test2[,j])
}
ind2=(apply(V,1,sum)<1)


ind=ind1                  ###########ind=ind2
X_Test=test1[ind,];      ###########X_Test=test2[ind,];
nrow(X_Test)/nrow(test1)  ###########nrow(X_Test)/nrow(test2)
#=>95%




#################################################################
#                                                              #
#                          MODELISATION                #
#                                                              #
#################################################################
yApp=matrix(X_App[,1])
yTest=matrix(X_Test[,1])
model=lm(yApp~Q+P1+P2+T1+T2,data=X_App)
modelTT=lm(yApp~Q+T1+T2,data=X_App)
modelPT=lm(yApp~Q+P2+T2,data=X_App)

#x11()
#par(mfrow=c(2,2))
#plot(model)
#summary(model)
#x11()
#par(mfrow=c(2,2))
#plot(modelTT)
#summary(modelTT)
#x11()
#par(mfrow=c(2,2))
#plot(modelPT)
#summary(modelPT)



#Estimation des valeurs de sorties pour la base de test
ypred=predict(model,X_Test)
ypredTT=predict(modelTT,X_Test)
ypredPT=predict(modelPT,X_Test)

#Estimation des valeurs de sorties pour la base d'APP
ypredApp=predict(model,X_App)
ypredAppTT=predict(modelTT,X_App)
ypredAppPT=predict(modelPT,X_App)




###Erreur quadratique moyenne###
Eq=(t(yTest-ypred)%*%(yTest-ypred))/nrow(X_Test)
sqrt(Eq); mean(yTest);sqrt(Eq)/ mean(yTest)
1-summary(model)$r.squared

#2e idee
Eq=(t(yTest-ypredTT)%*%(yTest-ypredTT))/nrow(X_Test)
sqrt(Eq); mean(yTest);sqrt(Eq)/ mean(yTest)
1-summary(modelTT)$r.squared

#3e idee
Eq=(t(yTest-ypredPT)%*%(yTest-ypredPT))/nrow(X_Test)
sqrt(Eq); mean(yTest);sqrt(Eq)/ mean(yTest)
1-summary(modelPT)$r.squared



#################################################################
#                                                              #
#                      ECART AU MODELE                        #
#              SURVEILLANCE DE L'EQUIPEMENT                    #
#                                                              #
#################################################################
#Surveillance horaire de l'equipement
Z=yTest-ypredTT
ZApp=yApp-ypredAppTT
summary(ZApp)
sd(ZApp)


#Affichage de la fonction de repartition des Z
x11()
par(mfrow=c(2,1))
xplot=sort(Z)
yplot=seq(0,1,by=1/(length(Z)-1))
plot(xplot,yplot,xlab="z(t)",ylab="P(Z<z(t))",
    main="Fonction de repartition empirique
          base de test"); grid()
xplot=sort(ZApp)
yplot=seq(0,1,by=1/(length(ZApp)-1))
plot(xplot,yplot,xlab="z(t)",ylab="P(ZApp<z(t))",
    main="Fonction de repartition empirique
          base d'apprentissage"); grid()


# Calcul des p-values Approche non paramétrique
compute_pvalue_ModEmp<-function(ent,ZApp)
{
pv=0*ent;
for (i in 1:length(ent))
   {
  pv[i]=(sum(ZApp>ent[i]))/length(ZApp);
   };
  return(pv)
}

# Calcul des p-values Approche paramétrique
compute_pvalue_ModGauss<-function(ent,ZApp)
{
  pvalue=pnorm(ent,mean(ZApp),sd(ZApp))
  return(pvalue)
}


#Comparaison des fonctions de repartition empirique et parametrique
x11()
xplot=sort(Z)
yplot=seq(0,1,by=1/(length(Z)-1))
plot(xplot,yplot,xlab="z(t)",ylab="P(Z<z(t))",col='blue',
    main="Comparaison des fonctions empirique
          et parametrique"); grid()
xplot=sort(Z)
yplot=pnorm(sort(Z),mean(ZApp),sd(ZApp))
points(xplot,yplot,col='red',pch=1)
lines(xplot,yplot,col='red',pch=1)
legend(-4000,0.7, " empirique  ", pch=21, pt.bg="white", lty=1, col = "blue")
legend(-4000,0.6, "parametrique", pch=21, pt.bg="white", lty=1, col = "red")


#################################################################
#                                                              #
#                          TEST DE                              #
#                  DYSFONCTIONNEMENT HORAIRE                  #
#                                                              #
#################################################################

#seuil
alpha=0.01
#p-value empirique
pve=compute_pvalue_ModEmp(Z,ZApp);
x11()
plot(pve,main='Affichage des p-value inférieures au seuil
              empirique',col='green');
ones=seq(1,length(pve));
ind=ones[pve<alpha];
indval=pve[ind];
abline(alpha,0,col='red');
points(ind,indval,col='red');



#p-value parametrique
pvg=compute_pvalue_ModGauss(Z,ZApp);
x11()
plot(pvg,main='Affichage des p-value inférieures au seuil
                parametrique (gaussien)',col='green');
ones=seq(1,length(pvg));
indG=ones[pvg<alpha];
indvalG=pvg[indG];
abline(alpha,0,col='red');
points(indG,indvalG,col='red');


#################################################################
#                                                              #
#                        TEST DE DERIVE                        #
#                      DE  DYSFONCTIONNEMENT                    #
#                                                              #
#################################################################
#Empirique
#24 heures
T=24
#Alarmes horaires Za
Za=Z
Za[ind]=1
Za[-ind]=0


#moyennes mobiles Za
mmza=rep(0,length(Za))
for (i in T : length(Za))
 {
  mmza[i]=mean(Za[(i-T+1):i,1])
 }
x11()
par(mfrow=c(2,1))
plot(Za,main='alarme horaire empirique'); grid();
plot(mmza,col='red',main='moyenne mobile (24H)'); grid();

#seuil pour les alarmes = 75%
seuil=0.75
indices=seq(1,length(mmza));
alerte=indices[mmza>seuil]
abline(seuil,0,col='blue')



#Parametrique (gaussien)
#24 heures
T=24
#Alarmes horaires Za
Za=Z
Za[indG]=1
Za[-indG]=0


#moyennes mobiles Za
mmzaG=rep(0,length(Za))
for (i in T : length(Za))
 {
  mmzaG[i]=mean(Za[(i-T+1):i,1])
 }
x11()
par(mfrow=c(2,1))
plot(Za,main='alarme horaire parametrique'); grid();
plot(mmzaG,col='red',main='moyenne mobile (24H)'); grid();

#seuil pour les alarmes = 75%
#si
seuil=0.75
indicesG=seq(1,length(mmza));
alerteG=indices[mmzaG>seuil]
abline(seuil,0,col='blue')

length(alerte)
length(alerteG)
Revenir en haut Aller en bas
Voir le profil de l'utilisateur http://mastertwo.jeun.fr
 
TP
Voir le sujet précédent Voir le sujet suivant Revenir en haut 
Page 1 sur 1

Permission de ce forum:Vous ne pouvez pas répondre aux sujets dans ce forum
mon master2 ISIFAR :: 2e semestre :: Consulting-
Sauter vers: