Commit 4ee2051c authored by François Févotte's avatar François Févotte

Ex 03-3: résumé + correction taux d'accroissement

- passage à Plotly pour les figures
- correction du style
parent 7893ffdc
*~
\ No newline at end of file
Exercice.so
exercice.jl
\ No newline at end of file
...@@ -112,9 +112,9 @@ version = "2.0.0" ...@@ -112,9 +112,9 @@ version = "2.0.0"
[[DataStructures]] [[DataStructures]]
deps = ["InteractiveUtils", "OrderedCollections"] deps = ["InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "73eb18320fe3ba58790c8b8f6f89420f0a622773" git-tree-sha1 = "4dead20a1606a60292529023d6eac18a1ef6432e"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.17.11" version = "0.17.12"
[[DataValueInterfaces]] [[DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6" git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
...@@ -170,9 +170,9 @@ uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4" ...@@ -170,9 +170,9 @@ uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
version = "0.24.7" version = "0.24.7"
[[ElasticArrays]] [[ElasticArrays]]
git-tree-sha1 = "5b5b7cb8cba44bcf337b8af0a1f3e57c89468660" git-tree-sha1 = "6643de157ea3332d73e35a6a6ed2bdffe7792b12"
uuid = "fdbdab4c-e67f-52f5-8c3f-e7b388dad3d4" uuid = "fdbdab4c-e67f-52f5-8c3f-e7b388dad3d4"
version = "1.0.0" version = "1.1.0"
[[ElasticPDMats]] [[ElasticPDMats]]
deps = ["LinearAlgebra", "MacroTools", "PDMats", "Test"] deps = ["LinearAlgebra", "MacroTools", "PDMats", "Test"]
......
#!/bin/bash
#=
exec julia --color=yes --startup-file=no "${BASH_SOURCE[0]}" "$@"
=#
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()
using PackageCompiler
@time create_sysimage(Symbol.(keys(Pkg.project().dependencies));
sysimage_path=Pkg.project().name * ".so",
precompile_execution_file=joinpath(@__DIR__, "make.jl"))
# Local Variables:
# mode: julia
# End:
This source diff could not be displayed because it is too large. You can view the blob instead.
...@@ -4,8 +4,61 @@ author : François Févotte ...@@ -4,8 +4,61 @@ author : François Févotte
date: avril 2020 date: avril 2020
options: options:
css: skeleton_css.css css: skeleton_css.css
template: julia_html.tpl
--- ---
Ce document a été réalisé à l'occasion du [MOOC "Recherche
Reproductible"](https://www.fun-mooc.fr/courses/course-v1:inria+41016+self-paced/info). Il
a donc plus vocation à fournir des résultats *reproductibles* que des résultats
*intéressants*. N'étant pas un scientifique des données, il est tout à fait
possible que l'analyse réalisée ici soit caduque ; je plaide pour l'indulgence
du lecteur !
Ce document a été automatiquement généré à partir du fichier
[`exercice.jmd`](exercice.jmd), que le lecteur intéressé pourra consulter
directement. Une [annexe](repro.html) détaille toutes les étapes nécessaires
afin de reproduire l'analyse et re-générer intégralement le présent
document. Cette annexe, nommée [`repro.html`](repro.html), devrait se trouver
dans le même dépôt git et même répertoire que le présent document.
------
**Résumé:**
Depuis 1958, la concentration atmosphérique en CO2 est mesurée régulièrement à
l'observatoire de Mauna Loa, à Hawaï. Ces mesures, initiées par [Charles
D. Keeling](https://en.wikipedia.org/wiki/Charles_David_Keeling), ont été les
parmi les premières à mettre en évidence l'accroissement rapide de la
concentration en CO2 dans l'atmosphère, un constat qui a par la suite permis
d'attirer l'attention de la communauté scientifique (puis par la suite du grand
public) sur l'impact des activités humaines sur l'atmosphère.
![Courbe de
Keeling](https://upload.wikimedia.org/wikipedia/commons/thumb/c/c5/Mauna_Loa_CO2_monthly_mean_concentration.svg/480px-Mauna_Loa_CO2_monthly_mean_concentration.svg.png)
```julia; echo=false;
monroe2015() = print("[[monroe2015](https://scripps.ucsd.edu/programs/keelingcurve/2015/02/12/is-the-rate-of-co2-growth-slowing-or-speeding-up/)]");
```
La [courbe de Keeling](https://en.wikipedia.org/wiki/Keeling_Curve) fait
apparaître une augmentation de la concentration atmosphérique de l'ordre de 1 à
2 ppm/an `j monroe2015()`. Elle fait aussi apparaître une variation
cyclique d'environ 5 ppm chaque année, qui rend l'analyse de la tendance
systématique haussière plus délicate.
Nous nous proposons ici de refaire pas à pas l'étude de ce jeu de données. Nous
séparons la composante cyclique de la composante tendancielle, afin de
caractériser cette dernière. Les paramètres de notre modèle correspondent à un
taux d'augmentation du CO2 atmosphérique de 0.8 ppm/an en 1959, et 1.5 ppm/an en
2015 (cette dernière valeur étant un peu sous-estimée par rapport aux 2.25
ppm/an annoncés dans `j monroe2015()`).
Notre modèle est construit en utilisant les données de 1959 à 2014. Il est
validé par comparaison aux données de 2015 à nos jours. Une extrapolation sur la
période 2020-2025 est aussi réalisée afin de permettre une comparaison à de
futures mesures.
# Gestion des dépendances # Gestion des dépendances
...@@ -36,6 +89,7 @@ analyse. ...@@ -36,6 +89,7 @@ analyse.
using Pkg using Pkg
Pkg.activate(@__DIR__) Pkg.activate(@__DIR__)
Pkg.instantiate() Pkg.instantiate()
Pkg.status()
``` ```
...@@ -56,7 +110,7 @@ using Printf ...@@ -56,7 +110,7 @@ using Printf
using Dates using Dates
using DataFrames using DataFrames
using Statistics using Statistics
using Plots; gr() using Plots; plotly()
``` ```
```julia; echo=false; results="hidden" ```julia; echo=false; results="hidden"
...@@ -133,10 +187,10 @@ end ...@@ -133,10 +187,10 @@ end
Le fichier est structuré en deux colonnes : Le fichier est structuré en deux colonnes :
- `date` : date de la mesure - `date` : date de la mesure
- `val` : concentration en CO2 (en ppm molaires) - `co2` : concentration en CO2 (en ppm molaires)
```julia; results="hidden" ```julia; results="hidden"
data_raw = CSV.read(data_file; skipto=skip+1, header=[:date, :val]) data_raw = CSV.read(data_file; skipto=skip+1, header=[:date, :co2])
``` ```
L'examen des premières et dernières lignes de données révèle qu'elles couvrent L'examen des premières et dernières lignes de données révèle qu'elles couvrent
...@@ -179,17 +233,19 @@ particulier ayant duré 19 semaines en 1964. Le traitement devra en tenir compte ...@@ -179,17 +233,19 @@ particulier ayant duré 19 semaines en 1964. Le traitement devra en tenir compte
Une visualisation de l'ensemble des données semble montrer une augmentation Une visualisation de l'ensemble des données semble montrer une augmentation
tendancielle de la concentration en CO2, à laquelle se superpose une oscillation tendancielle de la concentration en CO2, à laquelle se superpose une oscillation
à plus haute fréquence. à plus haute fréquence.
```julia; echo=false ```julia; echo=false; results="raw"
plot(xlabel="Date", ylabel="CO2 [ppm]", legend=:bottomright) plot(xlabel="Date", ylabel="CO2 [ppm]")
plot!(data_raw.date, data_raw.val, label="mesures") plot!(data_raw.date, data_raw.co2, label="mesures")
disp()
``` ```
Un zoom sur les dernières années laisse penser que le cycle court se produit sur Un zoom sur les dernières années laisse penser que le cycle court se produit sur
une période annuelle, avec un minimum local atteint chaque année autour du mois une période annuelle, avec un minimum local atteint chaque année autour du mois
d'octobre. d'octobre.
```julia; echo=false ```julia; echo=false; results="raw"
plot(xlabel="Date", ylabel="CO2 [ppm]", legend=:bottomright) plot(xlabel="Date", ylabel="CO2 [ppm]")
plot!(data_raw.date[end-150:end], data_raw.val[end-150:end], label="mesures") plot!(data_raw.date[end-150:end], data_raw.co2[end-150:end], label="mesures")
disp()
``` ```
# Analyse # Analyse
...@@ -277,7 +333,7 @@ de la "date numérique". C'est le paquet Julia ...@@ -277,7 +333,7 @@ de la "date numérique". C'est le paquet Julia
charge d'effectuer le gros du travail. charge d'effectuer le gros du travail.
```julia ```julia
interp = DI.LinearInterpolation(data_raw.val, date2num.(data_raw.date)); interp = DI.LinearInterpolation(data_raw.co2, date2num.(data_raw.date));
``` ```
Nous allons profiter de la construction des données interpolées pour gérer le Nous allons profiter de la construction des données interpolées pour gérer le
...@@ -286,13 +342,13 @@ problème des données manquantes : nous n'interpolerons aucune donnée dans les ...@@ -286,13 +342,13 @@ problème des données manquantes : nous n'interpolerons aucune donnée dans les
interpolées est aussi l'occasion d'enrichir les formats de représentation des interpolées est aussi l'occasion d'enrichir les formats de représentation des
dates. Nous avons maintenant 5 colonnes dans notre jeu de données interpolé : dates. Nous avons maintenant 5 colonnes dans notre jeu de données interpolé :
- `date` : date identifiant le jour de la mesure (ou de la valeur interpolée) - `date` : date identifiant le jour de la mesure (ou de la valeur interpolée)
- `val` : valeur de la mesure de CO2 (ou de l'interpolation) - `co2` : valeur de la mesure de CO2 (ou de l'interpolation)
- `date_num` : date convertie en nombre - `date_num` : date convertie en nombre
- `year` : composante de la date identifiant l'année - `year` : composante de la date identifiant l'année
- `day` : composante de la date identifiant le jour dans l'année (entre 0 et 365) - `day` : composante de la date identifiant le jour dans l'année (entre 0 et 365)
```julia; results="raw" ```julia; results="raw"
data_interp = DataFrame(date_num=Int[], val=Float64[], date=Date[], year=Int[], day=Int[]) data_interp = DataFrame(date_num=Int[], co2=Float64[], date=Date[], year=Int[], day=Int[])
for i in 2:length(dates) for i in 2:length(dates)
if dates[i]-dates[i-1] > 14Days if dates[i]-dates[i-1] > 14Days
# pas d'interpolation : seule la date de gauche est inclue # pas d'interpolation : seule la date de gauche est inclue
...@@ -312,7 +368,7 @@ for i in 2:length(dates) ...@@ -312,7 +368,7 @@ for i in 2:length(dates)
date = date, date = date,
year = year(date), year = year(date),
day = dayinyear(date), day = dayinyear(date),
val = interp(date_num))) co2 = interp(date_num)))
end end
end end
...@@ -324,16 +380,17 @@ manquantes, on observe bien le résultat attendu : une interpolation linéaire ...@@ -324,16 +380,17 @@ manquantes, on observe bien le résultat attendu : une interpolation linéaire
journalière lorsque les données sont disponibles, mais aucune interpolation journalière lorsque les données sont disponibles, mais aucune interpolation
lorsque les données sont manquantes. lorsque les données sont manquantes.
```julia; echo=false ```julia; echo=false; results="raw"
plot(xlabel="Date", ylabel="CO2 [ppm]", legend=:topleft) plot(xlabel="Date", ylabel="CO2 [ppm]")
idx1 = (data_interp.date .> Date("2006-01-20")) .& (data_interp.date .< Date("2006-03-15")) idx1 = (data_interp.date .> Date("2006-01-20")) .& (data_interp.date .< Date("2006-03-15"))
plot!(data_interp.date[idx1], data_interp.val[idx1], seriestype=:scatter, plot!(data_interp.date[idx1], data_interp.co2[idx1], seriestype=:scatter,
label="valeurs interpolées") label="valeurs interpolées")
idx2 = (data_raw.date .> Date("2006-01-20")) .& (data_raw.date .< Date("2006-03-15")) idx2 = (data_raw.date .> Date("2006-01-20")) .& (data_raw.date .< Date("2006-03-15"))
plot!(data_raw.date[idx2], data_raw.val[idx2], seriestype=:scatter, plot!(data_raw.date[idx2], data_raw.co2[idx2], seriestype=:scatter,
label="mesures") label="mesures")
disp()
``` ```
Enfin, nous ne gardons pour notre analyse que des années complètes. D'ailleurs, Enfin, nous ne gardons pour notre analyse que des années complètes. D'ailleurs,
...@@ -355,9 +412,10 @@ Sur ces années complètes, la composante `day` de la date devrait être ...@@ -355,9 +412,10 @@ Sur ces années complètes, la composante `day` de la date devrait être
équirépartie entre 0 et 365, ce qui est globalement le cas. Les données équirépartie entre 0 et 365, ce qui est globalement le cas. Les données
manquantes n'ont donc pas d'impact significatif de ce point de vue là. manquantes n'ont donc pas d'impact significatif de ce point de vue là.
```julia; echo=false ```julia; echo=false; results="raw"
plot(xlabel="Jour de l'année", ylabel="Nombre de valeurs") plot(xlabel="Jour de l'année", ylabel="Nombre de valeurs")
histogram!(data.day, bins=0:31:365, label=nothing) histogram!(data.day, bins=0:31:365, label=nothing)
disp()
``` ```
## Analyse des variations annuelles ## Analyse des variations annuelles
...@@ -407,13 +465,13 @@ plutôt des valeurs (notées $C_0$ et $C_1$ dans le code) moyennées sur les 7 ...@@ -407,13 +465,13 @@ plutôt des valeurs (notées $C_0$ et $C_1$ dans le code) moyennées sur les 7
premiers et 7 derniers jours de l'année. premiers et 7 derniers jours de l'année.
```julia ```julia
data.phi = copy(data.val) data.phi = copy(data.co2)
data.alpha = zero(data.val) data.alpha = zero(data.co2)
data.beta = zero(data.val) data.beta = zero(data.co2)
for year in unique(data.year) for year in unique(data.year)
idx = (data.year .== year) idx = (data.year .== year)
C = data.val[idx .& (data.day .< 7) ] |> mean C₀ = data.co2[idx .& (data.day .< 7) ] |> mean
C = data.val[idx .& (data.day .> 358)] |> mean C₁ = data.co2[idx .& (data.day .> 358)] |> mean
α = C₀ α = C₀
β = (C₁ - C₀) / 365 β = (C₁ - C₀) / 365
...@@ -428,25 +486,29 @@ l'année 1982. On voit, sur la figure du haut, les mesures brutes ainsi que la ...@@ -428,25 +486,29 @@ l'année 1982. On voit, sur la figure du haut, les mesures brutes ainsi que la
tendance locale (affine). Sur la figure du bas, la composante périodique locale tendance locale (affine). Sur la figure du bas, la composante périodique locale
vérifie bien les contraintes demandées, avec ses valeurs nulles aux bords. vérifie bien les contraintes demandées, avec ses valeurs nulles aux bords.
```julia; echo=false ```julia; echo=false; results="raw"
let df = data[data.year.==1982, :] let df = data[data.year.==1982, :]
p1 = plot(ylabel="CO2 [ppm]") p1 = plot(ylabel="CO2 [ppm]")
plot!(p1, df.date, df.val, label="mesures") plot!(p1, df.date, df.co2,
plot!(p1, df.date, df.alpha.+df.beta.*df.day, label="alpha + beta * day") label="mesures", linewidth=2)
plot!(p1, df.date, df.alpha.+df.beta.*df.day,
label="alpha + beta * day", linewidth=2)
p2 = plot(xlabel="date", ylabel="CO2 [ppm]") p2 = plot(xlabel="date", ylabel="CO2 [ppm]")
plot!(p2, df.date, df.phi, label="phi") plot!(p2, df.date, df.phi,
label="phi", linecolor="green", linewidth=2)
l = @layout [a ; b] l = @layout [a ; b]
plot(p1, p2, layout=l) plot(p1, p2, layout=l)
disp()
end end
``` ```
Pour que notre décomposition soit valide, il faut que les composantes Pour que notre décomposition soit valide, il faut que les composantes
périodiques locales de chaque année se ressemblent, au point de pouvoir être périodiques locales de chaque année se ressemblent, au point de pouvoir être
représentées par leur moyenne. Nous calculons donc cette moyenne pour toutes les représentées par leur moyenne. Nous calculons donc cette moyenne (et
années du jeu de données d'étude. l'écart-type associé) pour toutes les années du jeu de données d'étude.
```julia; results="raw" ```julia; results="raw"
avg = by(data, :day, :phi=>mean, :phi=>std) avg = by(data, :day, :phi=>mean, :phi=>std)
...@@ -455,13 +517,14 @@ info(avg) ...@@ -455,13 +517,14 @@ info(avg)
Et nous traçons l'ensemble des composantes oscillantes locales aux côtés de cette moyenne. Et nous traçons l'ensemble des composantes oscillantes locales aux côtés de cette moyenne.
```julia; echo=false ```julia; echo=false; results="raw"
plot(xlabel="date", ylabel="CO2 [ppm]") plot(xlabel="date", ylabel="CO2 [ppm]")
plot!(data.day, data.phi, seriestype=:scatter, label="phi") plot!(data.day, data.phi, seriestype=:scatter, label="phi")
plot!(avg.day, avg.phi_mean, linewidth=3, label="phi (moyenne)") plot!(avg.day, avg.phi_mean, linewidth=3, label="phi (moyenne)")
plot!(avg.day, avg.phi_mean.-avg.phi_std, linecolor="red", linewidth=1, label="phi (IC 67%)") plot!(avg.day, avg.phi_mean.-2avg.phi_std, linecolor="red", linewidth=2, label="phi (IC 95%)")
plot!(avg.day, avg.phi_mean.+avg.phi_std, linecolor="red", linewidth=1, label=nothing) plot!(avg.day, avg.phi_mean.+2avg.phi_std, linecolor="red", linewidth=2, label=nothing)
plot!(avg.day, avg.phi_std, linewidth=1, label="phi (écart-type)") plot!(avg.day, avg.phi_std, linewidth=2, label="phi (écart-type)")
disp()
``` ```
Comme on pouvait s'y attendre, il reste une forte variabilité d'année en Comme on pouvait s'y attendre, il reste une forte variabilité d'année en
...@@ -486,17 +549,18 @@ $$ ...@@ -486,17 +549,18 @@ $$
```julia ```julia
data = join(data, avg, on=:day) data = join(data, avg, on=:day)
data.theta = data.val .- data.phi_mean; data.theta = data.co2 .- data.phi_mean;
``` ```
Même s'il reste des oscillations locales, nous constatons tout de même que la Même s'il reste des oscillations locales, nous constatons tout de même que la
composante tendancielle est devenue suffisamment lisse pour récupérer une forme composante tendancielle est devenue suffisamment lisse pour récupérer une forme
de monotonie. de monotonie.
```julia; echo=false ```julia; echo=false; results="raw"
plot(xlabel="Date", ylabel="CO2 [ppm]", legend=:bottomright) plot(xlabel="Date", ylabel="CO2 [ppm]")
plot!(data.date, data.val, label="mesures") plot!(data.date, data.co2, label="mesures")
plot!(data.date, data.theta, label="theta") plot!(data.date, data.theta, label="theta")
disp()
``` ```
Nous allons maintenant tenter de caractériser la tendance sous-jacente. Au vu de Nous allons maintenant tenter de caractériser la tendance sous-jacente. Au vu de
...@@ -507,8 +571,11 @@ $$ ...@@ -507,8 +571,11 @@ $$
Le paquet Julia [`GLM`](https://github.com/JuliaStats/GLM.jl) permet de *fitter* Le paquet Julia [`GLM`](https://github.com/JuliaStats/GLM.jl) permet de *fitter*
des modèles linéaires (ou modèles linéaires généralisés). Nous l'utilisons ici des modèles linéaires (ou modèles linéaires généralisés). Nous l'utilisons ici
pour estimer les paramètres de notre modèle. L'estimation est réalisée sur un pour estimer les paramètres de notre modèle. Afin de ne pas considérer un nombre
sous-échantillonnage du jeu de données, comprenant un point (moyen) par an. artificiellement élevé de points de mesure, l'estimation est réalisée sur un
sous-échantillonnage du jeu de données, comprenant un point (moyen) par an, ce
qui correspond à l'échelle de temps que nous n'avons pas encore capturé dans la
composante oscillante.
```julia; wrap=false ```julia; wrap=false
sample = by(data, :year, theta = :theta=>mean, date_num = :date_num=>mean) sample = by(data, :year, theta = :theta=>mean, date_num = :date_num=>mean)
...@@ -519,30 +586,40 @@ Le modèle obtenu semble correspondre assez bien aux données, avec une ...@@ -519,30 +586,40 @@ Le modèle obtenu semble correspondre assez bien aux données, avec une
incertitude de l'ordre de quelques pourcents sur l'estimation des paramètres incertitude de l'ordre de quelques pourcents sur l'estimation des paramètres
$\alpha$ et $\beta$. $\alpha$ et $\beta$.
```julia ```julia; results="hidden"
α, β, γ = coef(model); # Paramètres du modèle et intervalle de confiance à 95%
α1, β1, γ1 = coef(model)-2*stderror(model); α, β, γ = coef(model)
α2, β2, γ2 = coef(model)+2*stderror(model); α1, β1, γ1 = coef(model)-2*stderror(model)
α2, β2, γ2 = coef(model)+2*stderror(model)
# Taux d'accroissement
τ(d) = (β + γ * date2num(d)) * 365
τ1 = τ(Date("1959"))
τ2 = τ(Date("2015"))
``` ```
En première approximation, on a en particulier une hausse tendancielle de la Ce modèle caractérise en particulier une hausse tendancielle de la concentration
concentration de CO2 atmosphérique de l'ordre de de CO2 atmosphérique de l'ordre de $\tau_1$ = `j @printf "%.2f" τ1` ppm/an en
$\beta$ = `j @printf "%.2e" β` ppm/jour ($\pm$`j @printf "%.2f" 100*(β2-β)/β`%), 1959 et $\tau_2$ = `j @printf "%.2f" τ2` ppm/an en 2015. Si la valeur de 1959
qui se traduit en une augmentation annuelle comprise entre correspond assez bien au taux de 0.75 ppm/an déterminé dans `j monroe2015()`, la
`j @printf "%.2f" 365*β1` et `j @printf "%.2f" 365*β2` ppm/an valeur que nous trouvons pour 2015 est sous-évaluée de plus de 30% par rapport
avec un niveau de confiance de 95%. aux 2.25 ppm/an calculés par Monroe.
On voit toutefois que l'incertitude sur les paramètres (en particulier $\beta$) Par ailleurs, s'il est clair que la tendance est à l'augmentation, on voit
est de nature à engendrer une perte de prédictibilité du modèle en temps long. toutefois que l'incertitude sur les paramètres n'est pas complètement
négligeable. L'incertitude sur $\beta$ est en particulier de nature à engendrer
une perte de prédictibilité du modèle en temps long.
```julia; echo=false ```julia; echo=false; results="raw"
data.theta_pred = predict(model, data) data.theta_pred = predict(model, data)
plot(xlabel="Date", ylabel="CO2 [ppm]", legend=:bottomright) plot(xlabel="Date", ylabel="CO2 [ppm]")
plot!(data_raw.date, data_raw.val, label="mesures") plot!(data_raw.date, data_raw.co2, label="mesures")
plot!(data.date, data.theta_pred, label="modèle theta", linewidth=3) plot!(data.date, data.theta_pred, label="modèle theta", linewidth=3)
plot!(data.date, α .+ β1 * data.date_num .+ γ*data.date_num.^2, linecolor="red", label="modèle theta (IC 95%)") plot!(data.date, α .+ β1 * data.date_num .+ γ*data.date_num.^2, linecolor="red", label="modèle theta (IC 95%)")
plot!(data.date, α .+ β2 * data.date_num .+ γ*data.date_num.^2, linecolor="red", label=nothing) plot!(data.date, α .+ β2 * data.date_num .+ γ*data.date_num.^2, linecolor="red", label=nothing)
disp()
``` ```
## Reconstruction du signal complet et prédiction ## Reconstruction du signal complet et prédiction
...@@ -570,14 +647,14 @@ prediction = DataFrame(date_num=date_num, ...@@ -570,14 +647,14 @@ prediction = DataFrame(date_num=date_num,
info(prediction) info(prediction)
``` ```
Le modèle `GLM` précédemment calibré est utilisé pour calculer $\theta(t)$. On y Le modèle `GLM` précédemment calibré est utilisé pour prédire $\theta(t)$. On y
ajoute la fonction de forme annuelle $\phi(t)$ par périodicité. ajoute la fonction de forme annuelle $\phi(t)$ par périodicité.
```julia; results="raw" ```julia; results="raw"
prediction.theta = predict(model, prediction) prediction.theta = predict(model, prediction)
prediction = join(prediction, avg, on=:day) prediction = join(prediction, avg, on=:day)
prediction.val = prediction.theta .+ prediction.phi_mean prediction.co2 = prediction.theta .+ prediction.phi_mean
info(prediction) info(prediction)
``` ```
...@@ -585,12 +662,14 @@ On voit que la forme annuelle semble bien reproduite sur les 5 premières année ...@@ -585,12 +662,14 @@ On voit que la forme annuelle semble bien reproduite sur les 5 premières année
pour lesquelles il est possible de comparer les prédictions avec les mesures pour lesquelles il est possible de comparer les prédictions avec les mesures
réelles. En revanche, la tendance ne colle que sur les deux premières années réelles. En revanche, la tendance ne colle que sur les deux premières années
d'extrapolation ; on observe un décalage significatif et croissant ensuite. d'extrapolation ; on observe un décalage significatif et croissant ensuite.
```julia; echo=false
```julia; echo=false; results="raw"
idx = data_raw.date .> Date(lastyear) idx = data_raw.date .> Date(lastyear)
plot(xlabel="Date", ylabel="CO2 [ppm]", legend=:bottomright) plot(xlabel="Date", ylabel="CO2 [ppm]")
plot!(data_raw.date[idx], data_raw.val[idx], label="mesures") plot!(data_raw.date[idx], data_raw.co2[idx], label="mesures")
plot!(prediction.date, prediction.theta, label="tendance") plot!(prediction.date, prediction.theta, label="tendance", linewidth=2)
plot!(prediction.date, prediction.val, label="prédiction") plot!(prediction.date, prediction.co2, label="prédiction", linewidth=2)
disp()
``` ```
## Pour aller plus loin ## Pour aller plus loin
...@@ -599,8 +678,8 @@ plot!(prediction.date, prediction.val, label="prédiction") ...@@ -599,8 +678,8 @@ plot!(prediction.date, prediction.val, label="prédiction")
@info "* Studying alternate model" @info "* Studying alternate model"
``` ```
En utilisant des modèles plus complexes, il est possible de mieux caractériser En utilisant des techniques de calage de modèle plus complexes, il est possible
la tendance en temps longs : de mieux caractériser la tendance en temps longs :
```julia; wrap=false ```julia; wrap=false
model2 = glm(@formula(theta ~ date_num + date_num^2), sample, model2 = glm(@formula(theta ~ date_num + date_num^2), sample,
...@@ -610,7 +689,7 @@ model2 = glm(@formula(theta ~ date_num + date_num^2), sample, ...@@ -610,7 +689,7 @@ model2 = glm(@formula(theta ~ date_num + date_num^2), sample,
En reprenant l'analyse précédente, ce nouveau modèle donne les prédictions En reprenant l'analyse précédente, ce nouveau modèle donne les prédictions
suivantes, qui collent quasi-parfaitement aux mesures dans la période de test : suivantes, qui collent quasi-parfaitement aux mesures dans la période de test :
```julia; echo=false ```julia; echo=false; results="raw"
date_num = date2num(Date(lastyear)):10:date2num(today()+5Years) date_num = date2num(Date(lastyear)):10:date2num(today()+5Years)
prediction2 = DataFrame(date_num=date_num, prediction2 = DataFrame(date_num=date_num,
date=num2date.(date_num), date=num2date.(date_num),
...@@ -618,13 +697,14 @@ prediction2 = DataFrame(date_num=date_num, ...@@ -618,13 +697,14 @@ prediction2 = DataFrame(date_num=date_num,
prediction2.theta = predict(model2, prediction2) prediction2.theta = predict(model2, prediction2)
prediction2 = join(prediction2, avg, on=:day) prediction2 = join(prediction2, avg, on=:day)
prediction2.val = prediction2.theta .+ prediction2.phi_mean prediction2.co2 = prediction2.theta .+ prediction2.phi_mean
idx = data_raw.date .> Date(lastyear) idx = data_raw.date .> Date(lastyear)
plot(xlabel="Date", ylabel="CO2 [ppm]", legend=:bottomright) plot(xlabel="Date", ylabel="CO2 [ppm]")
plot!(data_raw.date[idx], data_raw.val[idx], label="mesures") plot!(data_raw.date[idx], data_raw.co2[idx], label="mesures")
plot!(prediction2.date, prediction2.theta, label="tendance") plot!(prediction2.date, prediction2.theta, label="tendance", linewidth=2)
plot!(prediction2.date, prediction2.val, label="prédiction") plot!(prediction2.date, prediction2.co2, label="prédiction", linewidth=2)
disp()
``` ```
Je n'ai aucune idée de la signification qu'il est possible de donner à ces Je n'ai aucune idée de la signification qu'il est possible de donner à ces
...@@ -632,11 +712,6 @@ valeurs de paramètres, aussi est-il sans doute préférable de conserver, au mo ...@@ -632,11 +712,6 @@ valeurs de paramètres, aussi est-il sans doute préférable de conserver, au mo
dans un premier temps, l'interprétation donnée par le modèle quadratique simple dans un premier temps, l'interprétation donnée par le modèle quadratique simple
(quoi que donnant des prédictions plus éloignées des données). (quoi que donnant des prédictions plus éloignées des données).
```julia; echo=false; results="hidden"
# close plots
closeall()
```
<!-- Local Variables: --> <!-- Local Variables: -->
<!-- mode: markdown --> <!-- mode: markdown -->
<!-- ispell-local-dictionary: "french" --> <!-- ispell-local-dictionary: "french" -->
......
<!DOCTYPE html>
<HTML lang = "en">
<HEAD>
<meta charset="UTF-8"/>
<meta name="viewport" content="width=device-width, initial-scale=1.0, user-scalable=yes">
{{#:title}}<title>{{:title}}</title>{{/:title}}
{{{ :header_script }}}
<script src="plotly-latest.min.js"></script>
<script type="text/x-mathjax-config">
MathJax.Hub.Config({
tex2jax: {inlineMath: [['$','$'], ['\\(','\\)']]},
TeX: { equationNumbers: { autoNumber: "AMS" } }
});
</script>
<script type="text/javascript" async src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
</script>
{{{ :highlightcss }}}
<style type="text/css">
{{{ :themecss }}}
</style>
</HEAD>
<BODY>
<div class ="container">
<div class = "row">
<div class = "col-md-12 twelve columns">
<div class="title">
{{#:title}}<h1 class="title">{{:title}}</h1>{{/:title}}
{{#:author}}<h5>{{{:author}}}</h5>{{/:author}}
{{#:date}}<h5>{{{:date}}}</h5>{{/:date}}
</div>
{{{ :body }}}
<HR/>
<div class="footer"><p>
Published from <a href="{{{:source}}}">{{{:source}}}</a> using
<a href="http://github.com/mpastell/Weave.jl">Weave.jl</a>
{{:wversion}} on {{:wtime}}.
<p></div>
</div>
</div>
</div>
</BODY>
</HTML>
using Pkg #!/bin/bash
Pkg.activate(@__DIR__) #=
OPTS="--color=yes --startup-file=no"
[ -f Exercice.so ] && OPTS="$OPTS -J Exercice.so"
import Weave exec julia ${OPTS} "${BASH_SOURCE[0]}" "$@"
Weave.tangle(joinpath(@__DIR__, "exercice.jmd")) =#
using PackageCompiler using Weave
create_sysimage(Symbol.(keys(Pkg.project().dependencies)); let srcfile = joinpath(@__DIR__, "exercice.jmd")
sysimage_path=Pkg.project().name * ".so", tangle(srcfile)
precompile_execution_file=joinpath(@__DIR__, "exercice.jl")) Weave.weave() = @time weave(srcfile)
weave()
end
# Local Variables:
# mode: julia
# End:
This source diff could not be displayed because it is too large. You can view the blob instead.
# Permet de construire des durées avec une syntaxe comme :
# 7Days
const Years = Year(1)
const Days = Day(1)
# Affichage des DataFrames dans la sortie HTML # Affichage des DataFrames dans la sortie HTML
function info(df::DataFrame) function info(df::DataFrame)
try
WEAVE_ARGS
catch
display(df)
return
end
pretty(x) = x pretty(x) = x
pretty(x::Float64) = @sprintf "%.6f" x pretty(x::Float64) = @sprintf "%.6f" x
pretty(::Type{Union{T,Missing}}) where {T} = "$T?" pretty(::Type{Union{T,Missing}}) where {T} = "$T?"
function printrow(i) function printrow(i)
...@@ -55,7 +68,25 @@ function info(df::DataFrame) ...@@ -55,7 +68,25 @@ function info(df::DataFrame)
end end
# Permet de construire des durées avec une syntaxe comme :
# 7Days # Meilleure intégration entre Plotly et Weave
const Years = Year(1) function disp()
const Days = Day(1) try
WEAVE_ARGS
catch
display(plot!())
return
end
buf = IOBuffer()
show(buf, MIME("text/html"), plot!())
seekstart(buf)
str = String(read(buf))
m = match(r"<body>(.*)</body>"s, str)
return if m === nothing
str
else
m.captures[1]
end |> print
end
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment