L’idée de ce TP est de constater combien des modules/bibliothèques adaptés peuvent fournir des outils puissants et permettre un gain de temps gigantesque.
On va se placer dans un des champs les plus porteurs actuellement (et où python est très utilisé), l’analyse de données.

Cliquez sur cette invitation pour récupérer le repository du TP.

Exploration d’un jeu de données

Statistiques simples

import pandas as pd # bibliothèques dédiée au traitement de jeux de données
import matplotlib.pyplot as plt # bibliothèque graphique
import seaborn as sns # bibliothèque graphique reposant sur matplotlib et dédiée plus particulièrement à la représentation de jeux de données
import numpy as np # bibliothèque puissante permettant de gérer des tableaux multidimensionnels
import plotly.express as px # libraire permettant des graphes interactifs
import plotly.graph_objects as go # complémentaire à la première (seulement utile dans les cas complexes)

Pour pouvoir être importé, un module doit avoir été préalablement installé. Les plus importants sont installés par défaut dans certaines distributions (comme Anaconda).
Les gros modules sont généralement importés sous la forme import module as x où x est un raccourci pour le nom du module (np pour numpy ou plt pour matplotlib.pyplot). Se référer au cours Python pour les autres formes d’importation.
Pour obtenir de l’aide sur un module, on peut demander à Python (help(pd) par exemple pour avoir de l’aide sur pandas ou help(pd.read_csv) pour avoir de l’aide sur la fonction spécifique read_csv), mais il y a généralement beaucoup moins indigeste : l’aide en ligne des modules (pour Pandas par exemple).

# paramètres par défaut pour les graphes
plt.rcParams['figure.figsize'] = (15, 6)
plt.rcParams['font.family'] = "serif"
plt.rcParams['font.size'] = 13
sns.set_style("white")

Le premier jeu de données qu’on va utiliser est issu du World Happiness report (une publication annuelle de l’ONU mesurant le degrés de bonheur de la population mondiale par pays à partir de sondages).

url = "https://raw.githubusercontent.com/Info-TSI-Vieljeux/s1-tp3/main/2020.csv"
data_monde = pd.read_csv(url,sep=";",index_col=0) # data_monde est une dataframe Pandas
# Une dataframe est une sorte de dictionnaire dont les clés sont les en-têtes des colonnes et dont les lignes sont indexées.

data_monde

Précisions sur ces données :

• le score de bonheur est un score sur 10 correspondant à la moyenne des réponses des sondés (0 correspond à la pire vie possible et 10 à la meilleure)
• ce n’est pas le PIB par habitant mais son logarithme qui est utilisé pour ne pas avoir des valeurs sur des ordres de grandeur trop différents d’une colonne à l’autre
• entraide sociale : moyenne des réponses à la question binaire “en cas de difficultés, pouvez-vous compter sur de la famille ou des amis pour vous aider ?” (0 : non, 1 : oui)
• liberté des choix de vie : moyenne des réponses à la question binaire “êtes-vous satisfait ou non de votre liberté à choisir ce que vous voulez faire de votre vie ?” (0 : non, 1 : oui)
• générosité : moyenne des réponses à “Avez-vous donné à une association caritative le mois dernier ?” ajustée par rapport au PIB par habitant (valeur résiduelle)
• corruption perçue : moyenne des réponses à la question binaire “la corruption est-elle répandue dans le gouvernement ?” (0 : non, 1 : oui)

On simplifie un peu le jeu de données en retirant la colonne ‘Écart-type’ et ‘Score de bonheur en distopie’ (score minimal obtenu).

data_monde.drop(columns=['Écart-type','Score de bonheur en Distopie'], inplace=True)
data_monde.head(3)

data_monde.tail(3)

Traçons un histogramme brut du jeu de données complet pour y voir plus clair (la librairie Seaborn rend cela très simple).

sns.histplot(data=data_monde)

La méthode describe s’appliquant à des dataframe pandas retourne un résumé statistique très pratique des données de chaque colonne :

data_monde.describe()

Pour confirmer certaines des valeurs, vous allez construire différentes fonctions :

• une fonction decompte qui retourne le nombre d’éléments d’une liste,
• une fonction moyenne qui retourne la moyenne des éléments d’une liste,
• une fonction mediane qui retourne la médiane des éléments d’une liste triée en ordre croissant.

L’utilisation de fonctions statistiques déjà existantes est bien sûr prohibée.

def decompte(L):
    """
    decompte(L: liste) -> entier
    """
    # CODE

def moyenne(L):
    """
    decompte(L: liste) -> flottant
    """
    # CODE
    
def mediane(L):
    """
    decompte(L: liste) -> floattant ou entier (suivant les valeurs de L)
    """
    # CODE

Calculez, pour les 3 formes d’importation du module, l’écart-type des éléments de la liste Liste_scores en utilisant la fonction stdev du module statistics.
Il s’agit d’évaluer directement l’expresion (le nombre doit s’afficher sous la cellule sans utiliser de print).

import statistics
# CODE

from statistics import *
# CODE

import statistics as st
# CODE

Tracons maintenant un diagramme en batons des scores de bonheur des 60 premiers pays.

fig,ax = plt.subplots(figsize=(20,4))
sns.barplot(ax = ax,x = data_monde.index[:60], y = data_monde['Score de bonheur'].head(60))
plt.xticks(rotation=90)
ax.set_xlabel('')

On remarque que les pays sont classés par score de bonheur décroissant dans le jeu de données d’origine.
Mais on peut évidemment choisir un autre critère de classement si on le désire :

data_monde.sort_values(by="PIB par habitant (log)",ascending=True).head(10)

data_monde.sort_values(by="Corruption perçue",ascending=False).head()

data_monde.sort_values(by="Générosité",ascending=False).iloc[[45]]

D’après la cellule précédente, le 46^e (le 1^er est à l’indice 0) meilleur score de générosité appartient au Danemark.

Quel pays correspond à la 59^e plus courte espérance de vie en bonne santé ?

On peut aussi aisément filtrer le jeu de données en fonction de n’importe quel critère :

data_monde[(data_monde["Espérance de vie en bonne santé"]>60) & (data_monde["Espérance de vie en bonne santé"]<61)]
# Rq : pandas nécessite les opérateurs logiques bit à bit '&' (et) et '|' (ou) 
# plutôt que les opérateurs élément par élément 'and' et 'or' qui lèveraient une erreur.

Quel pays possède un score de bonheur inférieur à 5 malgré une valeur de corruption perçue inférieure à 0.5 ?

Pour récupérer l’ensemble des données d’un pays en particulier, on utilise :

data_monde.loc['France']

Région du monde                    Western Europe
Score de bonheur                           6.6638
PIB par habitant (log)                  10.584223
Entraide sociale                         0.937104
Espérance de vie en bonne santé         73.801933
Liberté des choix de vie                 0.825468
Générosité                              -0.130642
Corruption perçue                        0.583521
Name: France, dtype: object

Pour chaque variable mesurée (chaque colonne), on peut facilement tracer des histogrammes illustrant la répartition des valeurs.

sns.displot(data_monde, x="Score de bonheur", bins=20,  kde=True, height=4, aspect=3)
# bins contrôle le nombre de classes

On peut faciliter la lecture des graphes en les rendant interactif.

On utilise pour cela la bibliothèque Plotly express qui sait (comme seaborn) parler à une dataframe pandas.
On peut zoomer sur ces graphiques interactifs et obtenir des informations en survolant avec le curseur.

px.histogram(data_monde,'Corruption perçue',nbins=40,title="Corruption perçue")
# Cette fois-ci, le nombre de classes est désigné par nbins.

Modifez le graphe précedent pour répondre à cette question : combien la classe la plus peuplée de l’histogramme de l’espérence de vie en bonne santé compte-elle de valeurs si l’histogramme comporte 30 classes ?

Regroupement des données

On remarque que le jeu de données contient une colonne catégorielle : “Région du monde”.
Cela va nous permettre d’explorer de possibles dynamiques régionales : est-ce que les pays d’une même zone ont des indicateurs semblables ?

pd.unique(data_monde["Région du monde"]) # permet d'afficher une seule fois chacune des valeurs différentes de la colonne

array(['Western Europe', 'North America and ANZ','Middle East and North Africa', 'Latin America and Caribbean','Central and Eastern Europe', 'East Asia', 'Southeast Asia','Commonwealth of Independent States', 'Sub-Saharan Africa','South Asia'], dtype=object)

Traçons des diagrammes en boîte à moustaches représentant les scores de bonheur pour chacune des régions.

À nouveau Seaborn rend cela très simple…

sns.set_style("white")
fig, ax = plt.subplots(figsize=(12,8))
sns.boxplot(ax = ax, x="Score de bonheur", y="Région du monde", palette="husl", data=data_monde)
sns.despine(offset=10, trim=True)
ax.set_ylabel('')

Traçons maintenant un graphe plus général représentant toutes les relations possibles entre deux axes du jeu de données pour voir si certaines combinaisons discriminent plus nettement les différentes régions.

# Un peu long à s'exécuter (environ 30 s)
g = sns.pairplot(data_monde, hue="Région du monde", corner=True)
g._legend.set_bbox_to_anchor((0.6, 0.8))

On constate que les groupes régionaux sont relativement homogènes pour la plupart des critères.

Zoomons sur un de ces graphes :

sns.set_style("whitegrid")
sns.jointplot(data=data_monde,x="PIB par habitant (log)", y="Score de bonheur", hue="Région du monde", kind='scatter', height=8, legend=False)

Une version interactive du même graphique permet de consulter les informations pour chaque point :

px.scatter(data_monde,x='PIB par habitant (log)', y='Score de bonheur', hover_name=data_monde.index, color='Région du monde')

Trouvez la région du monde représentée sur le graphe suivant (le graphe interactif permet de trouver la réponse facilement).

Allons maintenant au-delà de la proximité géographique pour regrouper les pays en 3 grands blocs socioéconomiques : “Nord”, “Sud”, “Intermédiaire”.

conditions = [(data_monde['Région du monde'] == 'Western Europe') | (data_monde['Région du monde'] == 'North America and ANZ'),(data_monde['Région du monde'] == 'South Asia') | (data_monde['Région du monde'] == 'Sub-Saharan Africa')]
choix = ['"Nord"', '"Sud"']
data_monde['Groupe'] = np.select(conditions, choix, default='Autres')
deux_gpes = data_monde[data_monde["Groupe"].isin(['"Nord"','"Sud"'])]

# Un peu long à s'exécuter (environ 30 s)
sns.set_style("white")
g = sns.PairGrid(data_monde, diag_sharey=False, hue="Groupe")
g.map_upper(sns.scatterplot)
g.map_lower(sns.kdeplot,common_norm=False)
g.map_diag(sns.histplot,bins=20,kde=True)
g.add_legend(title="Grands groupes",adjust_subtitles=True)

L’homogénéité de ces 3 groupes saute aux yeux.

Corrélations

Les graphiques précédents mettent en évidence des corrélations assez fortes entre certaines grandeurs.
Creusons un peu.

g = sns.PairGrid(data_monde, y_vars=["Score de bonheur"], x_vars=["PIB par habitant (log)", "Corruption perçue"], height=7, aspect=1.5)
g.map(sns.regplot)

On constate sur cet exemple que le score de bonheur est corrélé positivement avec le PIB par habitant et négativement avec le degré de corruption perçue.

Pour avoir un panorama complet, traçons la matrice de corrélation donnant, pour chaque couple de variable, la valeur du coefficient de corrélation $r$ (valeur entre -1 et 1 traduisant le degré de dépendance linéaire entre deux variables) :

fig, ax = plt.subplots(figsize=(12,10))   
cmap = sns.diverging_palette(0, 230, 90, 60, as_cmap=True).reversed() # choix de la palette de couleurs
sns.heatmap(data_monde.iloc[:,1:].corr(), cmap=cmap, center=0, annot=True, fmt=".2f", linewidth = 0.5, ax=ax)

Citez les deux variables les moins corrélées entre elles (donner les noms exacts tels qu’ils apparaissent dans les données, attention à la casse). L’ordre des variables n’est pas important.

Fin du TP3a

Un chouïa d’apprentissage automatique
(machine learning)

On a vu qu’un regroupement des données en 3 grands groupes “Nord”, “Sud” et “Intermédiaire” semble plutôt cohérent.
Mais pourquoi pas laisser un algorithme décider lui-même de qui va le mieux ensemble ? Ensuite nous pourrons vérifier si cela recoupe notre découpage fait à la main.
On appelle cela un apprentissage non supervisé.

Nous allons utiliser l’algorithme des k-moyennes pour partitionner automatiquement nos données.
Il consiste à placer chaque point de données dans un espace à $n$ dimensions où $n$ est le nombre de variables (les descripteurs) et chercher à les regrouper en clusters en fonction de leurs distances.
Chaque variable correspondant à un axe du repère.

Pour aider l’algorithme, on peut tenter de réduire la dimension de l’espace dans lequel chaque point de données est plongé en utilisant une analyse en composantes principales.
L’idée est de déterminer les combinaisons des différentes variables expliquant le mieux la variance des données. Chaque nouvel axe ainsi formé (les composantes principales) explique une part décroissante mais complémentaire de la variance (sur la deuxième composante, les données sont moins étalées que sur la première, mais elles s’étalent dans une direction orthogonale, et ainsi de suite).
Projeter les données sur les premières composantes permet de les étaler le plus possible. On peut ainsi réduire l’espace à n dimensions du départ à un espace de seulement 2 ou 3 dimensions expliquant la majorité de la variance des données.

L’animation suivante montre comment serait sélectionné l’axe de la composante principale dans un espace à deux dimensions : il correspond à la position de la droite pour laquelle la distance cumulée de tous les points à la droite est la plus grande.

La bibliothèque Scikit-learn, destinée à l’apprentissage automatique, contient tout ce qu’il nous faut :

from sklearn.decomposition import PCA # l'algorithme d'analyse en composantes principales (PCA en anglais)
from sklearn.preprocessing import StandardScaler # pour centrer-réduire les données
from sklearn.cluster import KMeans # l'algorithme des k-moyennes

variables = data_monde.columns.values[1:-1]
scaler = StandardScaler()
X = scaler.fit_transform(data_monde[variables]) 
# chaque vecteur correspondant à chacune des variables est maintenant centré-réduit
pca = PCA()
components = pca.fit_transform(X)

Quelle combinaison des variables de départ utilise la première composante ? Les quelqus lignes suivantes permettent de le déterminer.

data = data_monde.copy() # pour pouvoir revenir sur le graphe suivant même après ajout de colonnes à data_monde

n_c = 1 # numéro de la composante principale à décrire
px.bar(components.T, x=data.columns.values[1:-1], y=n_c-1, labels={f"{n_c-1}": f"Composante Principale (CP) {n_c}"})

Quelle est le nom de la variable participant le plus à la composante principale n°34 ?

Représentons le pourcentage de variance expliquée par chacune des composantes :

exp_var_cumul = np.cumsum(pca.explained_variance_ratio_)
fig = px.bar(x=range(1, exp_var_cumul.shape[0] + 1),y=pca.explained_variance_ratio_,labels={"x": "composante", "y": "% variance expliquée"})
fig.add_scatter(x=list(range(1, exp_var_cumul.shape[0] + 1)), y=exp_var_cumul, name="", showlegend=False)

Les trois premières composantes expliquent plus de 80% de la variance !

Plaçons les données dans un espace réduit à ces 3 dimensions :

px.scatter_3d(components, x=0, y=1, z=2, 
              color=data_monde['Groupe'],
              labels={'0': 'CP 1', '1': 'CP 2', '2': 'CP 3'},
              hover_name=data_monde.index)

On constate à nouveau que nos 3 groupes discriminent plutôt très bien nos données même si quelques chevauchements existent.

C’est le moment d’utiliser l’algorithme des k-moyennes pour essayer de former 3 groupes homogènes :

# on ne garde que les 3 premières composantes principales
pca = PCA(n_components = 3)
pca.fit(X)
score_pca = pca.transform(X)
kmeans_pca = KMeans(n_clusters=3,init='k-means++',random_state=42)
kmeans_pca.fit(score_pca)
data_monde["Cluster"]=kmeans_pca.labels_.astype(str)
data_monde.head(3)

fig = px.scatter_3d(components, x=0, y=1, z=2, 
              color=data_monde['Cluster'],
              labels={'0': 'CP 1', '1': 'CP 2', '2': 'CP 3'},
              color_discrete_sequence=px.colors.qualitative.Bold,
              hover_name=data_monde.index)
fig.update_layout(legend_title = "Cluster")

Les 3 clusters créés reproduisent à peu de chose près les 3 groupes “Nord”, “Sud”, “Intermédiaire” construits à la main.

À quel cluster correspondent approximativement les pays du groupe “Sud” ?

Mais l’accord n’est pas parfait !

Citez un pays qui appartient au groupe “Nord” mais qui n’appartient pas au cluster lui correspondant.

Nous allons voir dans la prochaine partie du TP comment représenter ces données sur une carte pour y voir plus clair.

Fin du TP3b

Un peu de géographie

Le module suivant va permettre d’ajouter à nos données le code à 3 lettres (SO 3166-1 alpha-3) de chaque pays.
Mais pourquoi donc ? plotly express permet de tracer la carte d’un pays directement à partir de ce petit code de 3 lettres !

import country_converter as coco

iso3 = coco.convert(names=data_monde.index, to='ISO3', not_found=None)
data_monde["code"] = iso3
data_monde.head()

fig = px.choropleth(data_monde,
                    locations = "code",
                    color = "Score de bonheur",
                    projection = "orthographic",
                    color_continuous_scale = "Spectral_r",
                    hover_name = data_monde.index,
                    hover_data = {"code" : False})

fig.update_geos(
    showland = True, landcolor = "LightGrey",
    showocean = True, oceancolor = "LightBlue",
    showlakes = True, lakecolor = "LightBlue",
    showframe = False)

fig.update_layout(margin={"r":0,"t":0,"l":0,"b":0})
fig.show()

On est maintenant paré pour représenter les 3 clusters obtenus par l’algo des k-moyennes du tp3b.

data_monde["Cluster"] = [f'n°{cluster}' for cluster in data_monde["Cluster"].astype('int64') if cluster != 'nan']
data_monde.head(1)

fig = px.choropleth(data_monde,
                    locations = "code",
                    color = "Cluster",
                    projection = "natural earth",
                    color_discrete_sequence = px.colors.qualitative.Set2,
                    hover_name = data_monde.index,
                    hover_data = {"code" : False}
                   )
#fig.update_geos(fitbounds="locations", visible=True)
fig.update_layout(margin = {"r":0,"t":0,"l":0,"b":0})
fig.update_geos(showframe = False)
fig.show()

Terminons en fabriquant une carte régionale.

for reg in pd.unique(data_monde["Région du monde"]):
    print(reg)

Western Europe
North America and ANZ
Middle East and North Africa
Latin America and Caribbean
Central and Eastern Europe
East Asia
Southeast Asia
Commonwealth of Independent States
Sub-Saharan Africa
South Asia

region = data_monde[data_monde["Région du monde"] == "Middle East and North Africa"]

fig = px.choropleth(region,
                    locations = "code",
                    color = "Score de bonheur",
                    projection = "natural earth",
                    color_continuous_scale = "Temps",
                    hover_name = region.index,
                    hover_data = {"code" : False}
                   )
fig.update_geos(fitbounds = "locations", visible = True)
fig.update_layout(margin = {"r":0,"t":0,"l":0,"b":0})
fig.update_geos(showframe = False, resolution = 50)
fig.show()

Modifiez les cellules qui précèdent pour que le graphique ci-dessus affiche la carte du score de générosité des pays d’Asie du sud-est.

De quelle couleur est le Vietnam sur cette carte ?

Fin du TP3c

Série temporelle

Utilisons un nouveau jeu de données comprenant des relevés de consommation électrique allemands entre 2006 et 2018 :

url = "http://cordier-phychi.toile-libre.org/Info/github/elec_allemagne.csv"
serie_temp = pd.read_csv(url,sep=",")
serie_temp.drop(columns="Wind+Solar",inplace=True)
serie_temp

Petit toilettage des données : on transforme les valeurs de la colonne des dates en un type date reconnu par pandas et on les utilise comme index.

serie_temp['Date'] = pd.to_datetime(serie_temp['Date'])
serie_temp = serie_temp.set_index('Date')
serie_temp.head()

On francise ensuite les noms de colonne…

serie_temp.columns = ["Consommation","Vent","Solaire"]
serie_temp.head()

Et enfin, on ajoute des colonnes “jour”, “mois” et “année”.

serie_temp['jour'] = serie_temp.index.day_name()
serie_temp['mois'] = serie_temp.index.month
serie_temp['année'] = serie_temp.index.year
serie_temp["date"] = serie_temp.index
serie_temp["date"] = serie_temp["date"].dt.date # pour aider Colab qui a des soucis avec les dates
serie_temp.head()

px.line(serie_temp[["Consommation","Vent","Solaire"]])

On constate d’importantes variations saisonnières.

zoom = serie_temp[serie_temp['année']==2016]
fig1 = px.line(zoom,'date','Consommation')
fig2 = px.scatter(zoom,'date','Consommation',color='jour')
fig = go.Figure()
fig.add_traces([fig1.data[0],*[fig2.data[i] for i in range(7)]])

Une variabilité hebdomadaire se superpose à la tendance saisonnière.

Grâce à la méthode des dataframe pandas groupby, on peut facilement grouper les donner de manière à obtenir les statistiques qui nous intéressent.
Exemple : trouvons combien d’électricité d’origine éolienne a été produite chaque mois en 2016.

serie_temp[serie_temp['année']==2016].groupby("mois")["Vent"].sum()

mois
1 9264.588
2 9814.294
3 6030.177
4 5910.504
5 6089.484
6 3369.069
7 4651.582
8 4742.343
9 4222.315
10 5585.248
11 8076.232
12 9252.290
Name: Vent, dtype: float64

Sur le modèle précédent, déterminez le jour de la semaine où l’Allemagne a consommé le plus d’électricité en moyenne en 2016 (vous pourrez utilisez la méthode mean à la place de sum).

Traçons une boîte à moustaches de la répartition des 3 variables mois par mois :

fig, axes = plt.subplots(3, 1, figsize=(15, 10), sharex=True)
for var, ax in zip(['Consommation', 'Solaire', 'Vent'], axes):
    sns.boxplot(data=serie_temp, x='mois', y=var, ax=ax)  
    ax.set_ylabel('GWh')
    ax.set_title(var)  
    if ax != axes[-1]:
        ax.set_xlabel('')

On observe que :

• les trois graphes présentent bien une variabilité saisonnière ; la consommation électrique est plus forte en hiver ainsi que la production éolienne (même si l’écart est moins marqué) et la production solaire est beaucoup plus importante en été.
• beaucoup de valeurs se retrouvent à l’extérieur des moustaches supérieures pour la production éolienne, ce qui est probablement dû à des périodes de fort vent.

Regardons maintenant jour par jour :

serie_temp["date"]=(serie_temp.index.strftime('%d %B'))
px.box(serie_temp,x='jour', y='Consommation',hover_data={"date"})

Pourquoi y a-t-il autant de points au-delà des moustaches les jours de semaine ?