Python pour le calcul scientifique/Régression et optimisation
Rappelons que dorénavant les programmes commencent tous par :
#!/usr/bin/python3
import numpy as np
import matplotlib.pyplot as plt
Régression
[modifier | modifier le wikicode]Régression linéaire
[modifier | modifier le wikicode]La fonction permettant la régression linéaire par la méthode des moindres carrés s'appelle est np.linalg.lstsq()
(linear algebra — least squares).
Nous disposons d'un nuage de n points (xi , yi )0 ≤ i ≤ n – 1. Nous voulons les modéliser par une loi affine :
- y = a⋅x + b
les paramètres a et b étant à déterminer. On peut représenter le problème sous la forme d'un système d'équations :
ce qui peut également s'écrire sous la forme d'une équation matricielle
- Y = XA + E
avec
La syntaxe est alors de la forme : resultat = np.linalg.lstsq(X, Y)
. La fonction renvoie une liste de valeurs :
resultat[0]
contient les valeurs (a, b) ;resultat[1]
contient la somme des carrés des résidus.
Par exemple
# **************
# * Paramètres *
# **************
nbpts = 20 # nombre de points
a = 5. # paramètres de la droite
b = 2.
var = 0.5 # variance du bruit
# *************
# * Fonctions *
# *************
def droite(A, x):
return A[0]*x + A[1]
# *******************
# * Nuage de points *
# *******************
x = np.linspace(0, 1, nbpts)
y = droite((a, b), x) + var*np.random.randn(nbpts)
# **************
# * Résolution *
# **************
X = np.vstack([x, np.ones(len(x))]).T
# X = [[x1, 1],
# [x2, 1],
# ...
# [xn, 1]]
resultat = np.linalg.lstsq(X, y)
aopt, bopt = resultat[0]
erreur = resultat[1][0]/nbpts
# ***************************
# * Affichage des résultats *
# ***************************
plt.plot(x, y, "b+") # nuage de points
plt.plot(x, droite((aopt, bopt), x)) # droite de régression
print("a =", aopt.round(2), ", b =", bopt.round(2),
"; χ² =", erreur.round(2))
Si nous voulons que la loi passe par l'origine (b = 0), le système d'équations devient :
soit
- Y = XA + E
avec
Il suffit donc de modifier le code source comme suit :
# **************
# * Résolution *
# **************
X = x.reshape(nbpts, 1)
# X = [[x1],
# [x2],
# ...
# [xn]]
resultat = np.linalg.lstsq(X, y)
[…]
print("a =", aopt.round(2), "; χ² =", erreur.round(2))
- Ressource
- (en) « numpy.linalg.lstsq », sur Numpy and Scipy Documentation (consulté le 17 avril 2019)
Régression polynomiale
[modifier | modifier le wikicode]NumPy propose un module polynomial.polynomial
. Une des manières de décrire le polynôme P(x) = a0 + a1x + a2x2 + … + anxn est d'utiliser un vecteur [a0, a1, …, an]
. Pour déterminer le polynôme de degré n le plus proche d'un jeu de données (x, y), nous disposons de la fonction numpy.polynomial.polynomial.polyfit(x, y, n)
. Pour avoir le résidus, il faut ajouter le paramètre full = True
.
Par exemple, avec cet outil, la régression linéaire consiste à faire une régression par un polynôme de degré 1. En reprenant l'exemple précédent, nous pouvons conserver les sections de code * Paramètres *
, * Fonctions *
, * Nuage de points *
et * Affichage des résultats *
, et nous ajoutons le code suivant :
import numpy.polynomial.polynomial as nppol
# **************
# * Résolution *
# **************
A, param = nppol.polyfit(x, y, 1, full=True)
aopt = A[1]
bopt = A[0]
erreur = param[0]/nbpts
Code complet
import numpy as np
import numpy.polynomial.polynomial as nppol
import matplotlib.pyplot as plt
# **************
# * Paramètres *
# **************
nbpts = 20 # nombre de points
a = 5. # paramètres de la droite
b = 2.
var = 0.5 # variance du bruit
# *************
# * Fonctions *
# *************
def droite(A, x):
return A[0]*x + A[1]
# *******************
# * Nuage de points *
# *******************
x = np.linspace(0, 1, nbpts)
y = droite((a, b), x) + var*np.random.randn(nbpts)
# **************
# * Résolution *
# **************
A, param = nppol.polyfit(x, y, 1, full=True)
aopt = A[1]
bopt = A[0]
erreur = param[0]/nbpts
# ***************************
# * Affichage des résultats *
# ***************************
plt.plot(x, y, "b+") # nuage de points
plt.plot(x, droite((aopt, bopt), x)) # droite de régression
print("a =", aopt.round(2), ", b =", bopt.round(2),
"; χ² =", erreur.round(2))
Si nous voulons avoir une loi qui passe par l'origine, il faut indiquer non pas le degré maximum du polynôme mais la liste des degrés. Comme nous ne voulons que le terme en x1 et pas le terme en x0, nous indiquons donc [1]
. Le code devient donc :
A, param = nppol.polyfit(x, y, [1], full=True)
- Pour plus de détails voir : Python pour le calcul scientifique/Polynômes#Régression polynomiale.
Régression non-linéaire
[modifier | modifier le wikicode]La régression non linéaire nécessite de charger le module optimize
de SciPy. Ce module propose alors plusieurs fonctions.
La première est scipy.optimize.least_squares()
. Cette méthode consiste à minimiser la fonction de résidus, qui est donc une fonction de coûts, par la méthode des moindres carrés ; il faut donc lui donner la fonction de résidus en paramètre. C'est une méthode itérative. Nous définissons donc :
- une fonction
residus(A, x, y)
oùA
est la liste des paramètres de la fonction ; - un jeu de paramètres initial
A0
(point de départ des itérations) ; - un nuage de points (x, y).
La syntaxe est alors resultat = scipy.optimize.least_squares(residus, A0, arg=(x, y))
. La fonction renvoie un dictionnaire contenant entre autres :
resultat.x
: paramètres optimisés Aopt ;resultat.cost
: valeur de la fonction de coûts.
Reprenons l'exemple précédent de la régression linéaire ; nous reprenons les sections * Paramètres *
, * Fonctions *
, * Nuage de points *
et * Affichage des résultats *
, et nous ajoutons le code suivant :
import scipy.optimize as opt
# **************
# * Paramètres *
# **************
A0 = (1, 0) # paramètres initiaux
# *************
# * Fonctions *
# *************
def residus(A, x, y):
return y - droite(A, x)
# **************
# * Résolution *
# **************
resultat = opt.least_squares(residus, A0, args=(x, y))
Aopt = resultat.x # paramètres à l'optimum
aopt = Aopt[0]
bopt = Aopt[1]
erreur = 2*resultat.cost/nbpts
Notez que :
- du point de vue de l'utilisation finale, la fonction que l'on cherche est une fonction de x, définie par un jeu de paramètres A ; on peut la noter ƒA(x) ;
- du point de vue de l'optimisation, on cherche le minimum de la fonction résidus, une fonction de A, définie par un jeu de paramètres (xi, yi)i ∈ [0 ; n – 1] qui sont les points expérimentaux ; on peut l'écrire résidus(xi, yi)i ∈ [0 ; n – 1](A) = ∑i [yi – ƒA(xi)].
Pour cette raison, lorsque l'on définit la fonction de résidus, c'est impérativement A qui vient en premier : residus(A, *args, **kwargs)
.
Code complet
import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize as opt
# **************
# * Paramètres *
# **************
nbpts = 20 # nombre de points
a = 5. # paramètres de la droite
b = 2.
var = 0.5 # variance
A0 = (1, 0) # paramètres initiaux
# *************
# * Fonctions *
# *************
def droite(A, x):
return A[0]*x + A[1]
def residus(A, x, y):
return y - droite(A, x)
# *******************
# * Nuage de points *
# *******************
x = np.linspace(0, 1, nbpts)
y = droite((a, b), x) + var*np.random.randn(nbpts)
# **************
# * Résolution *
# **************
resultat = opt.least_squares(residus, A0, args=(x, y))
Aopt = resultat.x # paramètres à l'optimum
aopt = Aopt[0]
bopt = Aopt[1]
erreur = 2*resultat.cost/nbpts
# ***************************
# * Affichage des résultats *
# ***************************
plt.plot(x, y, "b+") # nuage de points
plt.plot(x, droite((aopt, bopt), x)) # droite de régression
print("a =", aopt.round(2), "; b =", bopt.round(2),
"; χ² =", erreur.round(2))
Si nous voulons que la loi passe par l'origine (b = 0), il suffit de définir une loi purement linéaire :
def droite(x, A):
return A*x
def residus(A, x, y):
return y - droite(x, A)
[…]
Aopt = opt.least_squares(residus, A0, args=(x, y))
Le module optimize
propose également la fonction scipy.optimize.curve_fit()
qui permet d'éviter de définir la fonction de résidus. On lui indique directement la fonction ƒA(x). Les deux principales différences par rapport à .least_squares()
sont :
- lorsque l'on définit la fonction ƒA(x), on indique d'abord x puis ensuite la liste séparée des paramètres :
f(x, a0, a1, …, am)
; - la fonction
.curve_fit()
n'indique pas la somme quadratique des résidus (le khi carré, χ²), il faut la calculer à part.
La syntaxe est Aopt, Acov = scipy.optimize.curve_fit(f, x, y, p0=(a00, a01, …, a0m))
où A0 = (a00, a01, …, a0m) est le jeu de paramètres initial.
La variable de sortie Acov est la matrice de covariance du jeu de paramètres ; elle permet de connaître l'erreur standard sur la détermination des paramètres A.
Reprenons l'exemple précédent de la régression linéaire ; nous reprenons les sections * Paramètres *
, * Nuage de points *
et * Affichage des résultats *
, et nous ajoutons le code suivant :
import scipy.optimize as opt
# **************
# * Paramètres *
# **************
a00 = 1 # paramètres initiaux
a01 = 0
# *************
# * Fonctions *
# *************
def droite(x, a, b):
return a*x + b
# **************
# * Résolution *
# **************
Aopt, Acov = opt.curve_fit(droite, x, y, p0=(a00, a01))
aopt = Aopt[0]
bopt = Aopt[1]
res = y - droite(x, aopt, bopt)
erreur = (res*res).sum()/nbpts
Code complet
import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize as opt
# **************
# * Paramètres *
# **************
nbpts = 20 # nombre de points
a = 5. # paramètres de la droite
b = 2.
var = 0.5 # variance
a00 = 1 # paramètres initiaux
a01 = 0
# *************
# * Fonctions *
# *************
def droite(x, a, b):
return a*x + b
# *******************
# * Nuage de points *
# *******************
x = np.linspace(0, 1, nbpts)
y = droite(x, a, b) + var*np.random.randn(nbpts)
# **************
# * Résolution *
# **************
Aopt, Acov = opt.curve_fit(droite, x, y, p0=(a00, a01))
aopt = Aopt[0]
bopt = Aopt[1]
res = y - droite(x, aopt, bopt)
erreur = (res*res).sum()/nbpts
# ***************************
# * Affichage des résultats *
# ***************************
plt.plot(x, y, "b+") # nuage de points
plt.plot(x, droite(x, aopt, bopt)) # droite de régression
print("a =", aopt.round(2), "; b =", bopt.round(2),
"; χ² =", erreur.round(2))
Là encore, si nous voulons une loi passant par l'origine, le code devient :
def droite(x, A):
return A*x
[…]
Aopt, Acov = opt.curve_fit(droite, x, y, p0=A0)
- Ressources
- (en) « scipy.optimize.least_squares », sur Numpy and Scipy Documentation (consulté le 17 avril 2019)
- (en) « scipy.optimize.curve_fit », sur Numpy and Scipy Documentation (consulté le 17 avril 2019)
Régression sur un espace de dimension supérieure
[modifier | modifier le wikicode]Nous considérons maintenant une loi y = ƒA(x) où x est un vecteur de dimension m et y est un scalaire :
- ƒ : ℝn → ℝ ou bien ℂn → ℂ.
Régression linéaire
[modifier | modifier le wikicode]Le principe est le même, la matrice X est de dimension n × (m + 1) :
- Y = XA + E
avec
Par exemple, si nous considérons un espace de dimension 2 :
- ƒA : ℝ2 → ℝ
- (x, y) ↦ z = ƒA(x, y)
une loi bilinéaire z = ax + by + c décrit un plan, les matrices sont :
et le code devient :
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# **************
# * Paramètres *
# **************
nbpts = 5 # nombre de points par axe
a = 5. # paramètres du plan
b = 3.
c = 2.
sigma = 0.5 # écart type
# *************
# * Fonctions *
# *************
def plan(x, y, A):
return A[0]*x + A[1]*y + A[2]
# *******************
# * Nuage de points *
# *******************
x = np.linspace(0, 1, nbpts)
y = np.linspace(0, 1, nbpts)
grilleX, grilleY = np.meshgrid(x, y)
X = grilleX.flatten().reshape(nbpts*nbpts, 1)
Y = grilleY.flatten().reshape(nbpts*nbpts, 1)
Z = plan(X, Y, (a, b, c)) + sigma*np.random.randn(nbpts*nbpts, 1)
# **************
# * Résolution *
# **************
points = np.hstack((X, Y, np.ones_like(X)))
# points = [[x1, y1, 1]
# [x2, y2, 1]
# ...
# [xn, yn, 1]]
resultat = np.linalg.lstsq(points, Z)
aopt, bopt, copt = resultat[0]
erreur = resultat[1][0]/(nbpts*nbpts)
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection="3d")
ax.set_title("np.linalg.lstsq()")
ax.scatter(X, Y, Z, color="r")
ax.plot_surface(grilleX, grilleY, plan(grilleX, grilleY, (aopt, bopt, copt)), alpha=0.8)
print("""*****************
np.linalg.lstsq()""")
print("a =", aopt.round(2), ", b =", bopt.round(2), "; c =", copt.round(2),
"; χ² =", erreur.round(2))
Optimisation
[modifier | modifier le wikicode]L'optimisation consiste à trouver le minimum d'une fonction ƒ — le minimum ou le maximum, il suffit le cas échéant de changer de signe. Si la fonction ƒ est définie sur un ensemble E, nous cherchons donc
- Xopt ∈ E tel que ∀ X ∈ E, ƒ(X) ≥ ƒ(Xopt).
Notez qu'un tel Xopt n'existe pas forcément, et qu'il n'est pas forcément unique. Par exemple, la fonction ƒ(x) = x ne possède pas de minimum sur ℝ, et la fonction ƒ(x) = constante possède une infinité de minima.
On peut restreindre E, en particulier par des inégalités, par exemple trouver le minimum de ƒ sur ℝ avec x > 0.
Cas général
[modifier | modifier le wikicode]Dans le cas d'une fonction réelle, le module scipy.optimize
propose la fonction minimize_scalar()
. Par défaut, cette fonction utilise la méthode de Brent qui est une méthode itérative. La syntaxe est simplement resultat = optimize.minimise_scalar(f)
.
Le résultat est un objet de la classeOptimizeResult
: c'est un dictionnaire contenant entre autres :
x
: la valeur de x donnant le minimum ;success
: un booléen valantTrue
si la méthode a convergé,False
sinon ;fun
: valeur de la fonction au minimum ;jac
,hes
: valeur du jacobien et du hessien ;nit
: nombre d'itérations.
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
def fonction(x):
'''fonction présentant un minimum proche de zéro'''
return x*x + 5*np.sin(x)
x = np.linspace(-2*np.pi, 2*np.pi, 40)
y = fonction(x)
minimum = optimize.minimize_scalar(fonction)
xopt = minimum.x
ymin = minimum.fun
print(minimum)
print(round(xopt, 2), round(ymin, 2))
plt.plot(x, y)
plt.plot([xopt], [ymin], "+k")
Pour les fonction de plusieurs variables, nous disposons de la fonction minimize()
. Par défaut, cette fonction utilise la méthode BFGS (méthode de Broyden-Fletcher-Goldfarb-Shanno) qui est également une méthode itérative. Pour cette fonction, il faut indiquer un point de départ x0
:
optimize.minimize(f, x0)
Comme il s'agit d'un algorithme de descente, l'algorithme peut être piégé dans un minimum local.
Avec l'exemple précédent :
code source
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
def fonction(x):
'''fonction présentant un minimum proche de zéro'''
return x*x + 5*np.sin(x)
x = np.linspace(-2*np.pi, 2*np.pi, 40)
y = fonction(x)
minimum = optimize.minimize(fonction, 0)
xopt = minimum.x[0]
ymin = minimum.fun
print(minimum)
print(round(xopt, 2), round(ymin, 2))
plt.plot(x, y)
plt.plot([xopt], [ymin], "+k")
Notes et références
[modifier | modifier le wikicode]Interpolation, extrapolation et lissage < ↑ > Algèbre linéaire