Rappelons que dorénavant les programmes commencent tous par :
import numpy as np
import matplotlib.pyplot as plt
Régression.
Régression linéaire.
La fonction permettant la régression linéaire par la méthode des moindres carrés s'appelle est codice_1 "()".
Nous disposons d'un nuage de "n" points ("xi ", "yi ")0 ≤ "i" ≤ "n" – 1. Nous voulons les modéliser par une loi affine :
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
avec
La syntaxe est alors de la forme : codice_2. La fonction renvoie une liste de valeurs :
Par exemple
nbpts = 20 # nombre de points
a = 5. # paramètres de la droite
b = 2.
var = 0.5 # variance du bruit
def droite(A, x):
 return A[0]*x + A[1]
x = np.linspace(0, 1, nbpts)
y = droite((a, b), x) + var*np.random.randn(nbpts)
X = np.vstack([x, np.ones(len(x))]).T
resultat = np.linalg.lstsq(X, y)
aopt, bopt = resultat[0]
erreur = resultat[1][0]/nbpts
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
avec
Il suffit donc de modifier le code source comme suit :
X = x.reshape(nbpts, 1)
resultat = np.linalg.lstsq(X, y)
print("a =", aopt.round(2), "; χ² =", erreur.round(2))
Régression polynomiale.
NumPy propose un module codice_5. Une des manières de décrire le polynôme P("x") = "a"0 + "a"1"x" + "a"2"x"2 + … + "anxn" est d'utiliser un vecteur codice_6. 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 codice_7. Pour avoir le résidus, il faut ajouter le paramètre codice_8.
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 codice_9, codice_10, codice_11 et codice_12, et nous ajoutons le code suivant : 
import numpy.polynomial.polynomial as nppol
A, param = nppol.polyfit(x, y, 1, full=True)
aopt = A[1]
bopt = A[0]
erreur = param[0]/nbpts
import numpy as np
import numpy.polynomial.polynomial as nppol
import matplotlib.pyplot as plt
nbpts = 20 # nombre de points
a = 5. # paramètres de la droite
b = 2.
var = 0.5 # variance du bruit
def droite(A, x):
 return A[0]*x + A[1]
x = np.linspace(0, 1, nbpts)
y = droite((a, b), x) + var*np.random.randn(nbpts)
A, param = nppol.polyfit(x, y, 1, full=True)
aopt = A[1]
bopt = A[0]
erreur = param[0]/nbpts
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 "x"1 et pas le terme en "x"0, nous indiquons donc codice_13. Le code devient donc :
A, param = nppol.polyfit(x, y, [1], full=True)
Régression non-linéaire.
La régression non linéaire nécessite de charger le module codice_14 de SciPy. Ce module propose alors plusieurs fonctions.
La première est codice_15. 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 :
La syntaxe est alors codice_19. La fonction renvoie un dictionnaire contenant entre autres :
Reprenons l'exemple précédent de la régression linéaire ; nous reprenons les sections codice_9, codice_10, codice_11 et codice_12, et nous ajoutons le code suivant :
import scipy.optimize as opt
A0 = (1, 0) # paramètres initiaux
def residus(A, x, y):
 return y - droite(A, x)
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 :
Pour cette raison, lorsque l'on définit la fonction de résidus, c'est impérativement A qui vient en premier : codice_26.
import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize as opt
nbpts = 20 # nombre de points
a = 5. # paramètres de la droite
b = 2.
var = 0.5 # variance
A0 = (1, 0) # paramètres initiaux
def droite(A, x):
 return A[0]*x + A[1]
def residus(A, x, y):
 return y - droite(A, x)
x = np.linspace(0, 1, nbpts)
y = droite((a, b), x) + var*np.random.randn(nbpts)
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
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 codice_14 propose également la fonction codice_28 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 à codice_29 sont :
La syntaxe est codice_32 où A0 = ("a"00, "a"01, …, "a"0"m") 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 codice_9, codice_11 et codice_12, et nous ajoutons le code suivant :
import scipy.optimize as opt
a00 = 1 # paramètres initiaux
a01 = 0
def droite(x, a, b):
 return a*x + b
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
import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize as opt
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
def droite(x, a, b):
 return a*x + b
x = np.linspace(0, 1, nbpts)
y = droite(x, a, b) + var*np.random.randn(nbpts)
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
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)
Régression sur un espace de dimension supérieure.
Nous considérons maintenant une loi "y" = ƒA("x") où "x" est un vecteur de dimension "m" et "y" est un scalaire :
Régression linéaire.
Le principe est le même, la matrice X est de dimension "n" × ("m" + 1) :
avec
Par exemple, si nous considérons un espace de dimension 2 :
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
nbpts = 5 # nombre de points par axe
a = 5. # paramètres du plan
b = 3.
c = 2.
sigma = 0.5 # écart type
def plan(x, y, A):
 return A[0]*x + A[1]*y + A[2]
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)
points = np.hstack((X, Y, np.ones_like(X)))
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.
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
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.
Dans le cas d'une fonction réelle, le module codice_36 propose la fonction codice_37. Par défaut, cette fonction utilise la méthode de Brent qui est une méthode itérative. La syntaxe est simplement codice_38.
Le résultat est un objet de la classecodice_39 : c'est un dictionnaire contenant entre autres :
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 codice_48. 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 codice_49 :
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 :
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.
Interpolation, extrapolation et lissage < ↑ > Algèbre linéaire