TP analyse de données

Marc BUFFAT, dpt mécanique, Université Lyon 1

[1] inspiré par le cours "Engineering Computations" du Pr L. Barba (Washington Univ.)

temp_terre.png

Attention il faut exécuter la cellule vide suivante !!

In [42]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

from IPython.display import display, Markdown, clear_output
def printmd(string):
    display(Markdown(string))
# test si numero étudiant spécifier
try: NUMERO_ETUDIANT
except NameError: NUMERO_ETUDIANT = None 
if type(NUMERO_ETUDIANT) is not int :
    printmd("## ERREUR: numéro d'étudiant non spécifié!!!")
    NUMERO_ETUDIANT = 123456
    NOM    = "toto"
    PRENOM = "toto"
# parametres spécifiques
_uid_  = NUMERO_ETUDIANT
np.random.seed(_uid_)
printmd("## Etudiant {} {}  id={}".format(NOM,PRENOM,NUMERO_ETUDIANT))
# parametres
_i1 = 5   + np.random.randint(10)
_i2 = 146 - np.random.randint(10)
_an = 2030 + np.random.randint(10)
# generation du fichier de donnees
cde="'1,{}d;{},$d'".format(_i1,_i2)
#print("cde:",cde)
!sed $cde ./data/land_global_temperature_anomaly.csv > ./data/anomalie_temperature.csv
printmd("### Estimation de la hausse des température en {}".format(_an))

Etudiant toto toto id=123456

Estimation de la hausse des température en 2031

Evolution de la Température de la Terre

Dans cette leçon, nous appliquerons tout ce que nous avons appris (et plus) pour analyser des données d'évolution de la température de la Terre au cours du temps.

Les questions adressées, qui sont brûlante dans le monde d'aujourd'hui, sont:

  1. La température de la terre augmente-t-elle?
  2. Et de combien?

Les données sur l'évolution de la température terrestre sont disponibles à partir de plusieurs sources: la NASA, le National Climatic Data Center (NCDC) et l'Université d'East Anglia au Royaume-Uni. Consultez le site de University Corporation for Atmospheric Research (UCAR) pour une discussion approfondie.

Le NASA Goddard Space Flight Center est l'une des sources de données climatiques mondiales. Ils ont produit la vidéo ci-dessous montrant une carte en couleur de l'évolution de la surface globale anomalies de température de 1880 à 2015.

Le terme anomalie de température globale désigne la différence de température par rapport à une valeur de référence ou à une moyenne à long terme. C'est une manière très utile de regarder le problème et à bien des égards meilleur que la température absolue. Par exemple, un mois d'hiver peut être plus froid que la moyenne à Washington DC, ainsi qu'à Miami, mais les températures absolues seront différentes aux deux endroits.

Méthode d'analyse

Comment pourrions-nous comprendre les tendances à partir des données sur la température mondiale?

La première étape de l'analyse de données consiste à générer des graphiques simples en utilisant Matplotlib. Nous allons regarder l'historique des anomalies de température, contenu dans un fichier, et faire notre premier tracé pour explorer ces données.

Nous allons lisser les données, puis nous les ajusterons avec une droite pour trouver une tendance.

Commençons!

Compte rendu de la préparation

Pour calculer une prédiction d'élévation de la température en 2031, nous allons mettre en oeuvre les étapes suivantes:

  1. lire les données qui sont dans le fichier ./data/anomalie_temperature.csv, et qui sont écrites dans le format ouvert et standard CSV (Comma Separated Value). On visualise le résultat sur un graphe mise en forme avec un titre, des labels et une légende.

  2. Les données étant extrêmement fluctuantes, dans un premier, on va lisser les données avec une droite pour obtenir une première prédiction:

    • les coefficients de la droite sont obtenus en minimisant l'écart quadratique moyen entre les données brutes et la prédiction avec la méthode des moindres carrées,
    • pour calculer ces coefficients, on calcule les moyennes de plusieurs tableaux,
    • pour éviter de refaire les mêmes calculs, on écrit une fonction python générique permettant de calculer la moyenne pour un tableau quelconque, en validant la fonction,
    • a partir de cette fonction, on va écrire une autre fonction pour calculer les coefficients directeurs de la droite de lissage à partir des formules de la droite des moindres carrées,
    • à l'aide de cette fonction, on calcule la droite de lissage, que l'on trace avec les données, et qui permet d'obtenir une premier tendance

  1. En observant les données, on peut constater deux tendances sur l'évolution de l'anomalie de température terrestre: à partir des années 1970, on observe une croissance plus grande du réchauffement qu'avant 1970. On peut donc améliorer la prédiction en appliquant la méthode précédente d'abord sur les données avant 1970 et ensuite sur les données après 1970.

    </font>

Étape 1: lire les données

Les données sont issues de la page Web NOAA (National Oceanic and Atmospheric Administration). Nous avons utilisé des données entre 1880 et aujourd'hui, mais vous pouvez utiliser des données actualisées à partir de ce même site.

Dans le répertoire data, nous avons un fichier nommé anomalie_temperature.csv. Ce fichier contient l'année sur la première colonne, et les moyennes des anomalies de température du sol répertoriées séquentiellement sur la deuxième colonne, à partir de 1880 jusqu'à nos jours. Nous allons charger le fichier, puis faire un premier tracé pour voir à quoi il ressemble.

En utilisant la commande unix/linux head mon_fichier, on affiche les 10 premières lignes du fichier mon_fichier

Attention vous n'avez pas tous les mêmes données, donc ne vous contentez pas de comparer les résultats bruts, mais essayer de comprendre la méthode.

In [43]:
!head ./data/anomalie_temperature.csv

1881,-0.09
1882,0.08
1883,-0.28
1884,-0.24
1885,-0.45
1886,-0.18
1887,-0.52
1888,-0.33
1889,-0.03
1890,-0.34

lecture du fichier

Pour charger les données à partir du fichier, nous utiliserons la fonction numpy.loadtxt(), qui nous permet de charger directement les données lues dans des tableaux NumPy. (

Les données sont séparées par des , et on veut enregistrer directement les données dans les deux tableaux annee et anomalie_temp. On va donc fournir des arguments à la fonction np.loadtxt pour l'utiliser (voir la documentation de loadtxt)

    annee, anomalie_temp = np.loadtxt(fichier, delimiter=',', unpack=True)

Dans la cellule suivante, écrire le code python pour lire les données dans les deux tableaux annee et anomalie_temp

On aussi importer la bibliothèque numpy

In [44]:
import numpy as np
fichier = './data/anomalie_temperature.csv'
annee, anomalie_temp = np.loadtxt(fichier, delimiter=',', skiprows=5, unpack=True)
In [45]:
assert(type(annee) == np.ndarray)
assert(type(anomalie_temp) == np.ndarray)
printmd("### validation OK")

validation OK

Étape 2: tracer des données

Commençons par charger le module pyplot de Matplotlib appelé, pour créer des tracés 2D ainsi que la commande spéciale "magique", % matplotlib inline pour avoir les figures incluses dans le notebook.

On utilise ensuite la fonction plot() pour tracer l'anomalie de température en fonction de l'année

In [46]:
%matplotlib inline
import matplotlib.pyplot as plt
plt.plot(annee, anomalie_temp);

amélioration du tracé

Nous avons maintenant un tracé basique mais sans aucune information,car vous ne pouvez pas déterminer de quel type de données il s'agit! Nous avons besoin d'étiquettes sur les axes, d'un titre et pourquoi pas d'une meilleure couleur, police et taille des graduations.

Un tracé de qualité (au sens d'une publication scientifique) doit toujours être la norme. La façon dont vous présentez vos données permettra aux autres (et probablement à vous à l'avenir) de mieux comprendre votre travail.

Nous pouvons personnaliser le style de nos tracés en utilisant des paramètres pour les lignes, le texte, la police et d'autres options de tracé. Nous définissons des options de style qui s'appliquent à tous les tracés de la session courante avec pyplot.rc() Ici, nous allons créer la police d'un type et d'une taille spécifiques (empattement et 16 points). Vous pouvez également personnaliser d'autres paramètres comme la largeur de ligne, la couleur, etc. (consultez la documentation).

In [47]:
# police des titres
plt.rc('font', family='serif', size='18')
# augmente la taille de la figure
plt.figure(figsize=(10,5))
# tracer
plt.plot(annee, anomalie_temp, color='#2929a3', linestyle='-', linewidth=1) 
plt.title('Anomalie de température\n')
plt.xlabel('Année')
plt.ylabel('anomalie en [°C]')
plt.grid();

Etape 3: régression linéaire par moindres carrés

Afin d'avoir une idée du comportement général de ces données, nous allons chercher une courbe lisse qui correspond (approximativement) aux points de mesure. Nous recherchons une courbe simple (par exemple, un polynôme) qui permet de lisser le bruit toujours présent dans les données expérimentales.

Soit $f (x)$ la fonction que nous ajusterons aux $n+1$ points expérimentaux: $ (x_i, y_i) $, $ i = 0, 1, ..., n $:

La fonction dépend de $m+1$ paramêtres $a_0, ...,a_m$ avec $m<n$ $$ f (x) = f (x; a_0, a_1, ..., a_m) $$

Nous devons choisir la forme de $ f (x) $ à priori, en inspectant les données expérimentales et en utilisant notre connaissance du phénomène étudié. Ainsi, l'ajustement de courbe se compose de deux étapes:

  1. Choisir la forme de $ f (x) $.
  2. Calculer les paramètres $a_k$ qui donneront la "meilleure approximation" par rapport aux données.

Dans la cellule suivante, on donne les éléments de théories qui permettent de comprendre le principe.

Un peu de thèorie

Lorsque les coordonnées $y$ des points de données sont bruitées, il est courant d'utiliser un lissage par moindres carrés, qui consiste à minimiser l'erreur quadratique $Err$ entre la fonction $f$ et les points de mesure:

\begin{equation} Err(a_0, a_1, ..., a_m) = \sum_{i = 0}^ {n} [y_i - f (x_i)] ^ 2 \end{equation}

par rapport à chacun des paramètres $ a_j $. La valeur des paramètres est obtenue par résolution des équations suivantes:

\begin{equation} \frac{\partial{Err}} {\partial{a_k}} = 0, \quad k = 0, 1, ..., m. \label{eq1} \end{equation}

qui traduisent la condition de minimisation $\vec{\nabla} Err = \vec{0}$

Les termes $ r_i = y_i - f (x_i) $ sont appelés résidus: ils donnent l'écart en $ x_i $ entre les données et la fonction de lissage.

Nous voulons minimiser une expression qui est la somme des carrés des résidus. Les équations $\eqref{eq1}$ sont généralement non linéaires en $ a_j $ et peuvent être difficiles à résoudre. Par conséquent, la fonction de lissage est généralement choisie comme une combinaison linéaire de fonctions de base $ f_j (x) $, que l'on spécifie :

\begin{equation} f (x) = a_0 f_0 (x) + a_1 f_1 (x) + ... + a_m f_m (x) \end{equation}

ce qui conduit à des équations $\eqref{eq1}$ linéaires. Dans le cas où la fonction de lissage est polynomiale, on choissit $ f_0 (x) = 1, \; f_1 (x) = x, \; f_2 (x) = x ^ 2 $, et ainsi de suite.

Régression linéaire

La fonction de lissage la plus simple est la régression linéaire, qui ajuste une ligne droite aux données. Dans ce cas,

\begin{equation} f (x) = a_0 + a_1x \end{equation}

La fonction à minimiser s'écrit:

\begin{equation} Err(a_0, a_1) = \sum_{i = 0}^{n} [y_i - f (x_i)]^2 =\sum_{i = 0}^ {n} (y_i - a_0 - a_1x_i)^2 \end{equation}

Les équations de minimisation de $Err$ deviennent:

\begin{equation} \frac{\partial{Err}}{\partial{a_0}} = \sum_{i = 0}^ {n} -2 (y_i - a_0 - a_1x_i) = 2 \left[a_0 (n + 1) + a_1 \sum_{i = 0}^{n} x_i -\sum_ {i = 0}^{n} y_i \right] = 0 \label{eq2} \end{equation}

et

\begin{equation} \frac{\partial{Err}}{\partial{a_1}} = \sum_{i = 0}^{n} -2 (y_i - a_0 - a_1x_i) x_i = 2 \left[a_0 \sum_{i = 0 }^{n} x_i + a_1 \sum_{i = 0}^{n} x_ {i}^ 2 - \sum_{i = 0}^{n} x_iy_i \right] = 0 \label{eq3} \end{equation}

Divisons les deux équations par $ 2 (n + 1) $ et réorganisons les termes.

Réorganisation $\eqref{eq2}$:

\begin{eqnarray} 2 \left[a_0 (n + 1) + a_1 \sum_{i = 0}^{n} x_i - \sum_{i = 0}^{n} y_i \right] & = 0 \nonumber \\ a_0 \frac{(n + 1)} {n + 1} + a_1 \frac {\sum_{i = 0}^{n} x_i} {n + 1} - \frac{\sum_{i = 0}^{n} y_i} {n + 1} & = 0 \\ \end{eqnarray}

soit $$ a_0 = \bar {y} - a_1 \bar {x} \label{eq4} $$

où $\bar{x} = \frac{\sum_{i = 0}^{n} x_i}{n + 1}$ et $\bar{y} = \frac{\sum_{i = 0}^{ n} y_i} {n + 1}$ sont respectivement la moyenne des $\left( x_i \right)$ et des $\left( y_i \right)$

Réorganisation maintenant $\eqref{eq3}$:

\begin{eqnarray} 2 \left[a_0 \sum_{i = 0}^{n} x_i + a_1 \sum_{i = 0}^{n} x_ {i}^ 2 - \sum_{i = 0}^{n} x_iy_i \right] & = 0 \\ a_0 \sum_{i = 0}^{n} x_i + a_1 \sum_{i = 0}^ {n} x_ {i}^ 2 - \sum_{i = 0}^{n} x_iy_i & = 0 \label{eq5}\\ \end{eqnarray}

Maintenant, si nous remplaçons $a_0$ de l'équation $\eqref{eq4}$ dans $\eqref{eq5}$ et réorganisons les termes:

\begin{eqnarray*} (\bar{y} - a_1 \bar{x}) \sum_{i = 0}^{n} x_i + a_1 \sum_{i = 0}^ {n} x_{i}^2 - \sum_{i = 0}^{n} x_iy_i & = 0 \\ \end{eqnarray*}

$$

En remplaçant les définitions des valeurs moyennes dans l'équation,

\begin{eqnarray*} \left[\frac{1}{n + 1} \sum_{i = 0}^{n} y_i - \frac{a_1}{n + 1}\sum_{i = 0}^{n} x_i \right] \sum_{i = 0}^{n} x_i + a_1\sum_{i = 0}^{n}x_ {i}^ 2 - \sum_{i = 0}^{n} x_iy_i & = 0 \\ \frac{1}{n + 1} \sum_{i = 0}^{n} y_i \sum_{i = 0}^{n} x_i -\frac{a_1} {n + 1} \sum_{i = 0}^{n} x_i \sum_{i = 0}^{n} x_i + a_1 \sum_{i = 0}^{n} x_{i} ^ 2 - \sum_{i = 0}^{n} x_iy_i & = 0 \\ \end{eqnarray*}

En remplaçant $\frac{\sum_{i = 0}^{n} x_i}{n + 1} = \bar{x}$,

\begin{eqnarray*} \sum_{i = 0}^{n} y_i \bar{x} - a_1 \sum_{i = 0}^{n} x_i \bar{x} + a_1 \sum_{i = 0}^{n} x_{i}^2 - \sum_{i = 0}^{n} x_iy_i = 0 \end{eqnarray*}

en mettant les termes en $ a_1 $ à gauche:

\begin{eqnarray*} a_1 \left[\sum_{i = 0}^{n} x_ {i}^ 2 -\sum_{i = 0}^{n} x_i \bar{x} \right] & = \sum_{i = 0 }^{n} x_iy_i -\sum_{i = 0}^{n} y_i \bar{x} \\ a_1 \sum_{i = 0}^{n} (x_ {i}^ 2 - x_i \bar{x}) & = \sum_{i = 0}^ {n} (x_iy_i - y_i \bar{x}) \\ a_1 \sum_{i = 0}^{n} x_ {i} (x_ {i} - \bar{x}) & = \sum_{i = 0}^ {n} y_i (x_i - \bar{x} ) \end{eqnarray*}

on obtiens la valeur de $a_1$:

\begin{eqnarray} a_1 = \frac{\sum_{i = 0}^{n} y_{i} (x_i - \bar{x})} {\sum_{i = 0}^{n} x_i (x_i - \bar{x })} \end{eqnarray}

D'où l'expression des coefficients $a_1$ et $a_0$ sont:

\begin{eqnarray} a_1 = \frac{\sum_{i = 0}^{n} y_ {i} (x_i - \bar{x})} {\sum_{i = 0}^{n} x_i (x_i - \bar{x })} \quad, \quad a_0 = \bar {y} - a_1 \bar {x} \end{eqnarray}

Allons-y!

Lissons maintenant les données d'anomalie de température par une droite, pour en déduire l'évolution de la température terrestre. Nous utiliserons la régression linéaire par moindres carrés pour calculer la pente $a_1$ et la constante $a_0$

$$ y = a_1 x + a_0 $$

avec pour expression des coefficients $a_1$ et $a_0$:

\begin{eqnarray} a_1 = \frac{\sum_{i = 0}^{n} y_ {i} (x_i - \bar{x})} {\sum_{i = 0}^{n} x_i (x_i - \bar{x })} \quad, \quad a_0 = \bar {y} - a_1 \bar {x} \end{eqnarray}

Dans notre cas, les données x sont dans le tableau annee, et les données y sont dans le tableau anomalie_temp. Pour calculer nos coefficients avec la formule ci-dessus, nous avons besoin des valeurs moyennes des données. Comme nous devons calculer la moyenne à la fois pour x ety, il est utile d'écrire une fonction Python personnalisée qui calcule la moyenne pour n'importe quel tableau, que l'on peut ensuite réutiliser.

C'est une bonne pratique de programmation d' éviter de répéter les instructions et d'écrire du code qui soit réutilisable, non seulement parce qu'il conduit à moins de temps de développement, mais aussi parce qu'il réduit les erreurs. Si vous devez faire le même calcul plusieurs fois, il vaut mieux l'encapsuler dans une fonction.

Rappelez-vous un des concepts clé de la programmation que l'on a vu ensemble: *Une fonction est une collection compacte de code qui exécute une action sur ses arguments ou données et qui renvoie le résultat.

Une fois défini, vous pouvez appeler une fonction autant de fois que vous le souhaitez. Lorsque nous appelons une fonction, nous exécutons tout le code à l'intérieur de la fonction. Le résultat de l'exécution dépend de la définition de la fonction et des valeurs qui y sont passées en arguments. Les fonctions peuvent ou non renvoyer des valeurs lors de leur dernière opération.

La syntaxe pour définir ses propres fonctions est la suivante:

def nom_fonction (arg_1, arg_2, ...):
    '''
    docstring: description de la fonction
    '''
    <corps de la fonction>
    return valeur

Le docstring d'une fonction est un message du programmeur documentant ce qu'il a construit. Les docstrings doivent être descriptifs et concis. Ils sont importants car ils expliquent (ou rappellent) l'utilisation prévue de la fonction aux utilisateurs. Vous pouvez accéder ultérieurement à la docstring d'une fonction en utilisant la fonction help () et en passant le nom de la fonction. Si vous êtes dans un notebook, vous pouvez également ajouter un point d'interrogation ? avant le nom de la fonction et exécuter la cellule pour afficher les informations d'une fonction.

Essayez par exemple

   print?
     ou
   help(print)

question 1: définition

Définir une fonction valeur_moyenne qui calcule la valeur moyenne de n'importe quel tableau passé en argument. On devra utiliser une boucle sur les valeurs du tableau et ne pas utiliser la méthode .mean() de numpy qui calcule la moyenne arithmétique d'un tableau.

On testera la fonction dans la cellule suivante avant de passer le test de validation.

In [48]:
def valeur_moyenne(tableau):
    """ Calcul la valeur moyenne d'un tableau
    
    Argument
    ---------
    tableau: Numpy array 
    
    Retour
    -------    
    moyenne des valeurs du tableau
    """
    somme = 0
    for element in tableau:
            somme = somme + element
    
    moyenne = somme / len(tableau)
    return moyenne
In [49]:
# tester votre fonction ici
X = np.array([1.,2.,3.,4.,5.])
print("Valeur moyenne X=",valeur_moyenne(X))

Valeur moyenne X= 3.0
In [50]:
_X = np.linspace(0,1000000000,1000)
assert (np.abs(_X.mean()-valeur_moyenne(_X))> 1.e-9)
printmd("Bien la fonction valeur_moyenne n'utilise pas mean")
assert (np.abs(valeur_moyenne(np.linspace(1,3,11)) - 2.) < 1.e-10)
printmd("Validation OK")

Bien la fonction valeur_moyenne n'utilise pas mean

Validation OK

question 2: utilisation

Une fois que vous avez définie la fonction valeur_moyenne(), elle devient disponible pour être utilisée avec n'importe quel tableau de n'importe quelle dimension. Essayons la avec nos données pour calculer la valeur moyenne des tableaux annee et anomalie_temp et mettre le résultat dans 2 variables moyenne_annee et moyenne_anomalie.

In [51]:
moyenne_annee = valeur_moyenne(annee)
moyenne_anomalie = valeur_moyenne(anomalie_temp)
In [52]:
print("moyenne annee=",moyenne_annee)
assert (np.abs(moyenne_annee - annee.mean()) < 1.e-10)
print("moyenne anomalie temp.=",moyenne_anomalie)
assert (np.abs(moyenne_anomalie - anomalie_temp.mean()) < 1.e-10)
printmd("### Validation OK")

moyenne annee= 1951.5
moyenne anomalie temp.= 0.05280303030303036

Validation OK

question 3: calcul des coefficients

Maintenant que nous avons des valeurs moyennes, nous pouvons calculer nos coefficients en suivant les équations précédentes. Nous calculons d'abord $ a_1 $ puis utilisons cette valeur pour calculer $ a_0 $.

Nos coefficients sont donnés par les relations: $$ a_1 = \frac{ \sum_{i=0}^{n} y_{i} (x_i - \bar{x})}{\sum_{i=0}^{n} x_i (x_i - \bar{x})} \quad , \quad a_0 = \bar{y} - a_1\bar{x} $$

Nous avons déjà calculé les valeurs moyennes des tableaux de données, mais la formule nécessite deux sommes sur les nouveaux tableaux dérivés. Récrire la formule précédente pour faire apparaître uniquement les moyenne des x, des y et des moyennes de nouveaux tableaux calculés en fonction de x et y.

Créer 2 tableaux dérivées XY et X2 pour calculer a1 en utilisant cette nouvelle expression pour les tableaux X = annee et Y = moyenne_anomalie.

En déduire la valeur de $a_1$ et $a_0$.

In [53]:
X = annee
Y = anomalie_temp
# tableaux
X2 = X**2
XY = X*Y
Xm = valeur_moyenne(X)
Ym = valeur_moyenne(Y)
X2m = valeur_moyenne(X2)
XYm = valeur_moyenne(XY)
a1 = (XYm - Xm*Ym)/(X2m - Xm*Xm)
a0 = Ym - a1*Xm
print("a0=",a0,"a1=",a1)

a0= -15.811426814783074 a1= 0.008129249216031824
In [54]:
assert (np.abs(a1 - np.sum(X*(Y-np.ndarray.mean(Y)))/np.sum(X*(X-np.ndarray.mean(X)))) < 1.e-10)
printmd("Validation OK")

Validation OK

question 4: lissage

Écrivez une fonction qui calcule les coefficients $a_0$ et $a_1$, appelez la fonction pour les calculer et comparez le résultat avec les valeurs que nous avons obtenues auparavant. À titre indicatif, nous vous donnons la structure que vous devez suivre:

def coeff_lissage(x, y):
    """
    ecrire la description docstrings
    """

    a1 = 
    a0 = 

    return a1, a0
In [55]:
def coefficients_lissage(x, y):
    """
    calcul les coefficients de la droite de moindres carrés lissant les points x,y
    renvoie les valeurs a1,a0  (droite a1*x+a0)
    """
    x2 = x**2
    xy = x*y
    xm = valeur_moyenne(x)
    ym = valeur_moyenne(y)
    x2m = valeur_moyenne(x2)
    xym = valeur_moyenne(xy)
    a1 = (xym - xm*ym)/(x2m - xm*xm)
    a0 = ym - a1*xm
    
    return a1, a0
In [56]:
# test par rapport au calcul précédent
coefficients_lissage(annee,anomalie_temp)
Out[56]:
(0.008129249216031824, -15.811426814783074)
In [57]:
"""test si coefficients_lissage depends de valeur_moyenne."""
from nose.tools import assert_raises
orig_valeur_moyenne = valeur_moyenne
del valeur_moyenne
try:
    assert_raises(NameError, coefficients_lissage, annee, anomalie_temp)
except AssertionError:
    raise AssertionError("coefficients_lissage ne depends pas de valeur_moyenne")
finally:
    valeur_moyenne = orig_valeur_moyenne
printmd("### Validation OK")

Validation OK

Nous avons maintenant les coefficients d'une fonction linéaire qui lisse les données. Nous pouvons maintenant calculer les valeurs prédites d'anomalie de température $ f (x) $ avec $$ f(x) = a_0 + a_1 x$$

Appelons reg le tableau obtenu en évaluant $ f (x) $ pour toutes les valeurs du tableau année.

Calculer les valeurs de $a_0$ et $a_1$, puis les valeurs du tableau reg

In [58]:
a1, a0 = coefficients_lissage(annee,anomalie_temp)
reg = a0 + a1 * annee
In [59]:
assert(len(reg) == len(annee))
printmd("### Validation OK")

Validation OK

Avec les valeurs de notre régression linéaire, nous pouvons la tracer avec des données d'origine pour voir à quoi elles ressemblent. Etudiez le code ci-dessous.

In [60]:
plt.figure(figsize=(10, 5))

plt.plot(annee, anomalie_temp, color='#2929a3', linestyle='-', linewidth=1, alpha=0.5) 
plt.plot(annee, reg, 'k--', linewidth=2, label='Regression linéaire')
plt.xlabel('Annee')
plt.ylabel('Anomalie de temperature [°C]')
plt.legend(loc='best', fontsize=15)
plt.grid();

question 5: prédiction

On va utiliser les données précédentes pour obtenir une prédiction de la hausse des températures (si la tendance actuelle se confirme) dans l'année donnée ci dessous.

Dans la cellule suivante, calculer la valeur de cette prédiction dans la variable prediction et faite une petite analyse du résultat. On utilisera la fonction numpy np.round(val,2) qui arrondit une valeur val à 2 chiffres (on expliquera pourquoi).

In [61]:
printmd("### Estimation de la hausse des température en {}".format(_an))

Estimation de la hausse des température en 2031

In [62]:
prediction = 0

prediction = np.round(a0 + a1*2031,2)
In [63]:
print("prédiction de la hausse de temperature: {} degré en {}".format(prediction,_an))
assert (prediction*100 - np.trunc(prediction*100) == 0.)
printmd("### Validation OK")

prédiction de la hausse de temperature: 0.7 degré en 2031

Validation OK

question 6: analyse

En observant le graphique précédent, on constate que la droite de lissage est plutot en dessous des données pour les dernières années, et donc on a une sous-estimation de l'élévation de température.

En observant les données, on peut constater deux tendances sur l'évolution de l'anomalie de température terrestre: à partir des années 1970, on observe une croissance plus grande du réchauffement qu'avant 1970. On peut donc améliorer la prédiction en appliquant la méthode précédente d'abord sur les données avant 1970 et ensuite sur les données après 1970.

</font>

Etape 4: utilisation de Numpy

Ci-dessus, nous avons codé la régression linéaire à partir de zéro. Mais devinez quoi: NumPy a des fonctions intégrées qui font déja ce dont nous avons besoin!

Oui! Python et NumPy sont là pour vous aider! Avec la fonction polyfit(), on obtient les coefficients de la regression linéair: la pente $a_1$ et la constante $a_0$. Et avec poly1d(), nous pouvons construire la fonction linéaire à partir de $a_1$ et $a_0$.

  • a1,a0 = np.polyfit(X,Y,degre) lissage des points X,Y avec un polynome de degré 1
  • f_lineaire = np.poly1D(a1,a0) définition du polynôme de lissage $a_1*x+a_0$ comme fonction python, i.e. f_lineaire(x) calcul la valeur au point x

Calculer ces quantités dans la cellule suivante:

In [64]:
# calcul avec NumPy 
a1, a0 = np.polyfit(annee, anomalie_temp, 1)
f_lineaire = np.poly1d((a1, a0)) 
print(a1,a0)

0.00812924921603105 -15.811426814781568
In [65]:
assert (type(f_lineaire ) == np.poly1d)
printmd("### Validation OK")

Validation OK

In [66]:
plt.figure(figsize=(10, 5))

plt.plot(annee, anomalie_temp, color='#2929a3', linestyle='-', linewidth=1, alpha=0.5) 
plt.plot(annee, f_lineaire(annee), 'k--', linewidth=2, label='Regression linéaire')
plt.xlabel('Annee')
plt.ylabel('Anomalie de temperature [°C]')
plt.legend(loc='best', fontsize=15)
plt.grid();

Amélioration du modèle de prédiction

Si vous regardez le graphique ci-dessus, vous remarquerez peut-être que vers 1970, la température commence à augmenter plus rapidement que la tendance précédente. Peut-être qu'une seule droite ne nous donne pas un ajustement assez bon.

D'où l'idée de diviser les données en deux (avant et après 1970) et de faire une régression linéaire sur chaque partie?

Pour ce faire, nous devons d'abord trouver la position dans notre tableau annee où se trouve l'année 1970. Heureusement, NumPy a une fonction pour cela appelée numpy.where() qui peut nous aider. Nous passons une condition et numpy.where() nous indique où dans le tableau la condition est vrai (True). On rappelle que le test d'égalité en Python est noté ==.

Ensuite pour diviser les données, nous utilisons la notation de slicingavec : . N'oubliez pas qu'un deux-points entre deux indices indique une plage de valeurs d'un début à une fin. La règle est que [début:fin] inclut l'élément au début de l'index mais exclut celui à la fin de l'index.

Nous savons maintenant comment diviser nos données en deux partie, pour obtenir les deux droites de régression. Nous avons besoin de deux tranches des tableaux annee et anomalie_temp, que nous mettons dans de nouvelles variables. Puis, nous effectuons deux régressions linéaires à l'aide des fonctions NumPy apprises ci-dessus.

Déterminer l'indice de l'année 1970 et la mettre dans la variable i_1970. Vérifier le résultat. Calculer les 2 tranches des tableaux annee et anomalie_temp dans annee_1 , anomalie_temp_1 et annee_2 , anomalie_temp_2 respectivement.

In [67]:
# determiner la position de l'année 1970
i_1970 = np.where(annee == 1970)[0][0]

annee_1 , anomalie_temp_1 = annee[0:i_1970], anomalie_temp[0:i_1970]
annee_2 , anomalie_temp_2 = annee[i_1970:],  anomalie_temp[i_1970:]
In [68]:
print("Indice annee 1970 :",i_1970)
assert (len(annee_1) == i_1970 and len(anomalie_temp_1) == i_1970)
assert ((len(annee)-len(annee_2)) == i_1970 and (len(anomalie_temp)-len(anomalie_temp_2)) == i_1970)
printmd("### Validation OK")

Indice annee 1970 : 84

Validation OK

nouveau modèle

en utilisant les données précédentes, calculer les 2 fonctions f_lineaire_1 et f_lineaire_2 correspondant au nouveau lissage des données.

En déduire une nouvelle estimation de la hausse des températures dans la variable prediction_2 avec une précision de 2 chiffres

In [69]:
m1, b1 = np.polyfit(annee_1, anomalie_temp_1, 1)
m2, b2 = np.polyfit(annee_2, anomalie_temp_2, 1)

f_lineaire_1 = np.poly1d((m1, b1))
f_lineaire_2 = np.poly1d((m2, b2))
prediction_2 = np.round(f_lineaire_2(_an),2)
In [70]:
assert (type(f_lineaire_1 ) == np.poly1d)
assert (type(f_lineaire_2 ) == np.poly1d)
print("prédiction de la hausse de temperature: {} degré en {}".format(prediction_2,_an))
assert (prediction_2*100 - np.trunc(prediction_2*100) == 0.)
printmd("### Validation OK")

prédiction de la hausse de temperature: 1.03 degré en 2031

Validation OK

In [71]:
plt.figure(figsize=(10, 5))

plt.plot(annee, anomalie_temp, color='#2929a3', linestyle='-', linewidth=1, alpha=0.5) 
plt.plot(annee_1, f_lineaire_1(annee_1), 'g--', linewidth=2, label='1880-1969')
plt.plot(annee_2, f_lineaire_2(annee_2), 'r--', linewidth=2, label='1970-2020')

plt.xlabel('Annee')
plt.ylabel('Annomalie de temperature [°C]')
plt.legend(loc='best', fontsize=15)
plt.grid();

Conclusion

Nous avons obtenu deux courbes différentes pour mieux décrire nos données.

Ecrire votre analyse et votre conclusion en utilisation la cellule de texte suivante au format markdown, en utilisant le

       ## Description de la méthode d'analyse

       ## Résultat de l'analyse

       ## Conclusion

En particulier comparez vos résultats avec plusieurs de vos camarades. Que pouvez vous en conclure.

Discutez aussi la validité de l'analyse (regarder en particulier la continuité du lissage)

Analyse et conclusion

  1. En utilisant la méthode de lissage par la droite des moindres carrées, on a pu obtenir un premier modèle d'élévation de la température de la terre.
  2. Pour cela nous avons lu les données dans des tableaux numpy, puis créer des fonctions Python pour automatiser le calcul des coefficients de notre modèle.
  3. En analysant le résultat, on constate cependant qu'avec notre premier modèle la prédiction d'augmentation de la température de 0.7 degré en 2031 est fortement sous-estimée.
  4. En observant les données, on observe une croissance plus grande du réchauffement à partir des années 1970 qu'avant. On peut donc améliorer la prédiction en appliquant la méthode précédente d'abord sur les données avant 1970 et ensuite sur les données après 1970.
  5. En effectuant deux lissages sur les données avant et apres 1970, on obtient une nouvelle prédiction du réchauffement de 1.03 degré en 2031, qui est supérieure à notre première analyse.

  6. On peut enfin noter que notre second modèle avec deux droites de lissage présente cependant un petit défault, qui est la non continuité du modèle en particulier pour l'année 1970. Une amélioration du modèle passerait par la prise en compte de cette continuité

    </font>

FIN du TP

In [72]:
# version
from platform import python_version,uname,platform
print("Systeme       :",uname())
print("OS            :",platform())
print("version python:",python_version())
print("version numpy :",np.__version__)

Systeme       : uname_result(system='Linux', node='idefix2', release='5.4.0-65-generic', version='#73-Ubuntu SMP Mon Jan 18 17:25:17 UTC 2021', machine='x86_64', processor='x86_64')
OS            : Linux-5.4.0-65-generic-x86_64-with-glibc2.29
version python: 3.8.5
version numpy : 1.19.4