#+OPTIONS: ':nil *:t -:t ::t <:t H:3 \n:nil ^:nil arch:headline #+OPTIONS: author:t broken-links:nil c:nil creator:nil #+OPTIONS: d:(not "LOGBOOK") date:t e:t email:nil f:t inline:t num:t #+OPTIONS: p:nil pri:nil prop:nil stat:t tags:nil tasks:t tex:t #+OPTIONS: timestamp:t title:t toc:nil todo:nil |:t #+TITLE: Exercice 1, évaluation par les paires, module 3 #+DATE: <2018-09-24 lun.> #+AUTHOR: Christophe Pouzat #+EMAIL: christophe.pouzat@parisdescartes.fr #+LANGUAGE: fr #+SELECT_TAGS: export #+EXCLUDE_TAGS: noexport #+CREATOR: Emacs 25.3.1 (Org mode 9.0.9) #+LaTeX_CLASS: koma-article #+LaTeX_CLASS_OPTIONS: [koma,11pt] #+LaTeX_HEADER: \usepackage{cmbright} #+LaTeX_HEADER: \usepackage[round]{natbib} #+LaTeX_HEADER: \usepackage{alltt} #+LaTeX_HEADER: \usepackage[usenames,dvipsnames]{xcolor} #+LaTeX_HEADER: \renewenvironment{verbatim}{\begin{alltt} \scriptsize \color{Bittersweet} \vspace{0.2cm} }{\vspace{0.2cm} \end{alltt} \normalsize \color{black}} #+LaTeX_HEADER: \usepackage{listings} #+LaTeX_HEADER: \lstloadlanguages{C,Gnuplot,bash,sh,R} #+LaTeX_HEADER: \hypersetup{colorlinks=true,pagebackref=true} #+PROPERTY: header-args :eval no-export #+STARTUP: indent * Ajustement de variables d'environnement :noexport: Les variables nécessaires à la bonne exportation en LaTeX : #+NAME: org-latex-set-up #+BEGIN_SRC emacs-lisp :results silent :exports none (require 'ox-latex) (setq org-export-latex-listings t) (setq org-latex-listings 'listings) (setq org-latex-listings-options '(("frame" "lines") ("basicstyle" "\\footnotesize") ("numbers" "left") ("numberstyle" "\\tiny"))) (add-to-list 'org-latex-classes '("koma-article" "\\documentclass{scrartcl}" ("\\section{%s}" . "\\section*{%s}") ("\\subsection{%s}" . "\\subsection*{%s}") ("\\subsubsection{%s}" . "\\subsubsection*{%s}") ("\\paragraph{%s}" . "\\paragraph*{%s}") ("\\subparagraph{%s}" . "\\subparagraph*{%s}"))) (setq org-latex-pdf-process '("pdflatex -shell-escape -interaction nonstopmode -output-directory %o %f" "biber %b" "pdflatex -shell-escape -interaction nonstopmode -output-directory %o %f" "pdflatex -shell-escape -interaction nonstopmode -output-directory %o %f")) #+END_SRC Ajustement concernant =gnuplot= : #+NAME: ajuste-gnuplot-para #+BEGIN_SRC gnuplot :session *gnuplot* :results silent :eval no-export :exports none set terminal pngcairo size 1000,1000 #+END_SRC Redirection de la sortie =stderr= vers la sortie =stdout= afin que la première soit imprimée correctement (dans le rapport, j'ai trouvé comment faire ça sur le [[http://kitchingroup.cheme.cmu.edu/blog/2015/01/04/Redirecting-stderr-in-org-mode-shell-blocks/][Kitchin's blog]]) : #+NAME: stderr-redirection #+BEGIN_SRC emacs-lisp (setq org-babel-default-header-args:sh '((:prologue . "exec 2>&1") (:epilogue . ":")) ) (setq org-babel-use-quick-and-dirty-noweb-expansion t) #+END_SRC #+RESULTS: stderr-redirection : t * Téléchargement des données :export: Tout d'abord, je sauve la date et l'heure auxquelles cette analyse a été effectuée, avec le programme =date=, comme demandé dans l'énoncer : #+NAME: la-date-et-l-heure #+BEGIN_SRC sh :exports both :results output date #+END_SRC #+RESULTS: la-date-et-l-heure : mar. sept. 25 15:35:37 CEST 2018 Je récupère ensuite les données (hebdomadaires) à l'adresse indiquée sur la plateforme =FUN=, en utilisant la fonction =wget= : #+NAME: téléchargement-des-données-CO2-Mauna-Loa #+BEGIN_SRC sh :exports both :results output wget http://scrippsco2.ucsd.edu/assets/data/atmospheric\ /stations/in_situ_co2/weekly/weekly_in_situ_co2_mlo.csv #+END_SRC #+RESULTS: téléchargement-des-données-CO2-Mauna-Loa #+begin_example --2018-09-25 15:36:08-- http://scrippsco2.ucsd.edu/assets/data/atmospheric/stations/in_situ_co2/weekly/weekly_in_situ_co2_mlo.csv Résolution de scrippsco2.ucsd.edu (scrippsco2.ucsd.edu)… 169.228.224.138 Connexion à scrippsco2.ucsd.edu (scrippsco2.ucsd.edu)|169.228.224.138|:80… connecté. requête HTTP transmise, en attente de la réponse… 200 OK Taille : 62582 (61K) [text/csv] Sauvegarde en : « weekly_in_situ_co2_mlo.csv » 0K .......... .......... .......... .......... .......... 81% 138K 0s 50K .......... . 100% 870K=0,4s 2018-09-25 15:36:09 (162 KB/s) — « weekly_in_situ_co2_mlo.csv » sauvegardé [62582/62582] #+end_example Je regarde le contenu des 19 premières lignes du fichier avec : #+NAME: premières-19-lignes-du-fichier-de-données-CO2-Mauna-Loa #+BEGIN_SRC sh :exports both :results output head -n 19 weekly_in_situ_co2_mlo.csv #+END_SRC #+RESULTS: premières-19-lignes-du-fichier-de-données-CO2-Mauna-Loa #+begin_example "-------------------------------------------------------------------------------------------" " Atmospheric CO2 concentrations (ppm) derived from in situ air measurements " " at Mauna Loa, Observatory, Hawaii: Latitude 19.5°N Longitude 155.6°W Elevation 3397m " " " " Source: R. F. Keeling, S. J. Walker, S. C. Piper and A. F. Bollenbacher " " Scripps CO2 Program ( http://scrippsco2.ucsd.edu ) " " Scripps Institution of Oceanography (SIO) " " University of California " " La Jolla, California USA 92093-0244 " " " " Status of data and correspondence: " " " " These data are subject to revision based on recalibration of standard gases. Questions " " about the data should be directed to Dr. Ralph Keeling (rkeeling@ucsd.edu), Stephen Walker" " (sjwalker@ucsd.edu) and Stephen Piper (scpiper@ucsd.edu), Scripps CO2 Program. " " " " Baseline data in this file through 03-Sep-2018 from archive dated 04-Sep-2018 10:19:35 " " " "-------------------------------------------------------------------------------------------" #+end_example Je peux afficher les lignes allant de la 39e à la 44e -- afin d'avoir des informations sur l'organisation des données dans le fichier -- avec le programme =sed= (j'ai trouvé les « bonnes lignes » du fichier en explorant celui-ci semi-interactivement avec le [[https://fr.wikipedia.org/wiki/Visionneur_(Terminal)][visionneur]], ou /pager/, =less=) : #+NAME: lignes-39-à-44-du-fichier-de-données-CO2-Mauna-Loa #+BEGIN_SRC sh :exports both :results output sed -n '39,44p' < weekly_in_situ_co2_mlo.csv #+END_SRC #+RESULTS: lignes-39-à-44-du-fichier-de-données-CO2-Mauna-Loa : " " : " The data file below contains 2 columns indicaing the date and CO2 " : " concentrations in micro-mol CO2 per mole (ppm), reported on the 2008A " : " SIO manometric mole fraction scale. These weekly values have been " : " adjusted to 12:00 hours at middle day of each weekly period as " : " indicated by the date in the first column. " Les lignes 45 à 55 contiennent les 11 premières mesures : #+NAME: lignes-45-à-55-du-fichier-de-données-CO2-Mauna-Loa #+BEGIN_SRC sh :exports both :results output sed -n '45,55p' < weekly_in_situ_co2_mlo.csv #+END_SRC #+RESULTS: lignes-45-à-55-du-fichier-de-données-CO2-Mauna-Loa #+begin_example 1958-03-29, 316.19 1958-04-05, 317.31 1958-04-12, 317.69 1958-04-19, 317.58 1958-04-26, 316.48 1958-05-03, 316.95 1958-05-17, 317.56 1958-05-24, 317.99 1958-07-05, 315.85 1958-07-12, 315.85 1958-07-19, 315.46 #+end_example Enfin, la fonction =tail= affiche (par défaut) les dix dernières lignes du fichier : #+NAME: tail-du-fichier-de-données-CO2-Mauna-Loa #+BEGIN_SRC sh :exports both :results output tail weekly_in_situ_co2_mlo.csv #+END_SRC #+RESULTS: tail-du-fichier-de-données-CO2-Mauna-Loa #+begin_example 2018-06-23, 410.81 2018-06-30, 410.19 2018-07-07, 409.29 2018-07-14, 409.07 2018-07-21, 408.49 2018-07-28, 408.08 2018-08-04, 407.36 2018-08-11, 407.28 2018-08-18, 407.07 2018-08-25, 406.63 #+end_example ** Statistiques de bases et recherche de données aberrantes Je vais utiliser [[http://www.gnuplot.info/][gnuplot]] pour faire les figures et même un peu plus. Ici, il y a quelques « spécificités » du fichier de données dont je dois tenir compte : - les données proprement dites commencent en ligne 45, je vais donc utiliser =skip 44= lorsque la commande =plot= de =gnuplot= sera appelée ; - la virgule est utilisée pour séparer les colonnes du fichier, d'où le =set datafile separator ","= ci-dessous ; - l'abscisse correspond à une date ce que je précise par =set xdata time= ; - la date sera obtenue à partir de la première (année) et de la deuxième (mois) colonne, ce que je spécifie avec =set timefmt "%Y-%m-%d"= ; - je souhaite avoir les labels de l'abscisse sous forme "mois-année" ce que j'obtiens avec =set xtics format "%m-%Y"= ; - le reste devrait être assez clair. Avant d'utiliser concrètement =gnuplot=, j'inscris dans ce document la version utilisée afin d'améliorer la reproductibilité de cette analyse : #+NAME: version-de-gnuplot #+BEGIN_SRC sh :exports both :results output gnuplot --version #+END_SRC #+RESULTS: version-de-gnuplot : gnuplot 5.2 patchlevel 4 Ceci étant fait, je peux utiliser la fonction =stats= de =gnuplot= pour accéder à un résumé des données : #+NAME: co2_mlo-résumé #+BEGIN_SRC gnuplot :exports both :session *gnuplot* :eval no-export :results output set datafile separator "," stats "weekly_in_situ_co2_mlo.csv" skip 44 using 1:2 #+END_SRC #+RESULTS: co2_mlo-résumé #+BEGIN_EXAMPLE FILE: Records: 3082 Out of range: 0 Invalid: 0 Column headers: 0 Blank: 0 Data Blocks: 1 COLUMNS: Mean: 1988.3303 354.0582 Std Dev: 17.2668 27.2089 Sample StdDev: 17.2696 27.2134 Skewness: -0.0109 0.3316 Kurtosis: 1.8177 1.9420 Avg Dev: 14.9273 23.3342 Sum: 6.12803e+06 1.09121e+06 Sum Sq.: 1.21855e+10 3.88633e+08 Mean Err.: 0.3110 0.4901 Std Dev Err.: 0.2199 0.3466 Skewness Err.: 0.0441 0.0441 Kurtosis Err.: 0.0882 0.0882 Minimum: 1958.0000 [ 0] 313.0400 [ 60] Maximum: 2018.0000 [3081] 411.8100 [3068] Quartile: 1974.0000 329.4100 Median: 1988.0000 351.6000 Quartile: 2003.0000 375.7800 Linear Model: y = 1.559 x - 2745 Slope: 1.559 +- 0.004175 Intercept: -2745 +- 8.301 Correlation: r = 0.9891 Sum xy: 2.171e+09 #+END_EXAMPLE À ce stade il n'y a pas de signe de présence de données aberrantes. * Graphe global :export: La figure montrant la suite chronologique de la concentration de $CO_2$ en fonction de la date est obtenue avec : #+NAME: co2_mlo-fig #+HEADERS: :file co2_mlo.png #+BEGIN_SRC gnuplot :exports both :session *gnuplot* :eval no-export set datafile separator "," set xdata time set timefmt "%Y-%m-%d" set xtics format "%m-%Y" rotate set grid set xrange [*:*] set yrange [300:425] set title '[CO_{2}] atmosphérique (observatoire de Mauna Loa, Hawaï)' set xlabel 'Date' set ylabel '[CO_{2}] (ppm)' unset key plot "weekly_in_situ_co2_mlo.csv" skip 44 \ using 1:2 with lines linecolor 'red' #+END_SRC #+RESULTS: co2_mlo-fig [[file:co2_mlo.png]] * Estimer la tendance :export: ** Ajustement d'un modèle paramétrique Pour l'estimation de la tendance, je dois d'une manière ou d'une autre « supprimer » les oscillations périodiques annuelles dues à la pousse des feuilles des feuillus (phase décroissante) puis à leur chute et décomposition (phase croissante), comme expliqué dans la rubrique [[http://scrippsco2.ucsd.edu/history_legacy/early_keeling_curve][History & Legacy]] du site sur lequel se trouvent les données (c'est ce qui se passe dans l'hémisphère nord qui domine du fait de la plus grande surface de forêts à feuilles caducs qu'on y trouve). Une façon de « supprimer » des oscillations périodiques à fréquence relativement élevée (nous avons des oscillations dont la période est l'année alors que nous disposons de 60 ans d'observations) est d'intégrer le signal (l'intégrale d'une fonction sinus ou cosinus sur une période est nulle). C'est ce que je fais maintenant en prenant comme unité de temps l'année (52 semaines) et en soustrayant la « ligne de base » égale à 313 (d'après les statistiques obtenues ci-dessus) : #+NAME: co2_mlo-cumsum-fig-1 #+HEADERS: :file co2_mlo_cumsum_1.png #+BEGIN_SRC gnuplot :exports both :session *gnuplot* :eval no-export set datafile separator "," set grid set title '[CO_{2}] sommée (observatoire de Mauna Loa, Hawaï)' set xlabel 'Année' set ylabel '[CO_{2}] intégrée (ppm x temps)' unset key valeur_initiale = 313 s=0 plot [] [0:120000] "weekly_in_situ_co2_mlo.csv" skip 44 \ using ($0/52):(s=s+($2-valeur_initiale)) \ with lines lc 'red' lw 3 #+END_SRC #+RESULTS: co2_mlo-cumsum-fig-1 [[file:co2_mlo_cumsum_1.png]] Sur ces données intégrées, les oscillations ont disparu et la courbe a un aspect bien plus lisse que sur les données de départ. Je vais à présent utiliser la commande =fit= de =gnuplot= pour ajuster par la méthode des moindres carrés un polynôme de degré 4 aux données intégrées. J'ai choisi le degré du polynôme est les valeurs initiales des paramètres ci-dessous en répétant la procédure suivante en mode interactif: #+NAME: co2_mlo-cumsum-fit #+BEGIN_SRC gnuplot :exports both :session *gnuplot* :eval no-export set datafile separator "," s = 0 h(x) = a*x + b*x**2 +c*x**3+d*x**4 a=1 b=12.5 c=0.32 d=0.0005 fit h(x) "weekly_in_situ_co2_mlo.csv" skip 44 \ using ($0/52):(s=s+($2-valeur_initiale)) via a,b,c,d #+END_SRC #+RESULTS: co2_mlo-cumsum-fit #+BEGIN_EXAMPLE iter chisq delta/lim lambda a b c d 0 1.9126393665e+11 0.00e+00 1.60e+04 1.000000e+00 1.250000e+01 3.200000e-01 5.000000e-04 set output 1 1.8872569509e+08 -1.01e+08 1.60e+03 1.014841e+00 2.729204e+01 1.160925e-01 4.913335e-04 2 8.3513646197e+07 -1.26e+05 1.60e+02 1.341773e+00 2.943099e+01 1.887339e-02 1.551113e-03 3 6.3051460268e+07 -3.25e+04 1.60e+01 1.891717e+01 2.818995e+01 4.599273e-02 1.364171e-03 4 4.0624049973e+07 -5.52e+04 1.60e+00 6.742819e+01 2.389166e+01 1.620568e-01 3.848147e-04 5 4.0606923156e+07 -4.22e+01 1.60e-01 6.880592e+01 2.376957e+01 1.653537e-01 3.569931e-04 6 4.0606923154e+07 -3.40e-06 1.60e-02 6.880632e+01 2.376954e+01 1.653547e-01 3.569852e-04 iter chisq delta/lim lambda a b c d After 6 iterations the fit converged. final sum of squares of residuals : 4.06069e+07 rel. change during last iteration : -3.40215e-11 degrees of freedom (FIT_NDF) : 3078 rms of residuals (FIT_STDFIT) = sqrt(WSSR/ndf) : 114.859 variance of residuals (reduced chisquare) = WSSR/ndf : 13192.6 Final set of parameters Asymptotic Standard Error ======================= ========================== a = 68.8063 +/- 1.21 (1.758%) b = 23.7695 +/- 0.1107 (0.4656%) c = 0.165355 +/- 0.003157 (1.909%) d = 0.000356985 +/- 2.819e-05 (7.896%) correlation matrix of the fit parameters: a b c d a 1.000 b -0.968 1.000 c 0.917 -0.986 1.000 d -0.866 0.958 -0.992 1.000 #+END_EXAMPLE Je peux à présent faire un graphe avec la concentration intégrée (en rouge épais) et le modèle ajusté (noir fin) : #+NAME: co2_mlo-cumsum-fig #+HEADERS: :file co2_mlo_cumsum.png #+BEGIN_SRC gnuplot :exports both :session *gnuplot* :eval no-export set datafile separator "," set grid set title '[CO_{2}] intégrée et modèle polynomial ajusté' set xlabel 'Année' set ylabel '[CO_{2}] intégrée (ppm x temps)' set key s = 0 plot [] [0:120000] "weekly_in_situ_co2_mlo.csv" skip 44 \ using ($0/52):(s=s+($2-valeur_initiale)) \ with lines linecolor 'red' linewidth 2 title '[CO_{2}] intégrée', \ h(x) linecolor 'black' linewidth 1 dashtype '-' title 'Polynôme de degré 4' #+END_SRC #+RESULTS: co2_mlo-cumsum-fig [[file:co2_mlo_cumsum.png]] ** Estimation d'un modèle non-paramétrique Une autre façon d'estimer la tendance consiste à travailler année par année et, en utilisant le fait que la croissance / dégénérescence des feuilles caducs -- responsable des oscillations annuelles -- est de somme nulle sur l'année, à intégrer sur chaque année les concentrations. On peut alors approcher la vraie tendance par une fonction constante par morceaux -- la constante de chaque année étant la moyenne calculée sur cette même année -- ou par une fonction linéaire par morceaux -- la pente étant constante sur l'année. Pour faire ces calculs, j'ai besoin d'outils un petit peu plus sophistiqués que ceux mis à ma disposition par =gnuplot= et je vais employer =Python=. *** Lecture des données dans =Python= Je lis les données ligne pas ligne en faisant attention au fait que les 44 premières lignes du fichier contiennent des méta-données. Je crée trois [[https://en.wikipedia.org/wiki/Python_syntax_and_semantics#Lists,_tuples,_sets,_dictionaries][listes]] : - la_date :: contient la date sous forme de chaîne caractères, telle qu'elle se trouve dans le fichier ; - annee :: contient l'année (un entier) où chaque mesure a été effectuée (nécessaire pour calculer la moyenne annuelle) ; - CO2 :: contient la concentration (un réel) de CO$_{2}$. Pour « m'assurer » que l'importation s'est bien passée j'imprime les dix premières et les dix dernières valeurs de [CO$_{2}$] pour les comparer avec les valeurs précédentes (affichées avec =sed= et =tail=) : #+NAME: lecture-des-données-dans-Python #+BEGIN_SRC python :session *Mauna-Loa* :results output :exports both :eval no-export f = open("weekly_in_situ_co2_mlo.csv","r") titi = [f.readline() for i in range(44)] toto = f.readline() CO2 =[] annee = [] la_date = [] while toto: la_date.append(toto.split(',')[0]) CO2.append(float(toto.split(',')[1])) annee.append(int(toto.split('-')[0])) toto = f.readline() f.close() print("\nLes dix premières valeurs de [CO2] sont:\n") print(CO2[:10]) print("Les dix dernières valeurs de [CO2] sont:\n") print(CO2[-10:]) #+END_SRC #+RESULTS: lecture-des-données-dans-Python : : >>> >>> >>> >>> >>> ... ... ... ... ... >>> >>> : Les dix premières valeurs de [CO2] sont: : [316.19, 317.31, 317.69, 317.58, 316.48, 316.95, 317.56, 317.99, 315.85, 315.85] : Les dix dernières valeurs de [CO2] sont: : [410.81, 410.19, 409.29, 409.07, 408.49, 408.08, 407.36, 407.28, 407.07, 406.63] *** Calcul de la moyenne annuelle Je calcule maintenant la « moyenne », c'est-à-dire que pour chaque année je construis une liste à quatre éléments : - l'année ; - la première valeur de l'année ; - le nombre de semaines « dans l'année » (la première et la dernière années sont tronquées) ; - la moyenne sur l'année. Le résultat est stocké dans une liste de listes (=CO2_moyenne_annuelle=). Une fois cette liste créée, j'affiche les dix premières moyennes annuelles : #+NAME: moyenne-de-CO2-annuelle #+BEGIN_SRC python :session *Mauna-Loa* :results output :exports both :eval no-export CO2_moyenne_annuelle = list(range(1958,2019)) j = 0 for i in range(len(CO2_moyenne_annuelle)): cible = CO2_moyenne_annuelle[i] s = 0 n = 0 debut = CO2[j] while (j < len(annee)) and (annee[j] == cible): s += CO2[j] n += 1 j +=1 CO2_moyenne_annuelle[i] = [cible,debut,n,s/n] print("\nLes dix première moyennes annuelles sont:") print([round(l[3],1) for l in CO2_moyenne_annuelle[:10]]) #+END_SRC #+RESULTS: moyenne-de-CO2-annuelle : : >>> ... ... ... ... ... ... ... ... ... ... >>> : Les dix première moyennes annuelles sont: : [315.5, 315.9, 316.9, 317.6, 318.6, 319.0, 318.6, 320.0, 321.4, 322.2] Maintenant, j'écris dans les colonnes successives d'un fichier au format =csv= : - la date comme dans le fichier originale ; - la moyenne de la [CO$_2$] sur l'année en cours ; - la prédiction de la tendance de [CO$_2$] à partir de l'approximation linéaire ; - le résidu : [CO$_2$] - moyenne courante ; - le résidu : [CO$_2$] - approximation linéaire de la tendance. #+NAME: estimations-non-paramétriques-de-la-tendance #+BEGIN_SRC python :session *Mauna-Loa* :eval no-export f = open("co2_mlo_hebdo_non_para.csv","w") f.write("# date\t[CO2] moyenne\tTendance 'linéaire'\tRésidu (moyenne)\tRésidu (linéaire)\n") k=0 annee_courante = 1958 index = 0 for i in range(len(CO2)): if annee[i] > annee_courante: index = 0 annee_courante += 1 k += 1 moyenne_courante = CO2_moyenne_annuelle[k][3] debut = CO2_moyenne_annuelle[k][1] if annee_courante < 2018: fin = CO2_moyenne_annuelle[k+1][1] else: fin = CO2[-1] valeur_courante = debut + index*(fin-debut)/CO2_moyenne_annuelle[k][2] index += 1 tpl = (la_date[i], moyenne_courante, valeur_courante, CO2[i] - moyenne_courante, CO2[i] - valeur_courante) f.write("%s, %f, %f, %f, %f\n" % tpl) f.close() #+END_SRC #+BEGIN_SRC python :session *Mauna-Loa* :exports none :eval no-export f = open("co2_mlo_hebdo.csv","w") s = 0 a = 68.8063 b = 23.7695 c = 0.165355 d = 0.000356985 valeur_initiale = 313 k=0 annee_courante = 1958 for i in range(len(CO2)): if annee[i] > annee_courante: annee_courante += 1 k += 1 moyenne_courante = CO2_moyenne_annuelle[k] s += CO2[i]-valeur_initiale semaine = i/52. prevision = (a+((((4*d*semaine)+3*c)*semaine)+2*b)*semaine)/52. Prevision = (((((d*semaine)+c)*semaine)+b)*semaine+a)*semaine tpl = (la_date[i], semaine, CO2[i]-valeur_initiale, moyenne_courante, prevision, CO2[i]-valeur_initiale-prevision, CO2[i] - moyenne_courante, s, Prevision, s-Prevision) f.write("%s, %f, %f, %f, %f, %f, %f, %f, %f, %f\n" % tpl) f.close() #+END_SRC J'utilise =gnuplot= pour faire le graphe de la tendance basée sur l'approximation constante par morceau (moyenne annuelle) : #+NAME: co2_mlo-tendance-moyenne-fig #+HEADERS: :file co2_mlo_tendance-moyenne.png #+BEGIN_SRC gnuplot :exports both :session *gnuplot* :eval no-export set datafile separator "," set xdata time set timefmt "%Y-%m-%d" set xtics format "%m-%Y" rotate set grid set xlabel 'Date' set ylabel '[CO_{2}] moyenne sur une année (ppm)' unset key plot "co2_mlo_hebdo_non_para.csv" using 1:2 \ with lines lc 'red' lw 3 #+END_SRC #+RESULTS: co2_mlo-tendance-moyenne-fig [[file:co2_mlo_tendance-moyenne.png]] De même, je construis le graphe de la tendance basée sur l'approximation linéaire par morceaux (les estimations des deux années extrêmes ne sont pas fiables puisque les premiers mois de l'année 1958 manquent, de même que les derniers mois de l'année 2018) : #+NAME: co2_mlo-tendance-lineaire-fig #+HEADERS: :file co2_mlo_tendance-lineaire.png #+BEGIN_SRC gnuplot :exports both :session *gnuplot* :eval no-export set datafile separator "," set xdata time set timefmt "%Y-%m-%d" set xtics format "%m-%Y" rotate set grid set xlabel 'Date' set ylabel '[CO_{2}] approx. linéaire de tendance (ppm)' unset key plot "co2_mlo_hebdo_non_para.csv" using 1:3 \ with lines lc 'red' lw 3 #+END_SRC #+RESULTS: co2_mlo-tendance-lineaire-fig [[file:co2_mlo_tendance-lineaire.png]] Si nous comparons ces données avec la première figure de la page Wikipédia [[https://fr.wikipedia.org/wiki/Ressources_et_consommation_%C3%A9nerg%C3%A9tiques_mondiales][Ressources et consommation énergétiques mondiales]], tirée du [[https://www.bp.com/en/global/corporate/energy-economics/statistical-review-of-world-energy.html][BP Statistical Review of World Energy]], nous constatons que la relation entre [CO$_2$] (mesurée à Mauna Loa) et consommation d'énergie fossile au niveau mondiale est loin d'être simple (en principe il faudrait intégrer la consommation, si nous supposons que tout le carbone rejeté s'accumule dans l'atmosphère). Les doubles chocs pétroliers des années 70 (à la suite desquels la pente de la consommation de pétrole diminue nettement) ne sont suivis d'aucun effet sur la [CO$_2$] -- on devine à peine un léger coude après 1973. La forte accélération de la consommation de charbon après 2000 ne se voit pas non plus sur la [CO$_2$]. Il faut clairement prendre en compte le rôle de tampon / réserve de CO$_2$ joué par les océans (qui deviennent plus [[https://fr.wikipedia.org/wiki/Acidification_des_oc%C3%A9ans][acides]]), mais là, ça se complique... Enfin, « l'inconvénient » de ces méthodes non-paramétriques locales est qu'elles ne permettent pas de faire d'extrapolation (sensée). * Estimer les oscillations :export: ** Avec l'approximation constante par morceaux de la tendance Je construis le graphe des résidus de [CO$_2$] obtenu en soustrayant l'approximation constante par morceaux de la tendance : #+NAME: co2_mlo-oscillation-moyenne-fig #+HEADERS: :file co2_mlo_oscillation-moyenne.png #+BEGIN_SRC gnuplot :exports both :session *gnuplot* :eval no-export set datafile separator "," set xdata time set timefmt "%Y-%m-%d" set xtics format "%m-%Y" rotate set grid set title 'Variation annuelle de [CO_{2}] (tendance approx. par fct cst par morceaux)' set xlabel 'Date' set ylabel '[CO_{2}] (ppm)' unset key plot "co2_mlo_hebdo_non_para.csv" using 1:4 \ with lines lc 'red' lw 2 #+END_SRC #+RESULTS: co2_mlo-oscillation-moyenne-fig [[file:co2_mlo_oscillation-moyenne.png]] ** Avec l'approximation linéaire par morceaux de la tendance Je construis le graphe des résidus de [CO$_2$] obtenu en soustrayant l'approximation linéaire par morceaux de la tendance : #+NAME: co2_mlo-oscillation-lineaire-fig #+HEADERS: :file co2_mlo_oscillation-lineaire.png #+BEGIN_SRC gnuplot :exports both :session *gnuplot* :eval no-export set datafile separator "," set xdata time set timefmt "%Y-%m-%d" set xtics format "%m-%Y" rotate set grid set title 'Variation annuelle de [CO_{2}] (tendance approx. par fct linéaire par morceaux)' set xlabel 'Date' set ylabel '[CO_{2}] (ppm)' unset key plot "co2_mlo_hebdo_non_para.csv" using 1:5 \ with lines lc 'red' lw 2 #+END_SRC #+RESULTS: co2_mlo-oscillation-lineaire-fig [[file:co2_mlo_oscillation-lineaire.png]] Les deux graphes (surtout le second) suggèrent une augmentation de l'amplitude des oscillations.