Régressions polynômiales (et autres) avec Python

Comment obtenir des courbes de régression non linéaires à partir de données expérimentales ?

Voici plusieurs modules Python pour tracer ces courbes, obtenir leurs paramètres, ainsi que le coefficient de corrélation R2 afin de comparer différents modèles de régression.

Régressions polynômiales avec le module Numpy

import matplotlib.pyplot as plt
import numpy as np
from sklearn.metrics import r2_score#pour calculer le coeff R2

fig=plt.figure(figsize=(7, 9), tight_layout=True)
plt.suptitle("Exemples de régressions de polynômes")
ax1=plt.subplot(211)
plt.title("polynôme degré 3")
#1er exemple ---------------------------------------------------------------------------
x = [1,2,3,5,6,7,8,9,10,12,13,14,15,16,18,19,21,22]
y = [100,90,80,60,60,55,60,65,70,70,75,76,78,79,90,99,99,100]

mymodel = np.poly1d(np.polyfit(x, y, 3))
myline = np.linspace(1, 22, 100)

plt.scatter(x, y)
plt.plot(myline, mymodel(myline))
plt.text(12, 60, mymodel)
plt.text(12, 55, '{:.5f}'.format(r2_score(y, mymodel(x))))

ax2=plt.subplot(212)
plt.title("polynôme degré 2")
#2e exemple :
x=[0, 
0.1, 
0.2, 
0.3, 
0.4, 
0.5, 
0.6, 
0.7, 
0.8, 
0.9, 
1, 
1.1, 
1.2, 
1.3, 
1.4, 
]

y=[0, 
10.92, 
23.66, 
40.04, 
61.88, 
85.54, 
112.84, 
143.78, 
178.36, 
214.76, 
252.98, 
293.02, 
336.7, 
380.38, 
425.88, 
]

P=np.polyfit(x, y, 2)
#P est une liste des valeurs a,b et c du polynome ax**2+bx+c=0
monmodele=np.poly1d(P)
plt.scatter(x, y)
plt.plot(x, monmodele(x))
#la valeur de R2 se récupère avec r2_score(y,monmodele(x))
plt.text(0.5, 0.5, '{:.5f}'.format(r2_score(y, monmodele(x))))
plt.text(0.4, 300, monmodele)
#------------------------------------------------------------------------
fig.savefig("RegressionPolynomes2-3.png", dpi=250)
plt.show() 

Ici, l’outil r2_score du module sklearn est appelé pour calculer le coefficient de corrélation.

Le premier exemple (polynôme de degré 3) est issu du forum python stackoverflow.

Régression exponentielle

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

def func(x, a, b, c):
    return a * np.exp(-b * x) + c
#cette fonction pourrait être remplacée par une autre. Adapter le nombre de paramètres a, b, c,...

x=np.linspace(0, 10, 11, endpoint=True)

#Données y collées depuis tableur geogebra :
y=[200.000000000000, 
157.000000000000, 
137.000000000000, 
115.000000000000, 
104.000000000000, 
85.0000000000000, 
71.0000000000000, 
59.0000000000000, 
50.0000000000000, 
42.0000000000000, 
37.0000000000000, 
]

popt, pcov = curve_fit(func, x, y)
#popt contient les coeff a, b et c en liste
plt.figure()
plt.plot(x, y, 'ko', label="valeurs expérimentales")
plt.plot(x, func(x, *popt), 'r-', label="Modélisation")
plt.text(2, 50, r"N=%s$\times \exp^{-%s*t}$"%(round(popt[0], 2), round(popt[1], 3)))
plt.legend()
plt.show()

Ici, on utilise l’outil curve_fit du module scipy. Pour obtenir le coefficient R2, on peut procéder comme dans l’exemple précédent, avec r2_score.

Le code précédent peut être adapté à une fonction différente de l’exponentielle, comme dans l’exemple suivant :

import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from sklearn.metrics import r2_score#pour calculer le coeff R2

def func(x, a, b, c):
#    return a*x**2+b*x+c#polynome degré 2
    return a/(x+b)**c
#    return a * np.exp(-b * x) + c
#cette fonction pourrait être remplacée une autre expression

#Données expérimentales :
x=[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 20]
y=[1123, 622, 390, 255, 177, 128, 97, 77, 64, 52, 45, 39, 34, 31, 29, 21]

popt, pcov = curve_fit(func, x, y)
ymodele=func(x, *popt)
R2=r2_score(y, ymodele)#coeff.de corrélation comparant y et ymodele
#popt contient les coeff a, b et c en liste
fig, ax=plt.subplots()
plt.plot(x, y, 'ko', label="valeurs expérimentales")
plt.plot(x, ymodele, 'r-', label="Modélisation")
plt.text(0.95, 0.5, r"$y=\dfrac{%.3e}{(x+%s)^{%s}}$"%(round(popt[0], 2), round(popt[1], 3), round(popt[2], 3)), ha='right', transform=ax.transAxes)
plt.text(0.95, 0.4, r"$R^2=%.7f$"%(R2), ha='right', transform=ax.transAxes)
plt.legend()
fig.savefig("Regression1.png")
plt.show()

Les données de cet exemple sont issues du livre de Laurent Mathieu, Comment j’ai pesé la Terre avec un chronomètre (Ellipses).