X-Git-Url: https://git.auder.net/?p=talweg.git;a=blobdiff_plain;f=Etude_En_Cours_R_E_bis.R;fp=Etude_En_Cours_R_E_bis.R;h=0000000000000000000000000000000000000000;hp=56240d07e421f4b7a73f264ce2d4587fe47d4c1b;hb=096e9798a4241f14d12c3663f8035d8def43d7e3;hpb=f71b975b140342f7ea80275359b7f4f9aa75153a

diff --git a/Etude_En_Cours_R_E_bis.R b/Etude_En_Cours_R_E_bis.R
deleted file mode 100644
index 56240d0..0000000
--- a/Etude_En_Cours_R_E_bis.R
+++ /dev/null
@@ -1,175 +0,0 @@
-predi <- function(ij)
-{
-#setwd("/Users/bp/Desktop/CONTRATS_AirNormand/2016/RapportFinalBruno")
-#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)
-#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 <- 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 = 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
-Data = cbind(pm.h[1:(Nlignes -1), ], pm.h[2:Nlignes, ])
-
-# dates2 = dates du 2eme jour au dernier
-dates2 = dates[2:Nlignes]
-
-rownames(Data) = dates2
-# df contient l'ensemble des données.
-#df <- cbind(Data,varexp[,-1])
-df <- Data
-
-# Complétion des variables exogènes (du 2eme jour) par les PM10 moyen ce même jour
-PMjour <- apply(df[,25:48],1,mean,na.rm=T)
-dfexp <- cbind(VarExp,PMjour)
-
-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")
-
-RepFig = "FIGURES_Etude"
-
-ResDates = NULL
-
-nbvois=10
-j=1 # numéro de semaine
-#ij=6 # numéro du jour (0 = lundi)
-
-Err24 = NULL
-ErrPrev = NULL
-Kvois = NULL
-
-#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
-	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] # 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
-	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
-
-	# Second conditionnement : les jours avec PMjour +/- large
-	large = 1
-	bornes = mean(dataj[25:48]) + c(-large,large)
-	indcond = varexp[,"PMjour"]>=bornes[1] & varexp[,"PMjour"]<=bornes[2]
-	if (sum(indcond) < 10)
-	{
-		large = 2
-		bornes = mean(dataj[25:48]) + c(-large,large)
-		indcond = varexp[,"PMjour"]>=bornes[1] & varexp[,"PMjour"]<=bornes[2]
-	}
-	while (sum(indcond) < 10)
-	{
-		large = large + 3
-		bornes = mean(dataj[25:48]) + c(-large,large)
-		indcond = varexp[,"PMjour"]>=bornes[1] & varexp[,"PMjour"]<=bornes[2]
-	}
-	data = data[indcond,] #pollution du 2eme jour == pollution du jour courant +/- 1
-	varexp = varexp[indcond,]
-
-	D = rep(0,nrow(data))
-	for (k in 1:nrow(data))
-	{
-		#D[k] = sqrt(sum((1:H)*(dataj[L] - data[k,L])^2))
-		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)
-
-#print("Voisins + ij")
-#print(sort(D)[1:nbvois])
-#print(dateJPrev)
-#print(rownames(data)[sort(ind)])
-#print(data[sort(ind), 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) {
-		erreurPrev = NA
-	} else {
-		erreurPrev = mean(abs(dataj[(H+1):48] - JourMoy[(H+1):48]))
-	}
-
-	list(serie=dataj, prev=JourMoy, err=erreurPrev, line=(nl+ij), neighbs=ind, dates=data[ind,1])
-#	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)
-}