L'objectif de ce TP est de mettre en pratique le langage python pour réaliser une regression linéaire. L'idée est, dans un premier temps, de reprendre les éléments de base du langage (condition, boucles ...) pour créer un outil qui calcule les paramètres par la méthode des moindres carrés. Dans une deuxième partie les modules spécifiques tels que numpy ou matplotlib seront mis à profit pour réaliser la même opération avec des méthodes existantes.
Note : Ce notebook est compatible python2 et python3. Il est recommandé d'utiliser python3.
Le programme que nous allons écrire devra réaliser les opérations suivantes :
donnees.dat.La régression linéaire consiste à chercher les paramètres $a$ et $b$ définissant la droite $y=ax+b$ qui passe au plus près d'un ensemble de points $(x_k,y_k)$. Les paramètres $a$ et $b$ sont déterminés par la méthodes des moindres carrés qui consiste, dans le cas d'une régression linéaire, à minimiser la quantité :
\begin{equation} Q(a, b) = \sum_{k=1}^N (y_k - a x_k - b)^2 \end{equation}Le minimum de $Q(a,b)$ est obtenu lorsque ses dérivées par rapport à $a$ et $b$ sont nulles. Il faut donc résoudre le système à deux équations deux inconnues suivant :
\begin{align} & \begin{cases} \displaystyle\frac{\partial Q(a,b)}{\partial a} = 0 \\ \displaystyle\frac{\partial Q(a,b)}{\partial b} = 0 \end{cases} & \Leftrightarrow & & & \begin{cases} \displaystyle -2 \sum_{k=1}^N x_k \left(y_k - a x_k - b\right) = 0 \\ \displaystyle -2 \sum_{k=1}^N \left(y_k - a x_k - b\right) = 0 \end{cases} \end{align}Les solutions de ce système sont :
\begin{align} a & = \frac{\displaystyle N \sum_{k=1}^N x_k y_k - \sum_{k=1}^N x_k\sum_{k=1}^N y_k}{\displaystyle N\sum_{k=1}^N x_k^2 - \left(\sum_{k=1}^N x_k\right)^2} & b & = \frac{\displaystyle \sum_{k=1}^N x_k^2 \sum_{k=1}^N y_k - \sum_{k=1}^N x_k\sum_{k=1}^N x_k y_k}{\displaystyle N\sum_{k=1}^N x_k^2 - \left(\sum_{k=1}^N x_k\right)^2} \end{align}Le programme sera écrit de plusieurs façon différentes.
In [1]:
cat data/donnees.dat
Lecture du fichier ligne par ligne :
In [2]:
with open("data/donnees.dat", "r") as inp:
for line in inp:
xi, yi = line.split()
print(xi, yi)
print("type de xi : ", type(xi))
On va maintenant enregistrer les valeurs xi et yi dans des listes. Par ailleurs, on peut remarquer que les valeurs sont lues comme des chaînes de caractères. Il faut les convertir en nombre flottants avec la fonction float().
In [3]:
# création des listes
x = list()
y = list()
# lecture du fichier
with open("data/donnees.dat", "r") as inp:
for line in inp:
xi, yi = line.split()
x.append(float(xi))
y.append(float(yi))
Nous avons maintenant nos listes de valeurs de x et y :
In [4]:
print(x)
In [5]:
print(y)
Nous verrons par la suite comment lire ce type de fichier de façon efficace avec la méthode loadtxt() de numpy.
In [6]:
# initialisation
x_sum = 0
# calcul de la somme
for xi in x:
x_sum += xi
# affichage
print("somme des valeurs de x = ", x_sum)
On fait de même pour la somme des valeurs de $y$ et des valeurs de $x$ au carré.
In [7]:
y_sum = 0.
for yi in y:
y_sum += yi
print("somme des valeurs de y = ", y_sum)
In [8]:
x2_sum = 0.
for xi in x:
x2_sum += xi**2
print("somme des valeurs de x^2 = ", x2_sum)
Il reste à calculer la somme des produits $x\times y$. Pour parcourir à la fois la liste des valeurs de $x$ et de $y$ on utilise la fonction zip qui joint les deux listes :
In [9]:
for xi, yi in zip(x, y):
print("xi = ", xi, "\tyi = ", yi)
Mettons cela à profit pour calculer la somme :
In [10]:
xy_sum = 0.
for xi, yi in zip(x, y):
xy_sum += xi * yi
print("somme des valeurs de x*y = ", xy_sum)
Maintenant que nous disposons de toutes les valeurs nécessaires, il ne reste plus qu'à calculer $a$ et $b$. Nous avons encore besoin du nombre de points. La fonction len donne le nombre d'éléments d'une liste.
In [11]:
npoints = len(x)
print("Nombre de points = ", npoints)
D'après les équations présentées en introduction :
In [12]:
a = (npoints * xy_sum - x_sum * y_sum) / (npoints * x2_sum - x_sum**2)
print("a = ", a)
In [13]:
b = (x2_sum * y_sum - x_sum * xy_sum) / (npoints * x2_sum - x_sum**2)
print("b = ", b)
Dans cette fonction nous allons regrouper les différentes étapes permettant de réaliser la régression linéaire. La fonction prend comme arguments les listes des valeurs de $x$ et $y$ et retourne les valeurs des paramètres $a$ et $b$. Les entrées et sorties de la fonction doivent être explicités dans la docstring située en dessous de la définition.
In [14]:
def regLin(x, y):
"""
Ajuste une droite d'équation a*x + b sur les points (x, y) par la méthode
des moindres carrés.
Args :
* x (list): valeurs de x
* y (list): valeurs de y
Return:
* a (float): pente de la droite
* b (float): ordonnée à l'origine
"""
# initialisation des sommes
x_sum = 0.
x2_sum = 0.
y_sum = 0.
xy_sum = 0.
# calcul des sommes
for xi, yi in zip(x, y):
x_sum += xi
x2_sum += xi**2
y_sum += yi
xy_sum += xi * yi
# nombre de points
npoints = len(x)
# calcul des parametres
a = (npoints * xy_sum - x_sum * y_sum) / (npoints * x2_sum - x_sum**2)
b = (x2_sum * y_sum - x_sum * xy_sum) / (npoints * x2_sum - x_sum**2)
# renvoie des parametres
return a, b
Utilisons maintenant cette nouvelle fonction.
In [15]:
a, b = regLin(x, y)
print("a = ", a)
print("b = ", b)
Pour afficher les nombres flotant, il est possible d'utiliser un format. Deux syntaxes existent suivant la version de python :
In [16]:
# python 2.7 et superieur
print("a = {:8.3f}".format(a))
print("b = {:8.3f}".format(b))
In [17]:
# python 2.X
print("a = %8.3f" % a)
print("b = %8.3f" % b)
In [18]:
import numpy as np
La somme des valeurs d'un tableau peut être obtenue par la méthode sum().
In [19]:
a = np.array([1, 2, 3])
a.sum()
Out[19]:
La somme des carrés ou des produits peut également être obtenue aisément.
In [20]:
(a**2).sum()
Out[20]:
Pour calculer les produits entre deux arrays numpy il suffit de les multiplier :
In [21]:
b = np.array([3, 2, 1])
print("a : ", a)
print("b : ", b)
print("a * b :", a * b)
La somme se calcule alors de la même manière que précédemment pour le carré :
In [22]:
(a * b).sum()
Out[22]:
Nous pouvons maintenant simplifier la fonction regLin en utilisant les fonctions de numpy :
In [23]:
def regLin_np(x, y):
"""
Ajuste une droite d'équation a*x + b sur les points (x, y) par la méthode
des moindres carrés.
Args :
* x (list): valeurs de x
* y (list): valeurs de y
Return:
* a (float): pente de la droite
* b (float): ordonnée à l'origine
"""
# conversion en array numpy
x = np.array(x)
y = np.array(y)
# nombre de points
npoints = len(x)
# calculs des parametres a et b
a = (npoints * (x*y).sum() - x.sum()*y.sum()) / (npoints*(x**2).sum() - (x.sum())**2)
b = ((x**2).sum()*y.sum() - x.sum() * (x*y).sum()) / (npoints * (x**2).sum() - (x.sum())**2)
# renvoie des parametres
return a, b
La nouvelle fonction renvoie évidemment les mêmes résultats :
In [24]:
a, b = regLin_np(x, y)
print("a = {:8.3f}\nb = {:8.3f}".format(a, b)) # \n est le caractere de fin de ligne
Numpy et Scipy sont deux modules scientifiques de python qui regroupent de nombreuses fonctions. Nous allons utiliser la méthode loadtxt pour lire le fichier texte et les méthodes polyfit et linregress pour réaliser la régression linéaire. Le module numpy est totalement inclu dans scipy. Compte tenu du grand nombre de modules et bibliothèques python existants, il est important de savoir lire une documentation pour utiliser les méthodes disponibles. De plus, l'utilisation d'une méthode déjà existante accélère le travail de développement tout en évitant de refaire la même chose. Pour obtenir la documentation à l'intérieur de ipyhon, ajouter un ? après le nom de la méthode.
Commençons par lire le fichier donnees.dat avec la méthode loadtxt pour lire le fichier de données.
In [25]:
data = np.loadtxt("data/donnees.dat")
print(data)
Les valeurs de $x$ correspondent à la première colonne.
In [26]:
x = data[:,0]
print(x)
Les valeurs de $y$ correspondent à la deuxième colonne.
In [27]:
y = data[:,1]
print(y)
Le tout peut être fait en une seul ligne :
In [28]:
x, y = np.loadtxt("data/donnees.dat", unpack=True)
print(x)
print(y)
In [29]:
parametres = np.polyfit(x, y, 1)
print(parametres)
Les paramètres sont bien les mêmes que ceux que nous avions déterminés à la main.
In [30]:
a, b = parametres
print("a = {:8.3f}\nb = {:8.3f}".format(a, b))
In [31]:
from scipy.stats import linregress
Utilisons maintenant la méthode linregress.
In [32]:
a, b, r, p_value, std_err = linregress(x, y)
print("a ={:8.3f}\nb ={:8.3f}\nr^2 ={:8.5f}".format(a, b, r**2))
Python offre la possibilité de réaliser une réprésentation graphique en utilisant le module matplotlib. Voici un exemple d'utilisation de ce module pour représenter les points $(x, y)$ et la droite de régression précédemment déterminée.
Pour une utilisation dans ipython il faut d'abord préparer l'environnement pour intégrer matplotlib.
In [33]:
%matplotlib inline
Chargeons le module matplotlib.
In [34]:
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = 8, 8 # ajuste la taille des figures
Préparons maintenant le graphique.
In [35]:
plt.plot(x, y, "bo", label="donnees") # les points (x, y) representes par des points
plt.plot( # droite de regression
[x[0], x[-1]], # valeurs de x
[a * x[0] + b, a * x[-1] + b], # valeurs de y
"r-", # couleur rouge avec un trait continu
label="regression") # legende
plt.xlabel("x") # nom de l'axe x
plt.ylabel("y") # nom de l'axe y
plt.xlim(0, 11) # échelle axe x
plt.legend() # la legende
plt.title("Regression Lineaire") # titre de graphique
Out[35]:
Nous allons maintenant utiliser les méthodes vues précédemment sur un échantillon plus large. De plus, nous représenterons les résidus qui correspondent à la différence entre les points issus des données et la droite obtenue par la régression linéaire.
Commençons par lire les données dans le fichier regLinData.dat
In [36]:
x, y = np.loadtxt("data/regLinData.dat", unpack=True)
Utilisons maintenant la fonction linregress pour trouver les paramètres de la régression linéaire :
In [37]:
a, b, r, p_vall, std_err = linregress(x, y)
print("a = ", a, " b = ", b, " r^2 = ", r**2)
Nous allons maintenant visualiser les points et la droite de régression dans un premier graphique, puis, les résidus dans un second graphique. Pour ce faire, nous utiliserons la méthode subplot de matplotlib qui place des graphiques sur une grille. Les trois arguments de subplot sont le nombre de lignes, le nombre de colonnes et le numéro du graphique.
In [38]:
# graphique 1: les données et la droite de régression
plt.subplot(2, 1, 1)
plt.plot(x, y, "bo", label="data")
plt.plot(x, a * x + b, "r-", label="regression")
plt.xlabel("x")
plt.ylabel("y")
plt.title("Regression linéaire avec résidus")
plt.legend(loc="lower right")
# graphique 2: les résidus
plt.subplot(2, 1, 2)
plt.plot(x, y - (a * x + b), "g-", label="residus")
plt.xlabel("x")
plt.ylabel("Résidus")
plt.legend(loc="upper center")
Out[38]:
Les résidus sont un bon indicateur de la qualité du modèle. Ici, on dit que les résidus sont 'structurés'. Ils ne sont pas aléatoirement répartis autour de zéro, ils présentent une variation parabolique. Cette structure des résidus indique qu'une fonction affine n'est pas adaptée pour représenter les données.
Utilisons la méthode polyfit pour ajuster une parabolle :
In [39]:
p0, p1, p2 = np.polyfit(x, y, 2)
print("p0 = ", p0, " p1 = ", p1, " p2 = ", p2)
Reprennons les lignes précédentes pour représenter les données graphiquement.
In [40]:
# graphique 1: les données et la droite de régression
plt.subplot(2, 1, 1)
plt.plot(x, y, "bo", label="data")
plt.plot(x, p0*x**2 + p1*x + p2, "r-", label="parabolle")
plt.xlabel("x")
plt.ylabel("y")
plt.title("Regression linéaire avec résidus")
plt.legend(loc="lower right")
# graphique 2: les résidus
plt.subplot(2, 1, 2)
plt.plot(x, y - (p0*x**2 + p1*x + p2), "g-", label="residus")
plt.xlabel("x")
plt.ylabel("Résidus")
plt.legend(loc="upper center")
Out[40]:
Cette fois, les résidus sont bien répartis aléatoirement autour de l'origine. Calculons le coefficient de détermination selon \begin{equation} R^2 = \frac{\sum_k (y^{calc}_k - \overline{y})^2}{\sum_k (y_k - \overline{y})^2} \end{equation} où les $y_k$ sont les données, les $y^{calc}$ sont ceux calculés avec le polynômes et $\overline{y}$ est la moyenne des valeurs $y_k$.
In [41]:
R2 = ((p0*x**2 + p1*x + p2 - y.mean())**2).sum() / ((y - y.mean())**2).sum()
print(R2)
Comme attendu le coefficient de détermination est meilleur que celui obtenu pour la droite (0.9916). Pour information, la fonction utilisée pour produire les données du fichier regLinData.dat est :
\begin{equation}
f(x) = 0.08 x^2 + 2x + 3 + \hat{g}(X)
\end{equation}
où $\hat{g}(X)$ est un bruit gaussien obtenu par le tirage de nombres aléatoires dans une distribution gaussienne centrée en zéro et de largeur 2.