TP analyse de données avec numpy

5.3. TP analyse de données avec numpy#

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

temp_terre.png

Attention il faut exécuter toutes les cellules suivantes pour afficher son nom !!

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
# police des titres
plt.rc('font', family='serif', size='18')

from IPython.display import display, Markdown, clear_output
from validation.valide_markdown import test_markdown, test_code, test_algorithme
from validation.validation import info_etudiant

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 :
    NOM, PRENOM, NUMERO_ETUDIANT = info_etudiant()
# 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))

try:
    printmd("INITIALISATION OK!")
except:
    print("Erreur vous n'avez pas executée la cellule vide précédente !")
    print("Votre Notebook n'est pas initialisé correctement !")

5.3.1. Objectifs#

Apprendre à lire des données dans un fichier, les analyser et les traiter pour en déduire des prédictions.

La démarche consiste en

  1. lecture des données

  2. analyse des données

  3. détermination d’une loi simple

  4. prédiction à partir de cette loi

  5. analyse des résultats

  6. rédaction d’un compte rendu

application

  • problème du réchauffement climatique.

5.3.2. Étape 1: lecture des données#

En utilisant la commande unix/linux head mon_fichier, on affiche les 10 premières lignes du fichier ./data/anomalie_temperature.csv dans la cellule ci-dessous

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.

!head ./data/anomalie_temperature.csv

5.3.2.1. 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 doit aussi importer la bibliothèque numpy

import numpy as np
annee = 0
anomalie_temp = 0
### BEGIN SOLUTION
### END SOLUTION
print(annee, anomalie_temp)

assert(type(annee) == np.ndarray)
assert(type(anomalie_temp) == np.ndarray)

5.3.3. É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, et on ajoute un titre et des labels pour les axes pour obtenir un tracé de qualité (au sens d’une publication scientifique).

On sauvegarde ensuite le tracé pour pouvoir l’inclure dans le compte rendu avec l’instruction

 plt.savefig("macourbe.png")

ou macourbe.png est le nom du fichier image contenant la figure, que l’on inclu ensuite simplement dans le compte rendu avec la ligne suivante séparée par 2 lignes blanches:

 ![figure xx: titre de la figure](macourbe.png)

attention a bien sauter une ligne avec et après cette commande et à commenter chaque figure !

%matplotlib inline
import matplotlib.pyplot as plt
# police des titres
plt.rc('font', family='serif', size='18')
### BEGIN SOLUTION
# tracer
#plt.savefig("temperature.png")
### END SOLUTION

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

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\):

(5.11)#\[\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 valeur_moyenne qui calcule la moyenne pour n’importe quel tableau, que l’on peut ensuite réutiliser.

5.3.4.1. question 1: programmation#

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 écrira un test de la fonction dans la cellule suivante avant de passer le test de validation.

Dans le fichier de compte-rendu on décrira le principe de l’algorithme de calcul de la moyenne avec :

  • la formule mathématique

  • le principe de calcul avec une boucle

# fonction valeur_moyenne
### BEGIN SOLUTION
### END SOLUTION

# test de vérification de la fonction en l'appelant avec des arguments 
# pour calculer par exemple la moyenne de 11 nombres équi-répartits entre 1 et 3
### BEGIN SOLUTION
### END SOLUTION

# test de la cellule précédente
assert(test_code("TP_regression_lineaire.ipynb","cell-verif0","valeur_moyenne("))

# test validation de la fonction
assert (np.abs(np.linspace(0,1000000000,1000).mean()-valeur_moyenne(np.linspace(0,1000000000,1000)))> 1.e-9)
assert (np.abs(valeur_moyenne(np.linspace(1,3,11)) - 2.) < 1.e-10)

5.3.4.2. 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.

moyenne_annee = 0
moyenne_anomalie = 0
### BEGIN SOLUTION
### END SOLUTION

# test de validation
assert(test_code('TP_regression_lineaire.ipynb','cell-verif1','valeur_moyenne('))
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)

5.3.4.3. question 3: calcul des coefficients de la droite#

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 = anomalie_temp.

En déduire la valeur de \(a_1\) et \(a_0\)

X = annee
Y = anomalie_temp
X2 = 0
XY = 0
a0 = 0
a1 = 0
### BEGIN SOLUTION
### END SOLUTION

# test de validation du code précédent
print("a0=",a0,"a1=",a1)
assert(test_code('TP_regression_lineaire.ipynb','cell-verif2','valeur_moyenne('))
assert (np.abs(a1 - np.sum(X*(Y-np.ndarray.mean(Y)))/np.sum(X*(X-np.ndarray.mean(X)))) < 1.e-10)

5.3.4.4. question 4: lissage#

Écrivez une fonction coefficients_lissage 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 coefficients_lissage(x, y):
    """
    ecrire la description docstrings
    """

    a1 = xxx
    a0 = yyy
    
    return a1, a0

### BEGIN SOLUTION
### END SOLUTION    

# test de verification par rapport au calcul précédent
### BEGIN SOLUTION
### END SOLUTION

# test verification de la cellule précédente
assert(test_code("TP_regression_lineaire.ipynb","cell-verif02","coefficients_lissage("))

# test de validation de la fonction 
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("Erreurs: coefficients_lissage ne depends pas de valeur_moyenne")
finally:
    valeur_moyenne = orig_valeur_moyenne

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

a0 = 0
a1 = 0
reg = 0
### BEGIN SOLUTION
### END SOLUTION

# validation de la cellule précécente
print(a1,a0)
assert(test_code('TP_regression_lineaire.ipynb','cell-verif3','coefficients_lissage('))
assert((a1>0.0) and (a0 != 0.) and len(reg) == len(annee))

5.3.4.5. 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 en mettant dans la variable an l’année pour cette prédiction 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).

Dans la cellule « tracer du resultat » on tracera sur une figure la régression linéaire avec les données d’origine, ainsi que le point de prédiction. Pour ce dernier on pourra utiliser l’instruction

   plt.plot([an],[prediction],'*')

On sauvegardera la figure dans un fichier pour l’inclure ensuite dans le compte rendu, comme vu précédemment.

printmd("### Estimation de la hausse des température en {}".format(_an))

# calcul de la prediction
prediction = 0
an = 0
### BEGIN SOLUTION
### END SOLUTION

# test validation du code de la cellule précédente
print("prédiction de la hausse de temperature: {} degré en {}".format(prediction,_an))
assert (a1>0) and ((prediction*100 - np.round(prediction*100) == 0.))

# tracer la figure avec les donnees, le lissage et le point de prédiction
### BEGIN SOLUTION
#plt.savefig("prediction.png")
### END SOLUTION

5.3.5. 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 et comparer avec le calcul précédent

# calcul avec NumPy 
a1 = 0
a0 = 0
f_lineaire = 0
### BEGIN SOLUTION
### END SOLUTION

# validation du code de la cellule précédente
print(a1,a0)
assert((a1>0) and (a0!=0) and (type(f_lineaire ) == np.poly1d) and (f_lineaire.order == 1))

5.3.6. 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 slicing avec : . 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.

i_1970 = 0
annee_1 = 0
anomalie_temp_1 = 0
annee_2 = 0
anomalie_temp_2 = 0
### BEGIN SOLUTION
### END SOLUTION

# validation du code de la cellule précédente
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)

5.3.6.1. 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.

Enfin dans la cellule « tracer du résultat », on tracera sur une figure la régression linéaire avec les données d’origine, ainsi que le point de prédiction. On sauvegardera la figure dans un fichier pour l’inclure ensuite dans le compte rendu, comme vu précédemment.

# calcul de la prediction
f_lineaire_1 = 0
f_lineaire_2 = 0
prediction_2 = 0
an = 0
### BEGIN SOLUTION
### END SOLUTION

# test validation du code de la cellule précédente
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.)

# tracer du résultat avec les données, le lissage et le point de prédiction
### BEGIN SOLUTION
#plt.savefig("prediction2.png")
### END SOLUTION

5.3.7. Conclusion#

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

Ecrire votre analyse et votre conclusion dans le compte rendu en insistant sur

  1. Description de la méthode d’analyse

  2. Résultat de l’analyse

  3. Conclusion

Le compte rendu est à écrire dans le fichier CompteRendu.md :

  1. Génération de la version HTML du Compte Rendu (avec mise en page)

  • Exécution de la commande ci-dessous pour générer le fichier html

  1. Visualisation du Compte Rendu (version html)

  • Cliquez sur le lien suivant après exécution de la commande ci-dessous

  • Cliquez sur le bouton mise à jour du navigateur pour mettre à jour la page web

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)

# génération de la version html du CR
!genereTPhtml CompteRendu

# test sur le compte rendu
assert(test_markdown('CompteRendu.md',None,minm=200,maxe=0.2))

5.4. FIN du TP#

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