Découvrir Scilab/Régression

6. Régression

Présentation du problème

Le problème a des similarités avec celui de l'interpolation : on a des points expérimentaux (x_i, y_i), et l'on détermine une fonction décrivant ces points. Les différences majeures sont :

la fonction est déterminée sur l'ensemble de définition complet, à partir de tous les points, et non pas par morceaux entre deux points ;
la fonctions ne passe en général pas par les points expérimentaux mais près de ces points.

Par ailleurs, la forme de la fonction est souvent déterminée par le type de problème auquel on s'attaque, elle a un « sens physique », là où l'interpolation se contente de prendre des fonctions simples.

Considérons un nuage de n points « expérimentaux » (X, Y) : X et Y sont des vecteurs

X = [X₁ ; X₂ ; … ; X_n] ;
Y = [Y₁ ; Y₂ ; … ; Y_n].

Considérons une fonction « modèle » y = ƒ(A, x) où A est un vecteur ; les valeurs de A sont des paramètres de la fonction ƒ, par exemple :

ƒ(A, x) = A(1)·exp((x - A(2))²/A(3)) : une fonction gaussienne ;
ƒ(A, x) = A(1)·x³ + A(2)·x² + A(3)·x + A(4) : un polynôme de degré 3.

La régression consiste à trouver le vecteur A tel que la fonction ƒ « colle le mieux » au nuage de points. On utilise pour cela la méthode des moindres carrés (least square) :

on définit la fonction « résidus » r(A, x, y) entre le modèle et les points expérimentaux :
r(A, x, y) = y - ƒ(A, x) ;
on définit la somme S des carrés des résidus comme étant :
S = ∑_i r²(A, X_i, Y_i) ;

la régression par la méthode des moindres carrés consiste donc à trouver le vecteur A donnant la valeur minimale de S.

Régression linéaire

Régression linéaire simple

Dans le cas de la régression linéaire, la fonction modèle ƒ est une fonction affine :

ƒ(A, x) = A(1)·x + A(2).

Comme il n'y a que deux paramètres dans le modèle, on note ceux-ci a et b :

ƒ(a, b, x) = a·x + b.

D'une certaine manière, la régression consiste à résoudre un système d'équations surdéterminé (chaque point expérimental étant une équation, il y a plus d'équations que d'inconnues) :

{\begin{cases}ax_{1}+b=y_{1}\\ax_{2}+b=y_{2}\\\ldots \\ax_{n}+b=y_{n}\\\end{cases}}

ce qui peut s'écrire sous forme matricielle

A×X = Y

avec

A = (a b)

\mathrm {X} ={\begin{pmatrix}x_{1}&x_{2}&\ldots &x_{n}\\1&1&\ldots &1\\\end{pmatrix}}

Y = (y₁ y₂ … y_n)

Si l'on demande à Scilab de résoudre le système avec la division matricielle alors que la matrice X n'est pas carrée, il effectue automatiquement une régression linéaire des moindres carrés :

\\ génération des données de test
x = 1:10
y = x + 0.3*rand(x);

\\ régression
X = [x ; ones(x)];
A = y/X

Si l'on veut forcer l'ordonnée à l'origine à 0, il suffit de ne pas mettre la ligne de 1 :

A = y/x

On peut aussi travailler avec des vecteurs colonne, et utiliser l'autre diviseur matriciel A = x\y.

Mais la méthode de la division matricielle ne donne pas d'autre information que les coefficients, et en particulier ne donne pas le résidu.

Pour cela, on peut utiliser la fonction reglin(), avec la syntaxe :

// définition du nuage de points
X = […];
Y = […];

// régression
[a, b, sigma] = reglin(X, Y);

Les matrices X et Y doivent être des vecteurs ligne.

On peut aussi utiliser la commande regress(). Elle ne permet pas d'obtenir l'écart type, et ne fonctionne qu'avec deux variables, mais par contre, les vecteurs peuvent être indifféremment ligne ou colonne ; ils peuvent même être de types différents (un vecteur ligne et l'autre colonne).

coefs = regress(X, Y)

Exemple

Dans l'exemple suivant, on génère des points dispersés aléatoirement autour d'une droite y = pente*x + ordonnee, puis on fait la régression linéaire de ce nuage de points.

Pour clarifier le code source, le programme est scindé en trois fichiers situés dans le même répertoire (dossier) appelé ici monrepertoire/ :

le fichier fonction.sce qui contient la fonction affine ;
le fichier generation.sce qui crée un nuage de points et le met dans le fichier donnees.txt ;
le fichier regression.sce qui lit le fichier donnees.txt et effectue la régression.

Fichier fonction.sce

chdir("monrepertoire/");

// **********
// Constantes
// **********

// **********
// Fonctions
// **********

deff("[y] = affine(a, b, x)", "y = a*x + b")

Fichier generation.sce

chdir("monrepertoire/");

// **********
// Constantes
// **********

// paramètres de la loi affine
pente = 1;
ordonnee = 10;

// paramètres de la loi normale
var = 1; // variance

// nombre de points
nbpts = 20;

// **********
// Initialisation
// **********

X = zeros(nbpts)';
Y = zeros(nbpts)';

// **********
// Fonctions
// **********

exec("fonction.sce", -1) // charge le fichier (-1 : sans commentaire)

// **********
// Programme principal
// **********

// nuage de points
for i = 1:nbpts
    X(1,i) = i + var*rand(1, "normal");
    Y(1,i) = affine(pente, ordonnee, i) + var*rand(1, "normal");
end

// enregistrement

write("donnees.txt", [X,Y]);

Fichier regression.sce

clf;
chdir("monrepertoire/");

// **********
// Programme principal
// **********

// lecture des données
donnees = read("donnees.txt", -1, 2)
Xexp = donnees(:,1);
Yexp = donnees(:,2);


// regression
[aopt, bopt, sigma] = reglin(Xexp, Yexp)

// loi modèle
Yopt = affine(aopt, bopt, X);

// affichage
plot(Xexp, Yexp, "+b")
plot(Xexp, Yopt, "-r")
xstring(0, ordonnee, "a = "+string(aopt)+" pour "+string(pente)...
+" ; b = "+string(bopt)+" pour "+string(ordonnee)...
+" ; sigma = "+string(sigma))

Comme la forme de la fonction ƒ est connue, il est inutile de la définir. La fonction reglin() retourne un vecteur de trois composantes : les paramètres a et b du modèle, et l'écart type σ du résidu (qui est la racine carrée de la somme des carrés des résidus, σ = √S).

Régression multilinéaire

Notons que a, b, x et y peuvent être des vecteurs. Si par exemple on a un modèle linéaire à plusieurs variables :

ƒ(a₁, a₂, …, a_m, b, x₁, x₂, …x_n) = a₁·x₁ + a₂·x₂ + … + a_n·x_n + b

alors il suffit de définir X comme étant une matrice m×n

X = [X₁₁, X₁₂, … X_1n ;
     X₂₁, X₂₂, … X_2n ;
     …
     X_m1, X_m2, … X_mn]

où X_ij est la valeur de x_j pour le point de mesure i ; et de définir Y comme étant un vecteur ligne de n valeurs

Y = [Y₁, Y₂, … , Y_n]

et le paramètre a retourné par la division matricielle ou par la fonction reglin() est le vecteur [a₁, a₂, … , a_n].

Régression polynomiale

La régression polynomiale est un cas particulier de régression multilinéaire.

Pour plus de détails voir : wikipedia:fr:Régression polynomiale.

Voici par exemple un programme générant n points pour x allant de 0 à 5 (valeurs exactes), et y suivant le polynôme x³ - 2·x² + 3·x - 4, bruité :

// n points aléatoires approchant un polynôme de degré 3

// chdir("mon_chemin_d_acces"); // répertoire où enregistrer le fichier

n = 10;
x = linspace(0, 5, n);
P = %s^3 - 2*%s^2 - 3*%s + 4;

y = horner(P, x); // valeurs exactes du polynôme

clf;
plot(x, y); // tracé du polynôme

y = y + 3*rand(y, "normal"); // valeurs bruitées

plot(x, y, "+"); // tracé des points bruités

write("polynome_bruite.txt", [x' , y']);

Et voici un programme permettant d'effectuer la régression sur ces données :

// Régression polynomiale de degré 3

// chdir("mon_chemin_d_acces"); // répertoire où chercher le fichier


M = read("polynome_bruite.txt", -1, 2); // lecture des données

x = M(:, 1);
y = M(:, 2);

X = [x.^3, x.^2, x, ones(x)];

A = X\y; // régression linéaire par les moindres carrés

P = [%s^3, %s^2, %s, 1]*A; // définition du polynôme

yth = horner(P, x); // valeurs correspondant au polynôme

// tracé du résultat
clf;
plot(x, yth)
plot(x, y, "+")

Régression non-linéaire

Dans le cas de la régression non-linéaire, ƒ est une fonction quelconque. Pour effectuer la régression, on utilise la fonction leastsq() avec la syntaxe suivante :

// fonction modèle
function y = f(A, x)
   ...
endfunction 

// fonction résidus
function e = r(A, x, y)
   e = f(A, x) - y
endfunction

// initialisation des paramètres du modèle
Ainit = […];

// définition du nuage de points
X = […];
Y = […];

[S, Aopt] = leastsq(list(r, X, Y), Ainit)

Avec cette syntaxe, la fonction leastsq() retourne donc un vecteur :

le premier élément, S, est la valeur de la somme des carrés des écarts ;
le second élément, Aopt, est le vecteur contenant la valeur « optimisée » des paramètres du modèle.

Les matrices X et Y doivent être ici des vecteurs colonne. Par ailleurs, la fonction de calcul du résidu doit s'appliquer aux vecteurs X et Y, pour renvoyer un vecteur e.

Exemple 1

Nous reprenons l'exemple de la régression linéaire, mais en utilisant cette fois-ci la fonction leastsq(). Lorsque l'on enregistre les données générées, on utilise :

write("donnees.txt", [X', Y']);

afin d'avoir des vecteurs colonne.

clf;
chdir("monrepertoire/");

// **********
// Initialisation
// **********

Ainit=[1;1];

// **********
// Fonctions
// **********

deff("[y] = f(a, b, x)", "y = a*x + b")

deff("[e] = res(A, x, y)", "e = f(A(1),A(2),x) - y")

// **********
// Programme principal
// **********
 
// lecture des données
donnees = read("donnees.txt", -1, 2)
Xexp = donnees(:, 1);
Yexp = donnees(:, 2);
 
// regression
[S, Aopt] = leastsq(list(res, Xexp, Yexp),Ainit)
 
// loi modèle
Ymod = f(Aopt(1), Aopt(2), Xexp);
 
// affichage
plot(Xexp, Ymod, "-r")
plot(Xexp, Yexp, "+b")
xstring(0, ordonnee, "a = "+string(Aopt(1))+" pour "+string(pente)...
+" ; b = "+string(Aopt(2))+" pour "+string(ordonnee)...
+" ; variance = "+string(S))

Exemple 2

Nous travaillons ici avec un pic généré par la somme d'une fonction gaussienne et d'une fonction lorentzienne. Ce genre de profil est classiquement utilisé pour décrire la forme des pics obtenus par diffractométrie X. On fait ce que l'on appelle une désommation de pics (on parle souvent de « déconvolution », mais ce terme est abusif puisqu'il ne s'agit pas de produit de convolution).

clf()

// **********
// Constantes
// **********

paramlorentz(1) = 5; // hauteur de la courbe lorentzienne
paramlorentz(2) = 2; // largeur de la courbe lorentzienne
paramgauss(1) = 10; // hauteur de la courbe gaussienne
paramgauss(2) = 3; // largeur de la courbe gaussienne
var=0.5; // variance de la loi normale
nbpts = 100 // nombre de points
demielargeur = max(3*paramgauss(2), 3*paramlorentz(2)) // pour intervalle x
pas = 2*demielargeur/nbpts;
Ainit = [1;1;1;1]; // paramètres initiaux du modèle

// **********
// Initialisation
// **********

X = zeros(nbpts,1);
Y = zeros(nbpts,1);
Yopt = zeros(nbpts,1);

// **********
// Fonctions
// **********

// Fonction gaussienne centrée sur 0

function [y] = gauss(A, x)
    // A(1) : hauteur de pic
    // A(2) : "largeur" du pic
    y = A(1)*exp(-x.*x/A(2));
endfunction

// Fonction lorentzienne centrée sur 0

function [y] = lorentz(A, x)
    // A(1) : hauteur de pic
    // A(2) : "largeur" du pic
    foo = A(2)*A(2)
    y = foo*A(1)./(foo + 4*x.*x);
endfunction

function [y] = profil(A, x)
    // A(1) : hauteur de pic lorentzien
    // A(2) : "largeur" du pic lorentzien
    // A(3) : hauteur de pic gaussien
    // A(4) : "largeur" du pic gaussien
    L = A(1:2);
    G = A(3:4);
    y = lorentz(L, x) + gauss(G, x);
endfunction

// Résidus

function [e] = residus(A, x, y)
    e = profil(A, x) - y;

endfunction

// **********
// Nuage de points
// **********

i=(1:nbpts)';
X = pas*i - demielargeur;
Y = profil([paramlorentz;paramgauss], x) + var*rand(X, "normal");

// **********
// Programme principal
// **********

[S, Aopt] = leastsq(list(residus, X, Y), Ainit)
Yopt = profil(Aopt, X);
YLopt = lorentz(Aopt(1:2),X);
YGopt = gauss(Aopt(3:4),X);

// affichage

plot(X, Yopt, "-r")
plot(X, YLopt, "-c")
plot(X, YGopt, "-g")
plot(X, Y, "+b")

hauteur = paramlorentz(1)+paramgauss(1);

xstring(-demielargeur, hauteur,...
 "lorentzienne : Al = "+string(0.01*round(100*Aopt(1)))...
+" ; Hl = "+string(0.01*round(100*Aopt(2))))
xstring(-demielargeur, 3*hauteur/4,...
 "gaussienne : Ag = "+string(0.01*round(100*Aopt(3)))...
+" ; Hg = "+string(0.01*round(100*Aopt(4))))
xstring(-demielargeur, hauteur/2,...
 "variance : S = "+string(0.01*round(100*S))))

L'algorithme par défaut de leastsq() est l'algorithme de pseudo-Newton. On peut aussi utiliser l'algorithme du gradient conjugué en ajoutant "gc" aux paramètres, ou bien algorithme foncitonnant avec les fonctions non-différentiables avec "nd". On peut enfin utiliser la commande lsqrsolve(), qui utilise l'algorithme de Levenberg-Marquardt.

Notes

Interpolation, extrapolation et lissage < ↑ > Calcul différentiel et intégral