Commit 275baf11 authored by François Févotte's avatar François Févotte

Ex 03-3 (suite)

parent 1188bdbe
This source diff could not be displayed because it is too large. You can view the blob instead.
...@@ -7,7 +7,7 @@ date: avril 2020 ...@@ -7,7 +7,7 @@ date: avril 2020
# Gestion des dépendances # Gestion des dépendances
## Environnement ### Environnement
```julia; echo=false; results="hidden" ```julia; echo=false; results="hidden"
@info "Instantiating project" @info "Instantiating project"
...@@ -37,7 +37,7 @@ Pkg.instantiate() ...@@ -37,7 +37,7 @@ Pkg.instantiate()
``` ```
## Chargement des dépendances ### Chargement des dépendances
```julia; echo=false; results="hidden" ```julia; echo=false; results="hidden"
@info "Loading dependencies" @info "Loading dependencies"
...@@ -57,6 +57,9 @@ using Statistics ...@@ -57,6 +57,9 @@ using Statistics
using Plots; gr() using Plots; gr()
``` ```
```julia; echo=false; results="hidden"
include("utils.jl")
```
# Données d'entrée # Données d'entrée
...@@ -65,18 +68,20 @@ CO2](https://scrippsco2.ucsd.edu/data/atmospheric_co2/primary_mlo_co2_record.htm ...@@ -65,18 +68,20 @@ CO2](https://scrippsco2.ucsd.edu/data/atmospheric_co2/primary_mlo_co2_record.htm
fondons l'analyse sur le jeu de données contenant des observations hebdomadaires. fondons l'analyse sur le jeu de données contenant des observations hebdomadaires.
## Téléchargement ### Téléchargement
```julia; echo=false; results="hidden" ```julia; echo=false; results="hidden"
@info "Retrieving data" @info "Retrieving data"
``` ```
Le jeu de données est téléchargé une seule fois ; c'est une copie locale qui Le jeu de données est téléchargé une seule fois ; c'est une copie locale qui
sert à réaliser l'analyse. Ceci permet de garantir la version des données sert à réaliser l'analyse. Ceci permet de garantir la reproductibilité des
utilisées pour l'analyse, et stockée dans git aux côtés du présent document. données utilisées pour l'analyse, qui seront stockées dans git aux côtés du
présent document.
Il est possible de forcer le téléchargement en positionnant la variable Il est toutefois possible de forcer le téléchargement en positionnant la
`force_download=true` afin d'actualiser le jeu de données. variable `force_download=true`. Ceci permet notamment d'actualiser le jeu de
données.
```julia; results="hidden" ```julia; results="hidden"
const data_url = "https://scrippsco2.ucsd.edu/assets/data/atmospheric/stations/in_situ_co2/weekly/weekly_in_situ_co2_mlo.csv" const data_url = "https://scrippsco2.ucsd.edu/assets/data/atmospheric/stations/in_situ_co2/weekly/weekly_in_situ_co2_mlo.csv"
...@@ -107,7 +112,7 @@ end ...@@ -107,7 +112,7 @@ end
``` ```
## Lecture ### Lecture
```julia; echo=false; results="hidden" ```julia; echo=false; results="hidden"
@info "Parsing data" @info "Parsing data"
...@@ -134,15 +139,12 @@ data_raw = CSV.read(data_file; skipto=skip+1, header=[:date, :val]) ...@@ -134,15 +139,12 @@ data_raw = CSV.read(data_file; skipto=skip+1, header=[:date, :val])
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
la période de fin mars 1958 jusqu'à nos jours. la période de fin mars 1958 jusqu'à nos jours.
```julia; echo=false ```julia; results="raw"
info(df::DataFrame) = show(IOContext(stdout, :limit=>true, :displaysize=>(15, 100)), df);
```
```julia
info(data_raw) info(data_raw)
``` ```
## Vérification des données manquantes ### Vérification des données manquantes
```julia; echo=false; results="hidden" ```julia; echo=false; results="hidden"
@info "Checking for missing values" @info "Checking for missing values"
...@@ -155,7 +157,7 @@ faudra se poser des questions sur le traitement à apporter. ...@@ -155,7 +157,7 @@ faudra se poser des questions sur le traitement à apporter.
```julia ```julia
dates = data_raw.date dates = data_raw.date
for i in 2:length(dates) for i in 2:length(dates)
if dates[i]-dates[i-1] > Dates.Day(14) if dates[i]-dates[i-1] > 14Days
println("Missing data: ", println("Missing data: ",
dates[i-1], " - ", dates[i], dates[i-1], " - ", dates[i],
" (", dates[i]-dates[i-1], ")") " (", dates[i]-dates[i-1], ")")
...@@ -166,7 +168,7 @@ Il y a 12 périodes durant lesquelles les données sont manquantes, dont une en ...@@ -166,7 +168,7 @@ Il y a 12 périodes durant lesquelles les données sont manquantes, dont une en
particulier ayant duré 19 semaines en 1964. Le traitement devra en tenir compte. particulier ayant duré 19 semaines en 1964. Le traitement devra en tenir compte.
## Aperçu global des données ### Aperçu global des données
```julia; echo=false; results="hidden" ```julia; echo=false; results="hidden"
@info "Plotting raw data" @info "Plotting raw data"
...@@ -177,7 +179,7 @@ tendancielle de la concentration en CO2, à laquelle se superpose une oscillatio ...@@ -177,7 +179,7 @@ tendancielle de la concentration en CO2, à laquelle se superpose une oscillatio
à plus haute fréquence. à plus haute fréquence.
```julia; echo=false ```julia; echo=false
plot(xlabel="Date", ylabel="CO2 [ppm]", legend=:bottomright) plot(xlabel="Date", ylabel="CO2 [ppm]", legend=:bottomright)
plot!(data_raw.date, data_raw.val, label="measurements") plot!(data_raw.date, data_raw.val, label="mesures")
``` ```
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
...@@ -185,13 +187,22 @@ une période annuelle, avec un minimum local atteint chaque année autour du moi ...@@ -185,13 +187,22 @@ une période annuelle, avec un minimum local atteint chaque année autour du moi
d'octobre. d'octobre.
```julia; echo=false ```julia; echo=false
plot(xlabel="Date", ylabel="CO2 [ppm]", legend=:bottomright) plot(xlabel="Date", ylabel="CO2 [ppm]", legend=:bottomright)
plot!(data_raw.date[end-150:end], data_raw.val[end-150:end], label="measurements") plot!(data_raw.date[end-150:end], data_raw.val[end-150:end], label="mesures")
``` ```
# Analyse # Analyse
Dans cette analyse, nous allons tenter de séparer ces deux composantes : Dans cette analyse, nous allons tenter de séparer ces deux composantes :
composante tendancielle "lisse" et composante oscillante de période annuelle. composante tendancielle "lisse" et composante oscillante de période
annuelle. Plus précisément, en notant $C$ la concentration en CO2 et $t$ le
temps, nous cherchons une approximation des mesures sous la forme :
$$
C(t) \approx \theta(t) + \phi(t),
$$
où $\phi(t)$ est une fonction périodique, de période 1 an, donnant la *forme* de
la variation de la concentration en CO2 sur des échelles de temps courtes. $\theta(t)$ est
une fonction très régulière, idéalement un polynôme d'ordre bas, donnant la
*tendance* de la variation de concentration en CO2 en temps long.
La démarche que nous suivons est globalement la suivante : le jeu de données va La démarche que nous suivons est globalement la suivante : le jeu de données va
être découpé en périodes annuelles. Chacune de ces périodes annuelles sera être découpé en périodes annuelles. Chacune de ces périodes annuelles sera
...@@ -200,22 +211,28 @@ une composante périodique. ...@@ -200,22 +211,28 @@ une composante périodique.
Si notre hypothèse est correcte, les composantes périodiques de chaque année Si notre hypothèse est correcte, les composantes périodiques de chaque année
devraient être relativement comparables les unes aux autres, et pouvoir être devraient être relativement comparables les unes aux autres, et pouvoir être
approchées par leur moyenne. On peut ensuite obtenir la composante lisse approchées par leur moyenne : ce comportement annuel moyen constituera notre
tendancielle en éliminant la composante oscillatoire moyenne du signal d'origine. composante oscillante $\phi(t)$.
On pourra ensuite obtenir la composante lisse tendancielle en éliminant la
composante oscillatoire moyenne du signal d'origine :
$$
\theta(t) = C(t)-\phi(t).
$$
La périodicité des données (hebdomadaire) n'étant que peu adaptée à un découpage La périodicité des données (hebdomadaire) n'étant que peu adaptée à un découpage
annuel, nous allons commencer par interpoler les données à une maille annuel, nous allons commencer par interpoler les données à une maille
journalière. Ceci nous permettra de découper le jeu de données en années. journalière. Ceci nous permettra de découper le jeu de données en années.
## Travaux sur les dates ### Travaux sur les dates
Il est plus simple d'interpoler entre deux nombres qu'entre deux dates. Dans la Il est plus simple d'interpoler entre deux nombres qu'entre deux dates. Dans la
suite, nous adopterons une convention selon laquelle chaque date peut être suite, nous adopterons une convention selon laquelle chaque date peut être
représentée par le nombre de jours qui la sépare de la première mesure : représentée par le nombre de jours qui la sépare de la première mesure :
```julia; results="hidden" ```julia; results="hidden"
date2num(d::Date) = Dates.value(d - data_raw.date[1]) date2num(d::Date) = Dates.value(d - data_raw.date[1])
num2date(n::Int) = data_raw.date[1] + Day(n) num2date(n::Int) = data_raw.date[1] + n * Days
``` ```
Par exemple, pour les premières mesures : Par exemple, pour les premières mesures :
```julia; echo=false ```julia; echo=false
...@@ -246,7 +263,7 @@ présence d'années bissextiles, mais rien n'a été fait pour les traiter à pa ...@@ -246,7 +263,7 @@ présence d'années bissextiles, mais rien n'a été fait pour les traiter à pa
l'impact, de l'ordre de 1/365 une année sur 4, a été jugé négligeable a priori. l'impact, de l'ordre de 1/365 une année sur 4, a été jugé négligeable a priori.
## Interpolation à la maille journalière ### Interpolation à la maille journalière
```julia; echo=false; results="hidden" ```julia; echo=false; results="hidden"
@info "Interpolating daily data" @info "Interpolating daily data"
...@@ -272,17 +289,17 @@ dates. Nous avons maintenant 5 colonnes dans notre jeu de données interpolé : ...@@ -272,17 +289,17 @@ dates. Nous avons maintenant 5 colonnes dans notre jeu de données interpolé :
- `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 ```julia; results="raw"
data_interp = DataFrame(date_num=Int[], val=Float64[], date=Date[], year=Int[], day=Int[]) data_interp = DataFrame(date_num=Int[], val=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] > Dates.Day(14) 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
range = dates[i-1] => dates[i-1] range = dates[i-1] => dates[i-1]
else else
# Interpolation entre le début de la période et # Interpolation entre le début de la période et
# la fin (exclue car traitée en tant que début de # la fin (exclue car traitée en tant que début de
# la prochaine période) # la prochaine période)
range = dates[i-1] => dates[i]-Day(1) range = dates[i-1] => dates[i]- 1Days
end end
# Pour chaque jour dans la période considérée, # Pour chaque jour dans la période considérée,
...@@ -297,8 +314,7 @@ for i in 2:length(dates) ...@@ -297,8 +314,7 @@ for i in 2:length(dates)
end end
end end
# 3 dernières lignes, pour vérification info(data_interp)
last(data_interp, 3)
``` ```
En zoomant sur les données interpolées autour de l'une des périodes de données En zoomant sur les données interpolées autour de l'une des périodes de données
...@@ -311,69 +327,80 @@ plot(xlabel="Date", ylabel="CO2 [ppm]", legend=:topleft) ...@@ -311,69 +327,80 @@ plot(xlabel="Date", ylabel="CO2 [ppm]", legend=:topleft)
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.val[idx1], seriestype=:scatter,
label="interpolated values") 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.val[idx2], seriestype=:scatter,
label="measurements") label="mesures")
``` ```
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,
afin de tester la validité de nos résultats, nous n'allons réaliser l'analyse afin de tester la validité de nos résultats, nous n'allons réaliser l'analyse
que sur une fraction des données ; nous garderons les 5 dernières années comme que sur une fraction des données ; nous garderons les 5 dernières années comme
données de test. données de test. Les données sur lesquelles portera l'analyse couvrent donc la
période 1959-2014.
```julia ```julia; results="raw"
firstyear = minimum(data_interp.year)+1 firstyear = minimum(data_interp.year)+1
lastyear = maximum(data_interp.year)-6 lastyear = maximum(data_interp.year)-6
idx = (data_interp.year.>=firstyear) .& (data_interp.year.<=lastyear) idx = (data_interp.year.>=firstyear) .& (data_interp.year.<=lastyear)
data = data_interp[idx, :]; data = data_interp[idx, :];
```
Les données sur lesquelles portera l'analyse couvrent donc la période 1959-2014.
```julia
info(data) info(data)
``` ```
Sur ces années complètes, la composante `day` de la date devrait être Sur ces années complètes, la composante `day` de la date devrait être
globalement équirépartie entre 0 et 365, ce qui est globalement le cas. Les équirépartie entre 0 et 365, ce qui est globalement le cas. Les données
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
histogram(data.day, bins=0:31:365, label=nothing) plot(xlabel="Jour de l'année", ylabel="Nombre de valeurs")
histogram!(data.day, bins=0:31:365, label=nothing)
``` ```
## Analyse des variations annuelles ### Analyse des variations annuelles
```julia; echo=false; results="hidden" ```julia; echo=false; results="hidden"
@info "Fitting yearly oscillations" @info "Fitting yearly oscillations"
``` ```
Pour chaque année, on commence par tenter d'extraire la composante oscillante de Pour chaque année, on commence par tenter d'extraire la composante oscillante de
la mesure. En supposant que cette composante est périodique, il faut que la la mesure. Si l'on note $C_a$ la concentration en CO2 durant l'année $a$ et
valeur de cette composante en début d'année soit égale à celle de fin $d$ un jour de cette année, on cherche à écrire :
d'année. On va aussi supposer que la composante tendancielle varie linéairement
en cours d'année. Si on note $C$ la concentration en CO2 et $d$ le jour, on
cherche à écrire :
$$ $$
\forall d\in{0\ldots365}, \quad \forall d\in{0\ldots365}, \quad
C(d) C_a(d)
= \underbrace{\alpha + \beta d}_{\text{tendance locale}} = \underbrace{\theta_a(d)}_{\text{tendance locale}}
+ \underbrace{\phi(d)}_{\text{composante oscillante}} + \underbrace{\phi_a(d)}_{\text{forme locale}}
$$
ce qui correspond à une version locale de l'expression globale cherchée pour la
variation de concentration en CO2 : $\theta_a$ et $\phi_a$ donnent
respectivement la tendance et la forme de la concentration pour l'année $a$.
Cherchant une tendance d'ordre aussi bas que possible, nous supposons que
$\theta_a$ peut être approchée par un modèle affine à cette échelle de temps
courts :
$$
\theta_a(d) = \alpha_a + \beta_a \, d.
$$
Pour que la forme $\phi_a$ soit périodique, il faut que cette fonction prenne la
même valeur en début d'année qu'en fin d'année : $\phi(0)=\phi(365)$. Afin de
définir la constante $\alpha$ de manière unique, on fixe de plus :
$$
\phi(0)=\phi(365)=0.
$$ $$
avec $\phi(0)=\phi(365)$ pour garantir la périodicité de la composante
oscillante. Afin de définir la constante $\alpha$ de manière unique, on fixe de
plus $\phi(0)=\phi(365)=0$.
On obtient donc On obtient donc
$$ $$
\begin{align*} \begin{align*}
\alpha &= C(0), \\[1em] \alpha_a &= C_a(0), \\[1em]
\beta &= \frac{C(365)-C(0)}{365}, \\[1em] \beta_a &= \frac{C_a(365)-C_a(0)}{365}, \\[1em]
\phi(d) &= C(d) - \alpha - \beta d. \phi_a(d) &= C_a(d) - \alpha_a - \beta_a d.
\end{align*} \end{align*}
$$ $$
En pratique, plutôt que des valeurs ponctuelles $C(0)$ et $C(365)$, on prend En pratique, plutôt que des valeurs ponctuelles $C_a(0)$ et $C_a(365)$, on prend
plutôt des valeurs (notées $C_0$ et $C_1$ dans le code) moyennées sur les 7 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.
...@@ -402,7 +429,7 @@ vérifie bien les contraintes demandées, avec ses valeurs nulles aux bords. ...@@ -402,7 +429,7 @@ vérifie bien les contraintes demandées, avec ses valeurs nulles aux bords.
```julia; echo=false ```julia; echo=false
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="measurements") plot!(p1, df.date, df.val, label="mesures")
plot!(p1, df.date, df.alpha.+df.beta.*df.day, label="alpha + beta * day") plot!(p1, df.date, df.alpha.+df.beta.*df.day, label="alpha + beta * day")
...@@ -416,10 +443,12 @@ end ...@@ -416,10 +443,12 @@ 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 années du jeu de données d'étude. représentées par leur moyenne. Nous calculons donc cette moyenne pour toutes les
années du jeu de données d'étude.
```julia ```julia; results="raw"
avg = by(data, :day, :phi=>mean, :phi=>std); avg = by(data, :day, :phi=>mean, :phi=>std)
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.
...@@ -428,7 +457,9 @@ Et nous traçons l'ensemble des composantes oscillantes locales aux côtés de c ...@@ -428,7 +457,9 @@ Et nous traçons l'ensemble des composantes oscillantes locales aux côtés de c
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_std, linewidth=1, label="phi (écart-type)") plot!(avg.day, avg.phi_mean.-avg.phi_std, linecolor="red", linewidth=1, label="phi (IC 67%)")
plot!(avg.day, avg.phi_mean.+avg.phi_std, linecolor="red", linewidth=1, label=nothing)
plot!(avg.day, avg.phi_std, linewidth=1, label="phi (écart-type)")
``` ```
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
...@@ -438,7 +469,7 @@ moyen est suffisamment représentatif pour permettre d'extraire la composante ...@@ -438,7 +469,7 @@ moyen est suffisamment représentatif pour permettre d'extraire la composante
tendancielle lisse des données. tendancielle lisse des données.
## Analyse des variations tendancielles ### Analyse des variations tendancielles
```julia; echo=false; output="hidden" ```julia; echo=false; output="hidden"
@info "Fitting underlying trend" @info "Fitting underlying trend"
...@@ -446,11 +477,14 @@ tendancielle lisse des données. ...@@ -446,11 +477,14 @@ tendancielle lisse des données.
Nous sommes maintenant prêts à extraire la composante tendancielle des Nous sommes maintenant prêts à extraire la composante tendancielle des
mesures. Il suffit pour cela de retrancher aux données brutes la composante mesures. Il suffit pour cela de retrancher aux données brutes la composante
oscillante moyenne. oscillante moyenne : d'après notre modèle, on a en effet
$$
\theta(t) = C(t) - \phi(t).
$$
```julia ```julia
data = join(data, avg, on=:day) data = join(data, avg, on=:day)
data.smooth = data.val .- data.phi_mean; data.theta = data.val .- 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
...@@ -459,23 +493,29 @@ de monotonie. ...@@ -459,23 +493,29 @@ de monotonie.
```julia; echo=false ```julia; echo=false
plot(xlabel="Date", ylabel="CO2 [ppm]", legend=:bottomright) plot(xlabel="Date", ylabel="CO2 [ppm]", legend=:bottomright)
plot!(data.date, data.val, label="measurements") plot!(data.date, data.val, label="mesures")
plot!(data.date, data.smooth, label="composante tendancielle") plot!(data.date, data.theta, label="theta")
``` ```
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
la courbe (concave), nous proposons un modèle quadratique de la forme : la courbe (convexe), nous proposons un modèle quadratique de la forme :
$$ $$
C(t) = \alpha + \beta t + \gamma t^2 \theta(t) = \alpha + \beta t + \gamma t^2
$$ $$
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
pour estimer les paramètres de notre modèle. L'estimation est réalisée sur un
sous-échantillonnage du jeu de données, comprenant un point (moyen) par an.
```julia; wrap=false ```julia; wrap=false
model = lm(@formula(smooth ~ date_num + date_num^2), sample = by(data, :year, theta = :theta=>mean, date_num = :date_num=>mean)
by(data, :year, smooth = :smooth=>mean, date_num = :date_num=>mean)) model = GLM.lm(@formula(theta ~ date_num + date_num^2), sample)
``` ```
Ce genre de modèle semble correspondre assez bien aux données, avec un assez bon Le modèle obtenu semble correspondre assez bien aux données, avec une
niveau de confiance sur l'estimation des paramètres $\alpha$ et $\beta$. incertitude de l'ordre de quelques pourcents sur l'estimation des paramètres
$\alpha$ et $\beta$.
```julia ```julia
α, β, γ = coef(model); α, β, γ = coef(model);
...@@ -490,51 +530,108 @@ qui se traduit en une augmentation annuelle comprise entre ...@@ -490,51 +530,108 @@ qui se traduit en une augmentation annuelle comprise entre
`j @printf "%.2f" 365*β1` et `j @printf "%.2f" 365*β2` ppm/an `j @printf "%.2f" 365*β1` et `j @printf "%.2f" 365*β2` ppm/an
avec un niveau de confiance de 95%. avec un niveau de confiance de 95%.
On voit toutefois que l'incertitude sur $\beta$ en particulier est de nature à On voit toutefois que l'incertitude sur les paramètres (en particulier $\beta$)
engendrer une perte de prédictibilité du modèle en temps long. est de nature à engendrer une perte de prédictibilité du modèle en temps long.
```julia; echo=false ```julia; echo=false
data.smooth_pred = predict(model, data) data.theta_pred = predict(model, data)
plot(xlabel="Date", ylabel="CO2 [ppm]", legend=:bottomright) plot(xlabel="Date", ylabel="CO2 [ppm]", legend=:bottomright)
plot!(data_raw.date, data_raw.val, label="measurements") plot!(data_raw.date, data_raw.val, label="mesures")
plot!(data.date, data.smooth_pred, label="smooth model", 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="smooth model (CI 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)
``` ```
# Reconstruction du signal complet et prédiction ### Reconstruction du signal complet et prédiction
```julia; echo=false; output="hidden"
@info "Reconstructing and predicting complete signal"
```
Nous avons maintenant tous les éléments nécessaires afin de reconstruire Nous avons maintenant tous les éléments nécessaires afin de reconstruire
l'évolution de la concentration en CO2 selon l'évolution de la concentration en CO2 selon notre modèle :
```julia $$
date_num = date2num(Date(lastyear)):10:date2num(today()+Year(5)) C(t) = \theta(t) + \phi(t).
$$
Les 5 dernières années du jeu de données brut n'ont pas servi à calibrer notre
modèle ; nous allons les utiliser comme données de test. Nous extrapolons aussi
les valeurs des 5 prochaines années. Pour toute cette période, nous
reconstruisons les valeurs de concentration en CO2 tous les 10 jours.
```julia; results="raw"
date_num = date2num(Date(lastyear)):10:date2num(today()+5Years)
prediction = DataFrame(date_num=date_num, prediction = DataFrame(date_num=date_num,
date=num2date.(date_num), date=num2date.(date_num),
day=dayinyear.(date_num)) day=dayinyear.(date_num))
info(prediction)
```
Le modèle `GLM` précédemment calibré est utilisé pour calculer $\theta(t)$. On y
ajoute la fonction de forme annuelle $\phi(t)$ par périodicité.
prediction.smooth = predict(model, prediction) ```julia; results="raw"
prediction.theta = predict(model, prediction)
prediction = join(prediction, avg, on=:day) prediction = join(prediction, avg, on=:day)
prediction.val = prediction.smooth .+ prediction.phi_mean prediction.val = prediction.theta .+ prediction.phi_mean
info(prediction)
```
On voit que la forme annuelle semble bien reproduite sur les 5 premières années,
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
d'extrapolation ; on observe un décalage significatif et croissant ensuite.
```julia; echo=false
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]", legend=:bottomright)
plot!(data_raw.date[idx], data_raw.val[idx], label="measurements") plot!(data_raw.date[idx], data_raw.val[idx], label="mesures")
plot!(prediction.date, prediction.smooth, label="smooth model") plot!(prediction.date, prediction.theta, label="tendance")
plot!(prediction.date, prediction.val, label="predicted value") plot!(prediction.date, prediction.val, label="prédiction")
``` ```
## Pour aller plus loin ### Pour aller plus loin
```julia; echo=false
@info "Studying alternate model"
```
En utilisant des modèles plus complexes, il est possible de mieux caractériser
la tendance en temps longs :
```julia; wrap=false ```julia; wrap=false
model = glm(@formula(smooth ~ date_num), model2 = glm(@formula(theta ~ date_num + date_num^2), sample,
by(data, :year, smooth = :smooth=>mean, date_num = :date_num=>mean), InverseGaussian(), InverseSquareLink())
InverseGaussian(), InverseSquareLink())
``` ```
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 :
```julia ```julia; echo=false
date_num = date2num(Date(lastyear)):10:date2num(today()+5Years)
prediction2 = DataFrame(date_num=date_num,
date=num2date.(date_num),
day=dayinyear.(date_num))
prediction2.theta = predict(model2, prediction2)
prediction2 = join(prediction2, avg, on=:day)
prediction2.val = prediction2.theta .+ prediction2.phi_mean
idx = data_raw.date .> Date(lastyear)
plot(xlabel="Date", ylabel="CO2 [ppm]", legend=:bottomright)
plot!(data_raw.date[idx], data_raw.val[idx], label="mesures")
plot!(prediction2.date, prediction2.theta, label="tendance")
plot!(prediction2.date, prediction2.val, label="prédiction")
```
Je n'ai aucune idée de la signification qu'il est possible de donner à ces
valeurs de paramètres, aussi est-il sans doute préférable de conserver, au moins
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).
```julia; echo=false; results="hidden"
# close plots
closeall() closeall()
``` ```
......
# Affichage des DataFrames dans la sortie HTML
function info(df::DataFrame)
pretty(x) = x
pretty(x::Float64) = @sprintf "%.6f" x
pretty(::Type{Union{T,Missing}}) where {T} = "$T?"
function printrow(i)
print("<tr>")
print("<th>$i</th>")
for j in 1:ncol(df)
print("<td>$(pretty(df[i,j]))</td>")
end
print("</tr>")
end
function printrows()
i = 1
while i<=nrow(df) && i<=3
printrow(i)
i += 1
end
i>nrow(df) && return
if i < nrow(df)-2
print("<tr><td>...</td></tr>")
end
i = max(i, nrow(df)-2)
while i<=nrow(df)
printrow(i)
i += 1
end
end
println("$(nrow(df)) rows × $(ncol(df)) columns")
print("<table>")
print("<thead>")
print("<tr><th></th>")
for col in names(df)
print("<th>$col</th>")
end
print("</tr>")
print("<tr><th></th>")
for col in names(df)
print("<th>$(pretty(eltype(df[!,col])))</th>")
end
print("</tr>")
print("</thead>")
print("<tbody>")
printrows()
print("</tbody>")
print("</table>\n")
end
# Permet de construire des durées avec une syntaxe comme :
# 7Days
const Years = Year(1)
const Days = Day(1)
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