« Découvrir Scilab/Régression » : différence entre les versions

import de Calcul numérique
(ébauche)
 
(import de Calcul numérique)
<span style="font-size:25px;">6. Régression</span>
----
 
'''Présentation du problème'''
 
Le problème a des similarités avec celui de l'interpolation : on a des points expérimentaux (''x<sub>i</sub>'', ''y<sub>i</sub>''), 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<sub>1</sub> ; X<sub>2</sub> ; … ; X<sub>''n''</sub>] ;
* Y = [Y<sub>1</sub> ; Y<sub>2</sub> ; … ; Y<sub>''n''</sub>].
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))<sup>2</sup>/A(3)) : une fonction gaussienne ;
* ƒ(A, ''x'') = A(1)·''x''<sup>3</sup> + A(2)·''x''<sup>2</sup> + A(3)·''x'' + A(4) : un polynôme de degré 3.
La régression consiste à touver 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 = ∑<sub>''i''</sub> ''r''<sup>2</sup>(A, X<sub>''i'', Y<sub>''i''</sub>) ;
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) :
: <math>\begin{cases}
a x_1 + b = y_1 \\
a x_2 + b = y_2 \\
\ldots \\
a x_n + b = y_n \\
\end{cases}</math>
ce qui peut s'écrire sous forme matricielle
: A×X = Y
avec
: A = (''a'' ''b'')
: <math>\mathrm{X} = \begin{pmatrix}
x_1 & x_2 & \ldots & x_n \\
1 & 1 & \ldots & 1 \\
\end{pmatrix}</math>
: Y = (''y''<sub>1</sub> ''y''<sub>2</sub> … ''y<sub>n&thinsp;</sub>'')
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 :
<source lang="scilab">
\\ 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
</source>
 
Si l'on veut forcer l'ordonnée à l'origine à 0, il suffit de ne pas mettre la ligne de 1 :
<source lang="scilab">
A = y/x
</source>
 
On peut aussi travailler avec des vecteurs colonne, et utiliser l'autre diviseur matriciel <code>A = x\y</code>.
 
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ésidus.
 
Pour cela, on peut utiliser la fonction <code>reglin()</code>, avec la syntaxe :
 
<source lang = "scilab">
// définition du nuage de points
X = […];
Y = […];
 
// régression
[a, b, sigma] = reglin(X, Y);
</source>
 
{{note
| Les matrices X et Y doivent être des vecteurs ''ligne''.
}}
 
On peut aussi utiliser la commande <code>regress()</code>. 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).
 
<source lang = "scilab">
coefs = regress(X, Y)
</source>
 
'''Exemple'''
 
Dans l'exemple suivant, on génère des point dispersés aléatoirement autour d'une droite <code>y = pente*x + ordonnee</code>, 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 <code>monrepertoire/</code> :
* le fichier <code>fonction.sce</code> qui contient la fonction affine ;
* le fichier <code>generation.sce</code> qui crée un nuage de points et le met dans le fichier <code>donnees.txt</code> ;
* le fichier <code>regression.sce</code> qui lit le fichier <code>donnees.txt</code> et effectue la régression.
 
: Fichier <code>fonction.sce</code>
<source lang="scilab">
chdir("monrepertoire/");
 
// **********
// Constantes
// **********
 
// **********
// Fonctions
// **********
 
deff("[y] = affine(a, b, x)", "y = a*x + b")
</source>
 
: Fichier <code>generation.sce</code>
<source lang="scilab">
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]);
</source>
 
: Fichier <code>regression.sce</code>
<source lang="scilab">
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))
</source>
 
Comme la forme de la fonction ƒ est connue, il est inutile de la définir. La fonction <code>reglin()</code> 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''<sub>1</sub>, ''a''<sub>2</sub>, …, ''a''<sub>''m''</sub>, ''b'', ''x''<sub>1</sub>, ''x''<sub>2</sub>, …''x''<sub>''n''</sub>) = ''a''<sub>1</sub>·''x''<sub>1</sub> + ''a''<sub>2</sub>·''x''<sub>2</sub> + … + ''a''<sub>''n''</sub>·''x''<sub>''n''</sub> + ''b''
alors il suffit de définir X comme étant une matrice ''m''×''n''
X = [X<sub>11</sub>, X<sub>12</sub>, … X<sub>1''n''</sub> ;
X<sub>21</sub>, X<sub>22</sub>, … X<sub>2''n''</sub> ;
X<sub>''m''1</sub>, X<sub>''m''2</sub>, … X<sub>''mn''</sub>]
où X<sub>''ij''</sub> est la valeur de ''x<sub>j</sub>'' pour le point de mesure ''i'' ; et de définir Y comme étant un vecteur ligne de ''n'' valeurs
Y = [Y<sub>1</sub>, Y<sub>2</sub>, … , Y<sub>''n''</sub>]
et le paramètre ''a'' retourné par la division matricielle ou par la fonction <code>reglin()</code> est le vecteur [''a''<sub>1</sub>, ''a''<sub>2</sub>, … , ''a''<sub>''n''</sub>].
 
== Régression polynomiale ==
 
La régression polynomiale est un cas particulier de régression multilinéaire.
 
{{loupe|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''<sup>3</sup> - 2·''x''<sup>2</sup> + 3·''x'' - 4, bruité :
 
<source lang="scilab">
// 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']);
</source>
 
Et voici un programme permettant d'effectuer la régression sur ces données :
 
<source lang="scilab">
// 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, "+")
</source>
 
== 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 <code>leastsq()</code> avec la syntaxe suivante :
<source lang="scilab">
// 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)
</source>
Avec cette syntaxe, la fonction <code>leastsq()</code> retourne donc un vecteur :
* le premier élément, <code>S</code>, est la valeur de la somme des carrés des écarts ;
* le second élément, <code>Aopt</code>, est le vecteur contenant la valeur « optimisée » des paramètres du modèle.
 
{{note
| Les matrices X et Y doivent être ici des vecteurs ''colonne''. Par ailleurs, la fonction de calcul du résidus 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 <code>leastsq()</code>. Lorsque l'on enregistre les données générées, on utilise :
<source lang="scilab">
write("donnees.txt", [X', Y']);
</source>
afin d'avoir des vecteurs colonne.
 
<source lang="scilab">
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))
</source>
 
[[File:Desommation gaussienne lorentzienne.svg|vignette|upright=2|Résultat de la régression.]]
 
''' 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).
 
<source lang="scilab">
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))))
</source>
 
L'algorithme par défaut de <code>leastsq()</code> est l'algorithme de pseudo-Newton. On peut aussi utiliser l'[[w:Méthode du gradient conjugué|algorithme du gradient conjugué]] en ajoutant <code>"gc"</code> aux paramètres, ou bien algorithme foncitonnant avec les fonctions non-différentiables avec <code>"nd"</code>. On peut enfin utiliser la commande <code>lsqrsolve()</code>, qui utilise l'algorithme de Levenberg-Marquardt.
 
== Notes ==
13 633

modifications