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 ("xi", "yi"), et l'on détermine une fonction décrivant ces points. Les différences majeures sont :
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
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 :
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)" :
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 :
Comme il n'y a que deux paramètres dans le modèle, on note ceux-ci "a" et "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) :
ce qui peut s'écrire sous forme matricielle
avec
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 codice_1.
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 codice_2, avec la syntaxe :
// définition du nuage de points
X = […];
Y = […];
// régression
[a, b, sigma] = reglin(X, Y);
On peut aussi utiliser la commande codice_3. 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 codice_4, 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 codice_5 :
chdir("monrepertoire/");
// Constantes
// Fonctions
deff("[y] = affine(a, b, x)", "y = a*x + b")
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]);
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 codice_2 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 :
alors il suffit de définir X comme étant une matrice "m"×"n"
 X = [X11, X12, … X1"n" ;
 X21, X22, … X2"n" ;
 X"m"1, X"m"2, … X"mn"]
où X"ij" est la valeur de "xj" pour le point de mesure "i" ; et de définir Y comme étant un vecteur ligne de "n" valeurs
 Y = [Y1, Y2, … , Y"n"]
et le paramètre "a" retourné par la division matricielle ou par la fonction codice_2 est le vecteur ["a"1, "a"2, … , "a""n"].
Régression polynomiale.
La régression polynomiale est un cas particulier de régression multilinéaire.
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"3 - 2·"x"2 + 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 codice_16 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 codice_16 retourne donc un vecteur :
 Exemple 1
Nous reprenons l'exemple de la régression linéaire, mais en utilisant cette fois-ci la fonction codice_16. 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