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

Régression linéaire

modifier

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 ≤ in – 1. Nous voulons les modéliser par une loi affine :

y = ax + 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

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

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

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)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).

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

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

Régression sur un espace de dimension supérieure

modifier

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

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

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

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 valant True 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 :

Notes et références

modifier

Interpolation, extrapolation et lissage < > Algèbre linéaire