Régression linéaire

L'objectif est de tester la validité d'une loi physique par le tracé d'un nuage de points et sa modélisation par une droite.

Le document propose plusieurs script :

  • Régression linéaire SANS barres d'erreur,
  • Régression linéaire AVEC barres d'erreur,
  • Estimation de l'incertitude-type sur les paramètres de la régression linéaire.

Bibliothèques utilisées


La bibliothèque numpy est utilisée pour réaliser des calculs et éventuellement simuler un processus aléatoire (module random).

La bibliothèque matplolib permet de tracer des graphiques.

Régression linéaire SANS barres d'erreur


Les listes V1 et V2 contiennent les valeurs expérimentales.

On introduit X et Y, grandeurs calculées à partir des valeurs des listes V1 et V2, entre lesquelles on cherche à tester une relation affine.

In [2]:
#Importation des bibliothèques
import numpy as np
import matplotlib.pyplot as plt   

#Valeurs expérimentales
V1=[2.45 , 4.14 , 8.5 , 17.6 , 36.9 , 96.0]
V2=[3.56 , 5.85 , 11.5 , 22.3 , 43.5 , 104]

#Calcul des ordonnées et abscisses 
X = [(V2[i] - V1[i])/1E6 for i in range(len(V1))]   
Y = [(V2[i]-V1[i])/V1[i] for i in range(len(V1))]  

#Modélisation 
p=np.polyfit(X,Y,1)

#Tracé du nuage de points et de la droite modèle
plt.figure(1)
plt.plot(X,Y,'o')
plt.plot(X,np.polyval(p,X))
plt.grid()
plt.show()

#Calcul et affichage des résidus
plt.figure(2)
plt.plot(X , Y-np.polyval(p,X),'ro')
plt.title("Résidus")
plt.axhline()
plt.grid()
plt.show()

#Affichage
print(f"Coefficient directeur : {p[0]}")
print(f"Ordonnée à l'origine : {p[1]}")
Coefficient directeur : -51789.24405177459
Ordonnée à l'origine : 0.5082053777164175

Régression linéaire avec barres d'erreur


Le script est quasiment le même que le précédent, à la différence qu'il nécessite de renseigner des barres d'erreur pour les valeurs de X et de Y.

Cette méthode est plus riche : la loi est validée si les barres d'erreur chevauchent la droite de regression.

In [6]:
#Importation des bibliothèques
import numpy as np
import matplotlib.pyplot as plt   

#Valeurs expérimentales
V1=[2.45 , 4.14 , 8.5 , 17.6 , 36.9 , 96.0]
V2=[3.56 , 5.85 , 11.5 , 22.3 , 43.5 , 104]

#Calcul des ordonnées et abscisses 
X = [(V2[i] - V1[i])/1E6 for i in range(len(V1))]   
Y = [(V2[i]-V1[i])/V1[i] for i in range(len(V1))]  

#Saisie de l'amplitude des erreurs sur X et sur Y 
BEX = 1E-7    # valeur par défaut : 0
BEY = 6E-3    # valeur par défaut : 0

#Modélisation 
p=np.polyfit(X,Y,1)   #renvoie une liste [coef. dir., ord. origine]

#Tracé du nuage de points et de la droite modèle
plt.figure(1)
plt.errorbar(X , Y , xerr=BEX , yerr=BEY , fmt='.')
plt.plot(X,np.polyval(p,X))
plt.grid()
plt.show()

#Calcul et affichage des résidus
plt.figure(2)
plt.errorbar(X, Y-np.polyval(p,X) , xerr=BEX , yerr=BEY , fmt='ro')
plt.title("Résidus")
plt.axhline()
plt.grid()
plt.show()

#Affichage
print(f"Coefficient directeur : {p[0]}")
print(f"Ordonnée à l'origine : {p[1]}")
Coefficient directeur : -51789.24405177459
Ordonnée à l'origine : 0.5082053777164175

Obtention de l'incertitude-type sur les paramètres de la régression linéaire


La méthode consiste à simuler la répétition de l'expérience en proposant un tirage au sort de jeux de points de coordonnées (X,Y). La procédure de tirage au sort nécessite de founir un intervalle au sein duquel les valeurs sont tirées.

Par défaut, la procédure utilise la fonction np.random.normal(valeur centrale,écart-type) ce qui suppose de fournir un écart-type sur X ou sur Y.

Deux listes, A et B, accueillent les valeurs simulées du coefficient directeur et de l'ordonnée à l'origine.

L'incertitude-type sur ces deus grandeurs est l'écart-type dee l'ensemble des valeurs des listes A et B.

In [ ]:
# Nombre de simulations MC
N = 100000

# Listes des valeurs du coefficient directeur et de l'ordonnée à l'origine
A = []
B = []

# Saisie des incertitudes-type sur X et sur Y
u_X = 3E-7
u_Y = 6E-3

# Réalisation de la simulation
for n in range(0,N):
    Xsim,Ysim=[],[]
    for i in range(len(X)) :
      Xsim.append(np.random.normal(X[i],u_X))
      Ysim.append(np.random.normal(Y[i],u_X))
    p = np.polyfit(Xsim,Ysim,1)
    A.append(p[0])
    B.append(p[1])

# Calcul des incertitudes-type
u_A=np.std(A,ddof=1)
u_B=np.std(B,ddof=1)

# Affichage des résultats
print(f"Incertitude-type sur le coefficient directeur : {u_A}")
print(f"Incertitude-type sur l'ordonnée à l'origine : {u_B}")
Incertitude-type sur le coefficient directeur : 2512.7521947600785
Incertitude-type sur l'ordonnée à l'origine : 0.012231112009902435