X-Git-Url: https://git.auder.net/?p=talweg.git;a=blobdiff_plain;f=Etude_En_Cours_R_E_bis.R;h=aa0d4900e02ef048c7952ab4bcd878721402ce54;hp=1346125b3418463617d5c14643743cada1345c16;hb=794d99e9614315f8a280ac70a991619d4cf45f87;hpb=9e0f25f6d5c73008c0ad59ea7b7e097b342f26b3 diff --git a/Etude_En_Cours_R_E_bis.R b/Etude_En_Cours_R_E_bis.R index 1346125..aa0d490 100644 --- a/Etude_En_Cours_R_E_bis.R +++ b/Etude_En_Cours_R_E_bis.R @@ -2,25 +2,25 @@ rm(list=ls()) # Lecture des données: pm = dataframe 2 colonnes, date-time puis PM10 horaire -pm = read.table("DATA/mesures_horaires_hloc_pm10_a_filer.csv",sep=";",dec=".",header=T) +pm = read.table("DATA/mesures_horaires_hloc_pm10_a_filer.csv",sep=",",dec=".",header=T) #n = dim(pm)[1] #datedebut = "10/12/2008" # Chargement des données météo et indicateurs: VarExp matrice des données météo, # première colonne = date, première rangée = second jour -VarExp <- read.table("DATA/meteo_extra_jourMois.csv",sep=";",dec=".",header=T) +VarExp <- read.table("DATA/meteo_extra_jourMois.csv",sep=",",dec=".",header=T) VarExp <- VarExp[-1,] # Lecture des dates: ??? sans doute liste des jours de 10/12/2008 à 06/2016... -dates = read.table("DATA/Dates_jours.csv",header=F,as.is=T) -dates = dates[,1] +#dates = read.table("DATA/Dates_jours.csv",header=F,as.is=T) +dates = VarExp[,1] #dates[,1] # Contruction des matrices de données: pm.h = matrice des séries en ligne pm.h <- matrix(pm[,2],ncol=24,byrow=TRUE) Nlignes = nrow(pm.h) # Data = matrice des couples de jours (séries de 48 PM10), -# première ligne = jour 1/2, dernière = jour N-1/N +# première ligne = jour 1,2, dernière = jour N-1,N Data = cbind(pm.h[1:(Nlignes -1), ], pm.h[2:Nlignes, ]) # dates2 = dates du 2eme jour au dernier @@ -35,10 +35,10 @@ df <- Data PMjour <- apply(df[,25:48],1,mean,na.rm=T) dfexp <- cbind(VarExp,PMjour) -Dates = c( -"16/03/2015", -"19/01/2015", -"27/04/2015") +Dates = c( #Month/Day/Year, as in meteo file +"3/16/2015", +"1/19/2015", +"4/27/2015") Categorie = c("Epandage", "Chauffage", "Non Polluée") @@ -54,23 +54,24 @@ Err24 = NULL ErrPrev = NULL Kvois = NULL -for (Hc in 5:24) -{ +#for (Hc in 5:24) +#{ + Hc = 7 #predict at 7:00, from 8:00 H=24+Hc #Hc: dernier PM10 connu avant prédiction des 24-Hc restants, le jour 2 L = 1:H # Premier conditionnement : mois - indcond <- dfexp[,"Mois"] == 2 | dfexp[,"Mois"] == 3 | dfexp[,"Mois"] == 4 - | dfexp[,"Mois"] == 9 | dfexp[,"Mois"] == 10 + indcond <- dfexp[,"Mois"] == 2 | dfexp[,"Mois"] == 3 | dfexp[,"Mois"] == 4 | dfexp[,"Mois"] == 9 | dfexp[,"Mois"] == 10 data = df[indcond,] #restriction aux couples dont le mois du 2eme jour est 2,3,4,9 ou 10 varexp = dfexp[indcond,] #de même sur les exogènes (+ PM10 moyens) nl = (1:nrow(data))[rownames(data)==Dates[j]] #numéro de ligne où 2eme jour == 16/03/2015 - dateJPrev = rownames(data)[nl+ij] # ??? la date du 1er jour du couple de jours pile une semaine après le couple 'nl' + dateJPrev = rownames(data)[nl+ij] # translation : date du jour à prévoir (2eme colonne) dataj = as.numeric(data[nl+ij, 1:48]) #extraction de cette série dans dataj - data = data[-(nl + ij), ] #suppression de cette série dans data (????) - varexp = varexp[-(nl + ij), ] #idem dans les variables exogènes (date 16/03/2015 + 6 jours + data = data[-(nl + ij), ] #suppression de cette série dans data + varexp = varexp[-(nl + ij), ] #idem dans les variables exogènes indNA = attr(na.omit(data[, 1:48]),"na.action") +# print(indNA) data = data[-indNA,] #remove all pairs of days with NAs varexp = varexp[-indNA,] #and remove corresponding exogenous variables @@ -88,62 +89,63 @@ for (Hc in 5:24) D[k] = sqrt(sum((dataj[L] - data[k,L])^2)) } ind = order(D)[1:nbvois] - w = 1/(D[ind]^2) - w = w/sum(w) - W = w %o% rep(1,48) +# w = 1/(D[ind]^2) +# w = w/sum(w) +# W = w %o% rep(1,48) JourMoy = apply(data[ind, 1:48], 2, mean) #JourMoy = apply(W*data[ind, 1:48], 2, sum) NomFile = paste("Voisins_Epandage_PMjour_Hc_",Hc,".png",sep="") Titre = paste("Jour à prévoir : ",dateJPrev," - ", length(ind)," voisins",sep="") #erreur = sqrt(sum((dataj[25:48] - JourMoy[25:48])^2)) - if(Hc==24) + if(Hc==24) { erreurPrev = NA - else + } else { erreurPrev = mean(abs(dataj[(H+1):48] - JourMoy[(H+1):48])) - erreur24 = mean(abs(dataj[25:48] - JourMoy[25:48])) - #png(NomFile) - matplot(t(data[ind, 1:48]), type = "l", lwd=1.4, lty=1, col=1:length(ind), - cex.axis=1.4, cex.main = 1.7, cex.lab=1.5, - xlab="Heures locales", ylab=paste0("PM10 - Erreurs = ",round(erreur24,1), - " / ",round(erreurPrev,1)), main=Titre) - legend("top",rownames(data)[ind],ncol=2,lwd=1.4, lty=1, col=1:length(ind)) - lines(1:48, dataj, lwd=2.5) - lines(JourMoy, lty = 2, lwd=2) - abline(v=c(24.5,H+0.5), lty = 2, lwd=1.2) - #xx=dev.off() - - Err24 = c(Err24, erreur24) - ErrPrev = c(ErrPrev, erreurPrev) - ResDates = cbind(ResDates, rownames(data)[ind]) -} - -rownames(ResDates) = 1:10 - -Kvois = NULL -for (Col in ncol(ResDates):1) -{ - K = 0 - for (I in 1:10) - { - for (J in 1:10) - { - if (ResDates[I,1] == ResDates[J,Col]) - K = K +1 - } } - Kvois = c(Kvois, K) -} - -# pdf("Erreur_Epandage_PMjour.pdf") -ymin = min(na.omit(ErrPrev), Err24) -ymin = min(ymin, Kvois) -ymax = max(na.omit(ErrPrev), Err24) -ymax = max(ymax, Kvois) -plot(5:24,Err24, type = "l", lwd = 2, cex.axis=1.4, cex.lab = 1.5, - ylab = "MAE",xlab = "Heures de prévision", ylim= c(ymin,ymax)) -lines(5:24,ErrPrev, lwd=2, lty = 2) -legend("topright", legend=c("Erreur 24h", "Erreur prévision"), lty = c(1,2), lwd=2, cex=1.5) -points(5:24, Kvois, cex=1.8, pch = 19) -# xx = dev.off() - -#length(D) +# erreur24 = mean(abs(dataj[25:48] - JourMoy[25:48])) +# #png(NomFile) +# matplot(t(data[ind, 1:48]), type = "l", lwd=1.4, lty=1, col=1:length(ind), +# cex.axis=1.4, cex.main = 1.7, cex.lab=1.5, +# xlab="Heures locales", ylab=paste0("PM10 - Erreurs = ",round(erreur24,1), +# " / ",round(erreurPrev,1)), main=Titre) +# legend("top",rownames(data)[ind],ncol=2,lwd=1.4, lty=1, col=1:length(ind)) +# lines(1:48, dataj, lwd=2.5) +# lines(JourMoy, lty = 2, lwd=2) +# abline(v=c(24.5,H+0.5), lty = 2, lwd=1.2) +# #xx=dev.off() +# +# Err24 = c(Err24, erreur24) +# ErrPrev = c(ErrPrev, erreurPrev) +# ResDates = cbind(ResDates, rownames(data)[ind]) +#} + +#rownames(ResDates) = 1:10 +# +#Kvois = NULL +#for (Col in ncol(ResDates):1) +#{ +# K = 0 +# for (I in 1:10) +# { +# for (J in 1:10) +# { +# if (ResDates[I,1] == ResDates[J,Col]) +# K = K +1 +# } +# } +# Kvois = c(Kvois, K) +#} +# +## pdf("Erreur_Epandage_PMjour.pdf") +#ymin = min(na.omit(ErrPrev), Err24) +#ymin = min(ymin, Kvois) +#ymax = max(na.omit(ErrPrev), Err24) +#ymax = max(ymax, Kvois) +#plot(5:24,Err24, type = "l", lwd = 2, cex.axis=1.4, cex.lab = 1.5, +# ylab = "MAE",xlab = "Heures de prévision", ylim= c(ymin,ymax)) +#lines(5:24,ErrPrev, lwd=2, lty = 2) +#legend("topright", legend=c("Erreur 24h", "Erreur prévision"), lty = c(1,2), lwd=2, cex=1.5) +#points(5:24, Kvois, cex=1.8, pch = 19) +## xx = dev.off() +# +##length(D)