Description

Travaux pratiques 3 : R´egression

Le but de ce TP est de traiter du probl`eme suivant : ´etant donn´es n points (xi,yi) dans un plan (les donn´ees), comment trouver une fonction ou une courbe f dont le graphe approxime bien de ces points?
En ´equations, on a

ou` i est l’erreur que l’on fait en approximant yi par f(xi). On aimerait d´eterminer une fonction f qui minimise l’erreur totale.
Ainsi pos´e, le probl`eme n’est pas bien d´efini. En effet il existe une infinit´e de courbes passant exactement par ces n points, et ayant par cons´equent une erreur nulle. Il faut donc donner d’autres informations concernant la courbe f : par exemple d´efinir ce que l’on entend par ”approximer” et aussi pr´eciser les propri´et´es de f. Typiquement, on aimerait choisir f parmi une famille de fonctions (par exemple des fonctions polynomiales de degr´e donn´e, p´eriodiques, gaussiennes, etc…). On appelle ce probl`eme r´egression ou fitting. Matlab fournit un certain nombre d’outils puissants pour l’analyse de donn´ees, que nous allons introduire dans ce TP.
1 Importer des donn´ees
Matlab permet d’importer des donn´ees contenues dans des fichiers de donn´ees (*.dat). Pour ceci, cliquez sur Home → Import Data et choisissez le fichier de donn´ees. Alternativement, vous pouvez choisir le r´epertoire contenant les donn´ees dans la fenˆetre Current folder et double cliquer sur le fichier correspondant. Un assistant s’ouvre alors, qui permet de cr´eer dans le Workspace des vecteurs ou des matrices contenant ces donn´ees.
Testez l’importation de donn´ees en t´el´echargeant le fichier ”data1.dat” sur Moodle et en l’important dans une variable matricielle de taille n × 2 appell´ee data1, ou` n est le nombre de mesures. Notez que pour importer le fichier de donn´ees comme variable matricielle, vous devez choisir l’option Numeric matrix au lieu de l’option par d´efaut Column vectors.
2 R´egression lin´eaire : la m´ethode des moindres carr´es
2.1 Moindre carr´es
D´efinissons une quantit´e mesurant dans quelle mesure une fonction approche bien notre ensemble de points. Etant donn´e les´ n points (xi,yi), on va chercher une fonction f qui minimise n
Ω(f) = X(yi − f(xi))2 ≥ 0, (1)
i=1
c’est-`a-dire la somme des carr´es des distances entre les donn´ees et les valeurs de f. Il est clair que si aucune sp´ecification n’est requise sur f, alors il existe une infinit´e de fonctions annulant Ω pour un nombre fini de points donn´es. Il est donc indispensable de d´efinir le type de fonctions approximant nos donn´ees.
2.2 R´egression par droites affines
Consid´erons d’abord le cas ou` l’on demande que la courbe f soit une droite, f(x) = ax+b, a,b ∈ R. Dans ce cas, on peut faire ce calcul explicitement. On r´e´ecrit d’abord (1)
. (2)
Pour minimiser la fonction Ω(a,b) on doit trouver les param`etres a et b tels que
. (3)
En calculant ces deriv´ees explicitement, on obtient le syst`eme d’´equations suivant :
. (4)
Exercice 1 R´egression affine. Commencez par v´erifier la formule (4).
1. Ecrivez une fonction RegressionAffine prenant comme param`etres d’entr´ee une matrice n×2 contenant les points (xi,yi) et retournant les param`etres a et b r´esolvant le syst`eme d’´equations (4).
2. Cette fonction doit dessiner sur un mˆeme graphe la droite de r´egression d´etermin´eepar (a,b) ainsi que les donn´ees de d´epart.
Indications :
— N’oubliez pas d’inclure la commande hold on pour pouvoir dessiner plusieurs objets sur le mˆeme graphique.
— Pour dessiner les donn´ees, fournissez l’argument suppl´ementaire ’ob’ `a la commande plot, de mani`ere `a dessiner un cercle bleu pour chaque point de donn´ee.
— Pour dessiner la droite, utilisez l’argument suppl´ementaire ’-r’ de mani`ere `a lier les points par une ligne rouge.
3. Une fois la fonction programm´ee, t´el´echargez les fichiers de donn´ees ”data1.dat”,”data2.dat”, ”data3.dat”, ”data4.dat” et ”data5.dat” sur Moodle, importez-les dans Matlab et utilisez RegressionAffine sur ces donn´ees dans votre script. Que constatezvous en examinant vos graphiques?
4. Pour d´eterminer le minimum global d’une fonction, il n’est en g´en´eral pas suffisant decalculer simplement les points critiques (c’est-`a-dire un point ou` la d´eriv´ee s’annule). Un point critique peut aussi correspondre `a un maximum local ou global, `a un minimum local, `a un point d’inflexion, ou `a un point-selle. Dans notre cas, pourquoi n’est-il pas n´ecessaire d’ˆetre plus prudent et de v´erifier que le point critique que l’on a trouv´e est bien un minimum?
5. Il existe une fa¸con plus simple d’implementer la fonction RegressionAffine. En effet, Matlab impl´emente la m´ethode des moindre carr´es pour la r´egression lin´eaire par l’op´eration [], d´ej`a utilis´ee au TP1. On rappelle que X Y est une matrice A telle que XA = Y . Dans le cas d’une r´egression, ce syst`eme n’est pas r´esoluble exactement (il n’y a que deux variables ind´ependentes, a et b, et une contrainte pour chaque paire (xi,yi)). Au lieu de retourner une erreur, Matlab retourne la meilleure approximation selon la m´ethode des moindre carr´es.
Impl´ementez une fonction RegressionAffine2 faisant le mˆeme travail que RegressionAffine, mais utilisant []. V´erifiez que vous obtenez les mˆemes droites de r´egression sur les donn´ees.
Indication : Une petite difficult´e est qu’il faut transformer notre r´egression affine en un probl`eme lin´eaire. Ceci est possible en ajoutant une colonne `a la matrice X ne comportant que des 1s. Quelle forme a la matrice A?
2.3 R´egression par cercles
Il arrive que des probl`emes apparemment non lin´eaires puissent ˆetre n´eanmoins lin´earis´es, comme on va le voir dans l’exemple suivant.
Dans les fichiers ”data4.dat” et ”data5.dat”, on a vu des dispositions de points ressemblant plus `a des ellipses ou `a des cercles qu’`a une droite. Notre prochain but est de d´evelopper une m´ethode lin´eaire de r´egression pour trouver un cercle approximant au mieux un nuage de points.
Un cercle est d´etermin´e par son centre de coordonn´ees (a,b) et son rayon r. Un point P de coordonn´ees (x,y) est sur le cercle si seulement si (x − a)2 + (y − b)2 = r2, ou de mani`ere
´equivalente, si
Φ(a,b,r,x,y) := ((x − a)2 + (y − b)2 − r2)2 (5)
s’annule. On remarque que la fonction Φ est positive, diff´erentiable et qu’elle s’annule sur les points du cercle, il est donc naturel d’essayer de la minimiser par rapport aux param`etres
(a,b,r).
Si on a n points de coordonn´ees (xi,yi) pour i variant de 1 `a n, on pose
) (6)
(7)
(8)
. (9)
L’astuce : On pose u = a2+b2−r2 et on consid`ere une nouvelle fonction en les trois variables ind´ependantes a, b et u d´efinie par
n
Ω(a,b,u) = X(x2i − 2xi a + yi2 − 2yi b + u)2 (10)
i=1
On d´etermine a, b et u en minimisant cette fonction, et obtient le syst`eme lin´eaire
, , = 0 (11)
que l’on r´esout comme pr´ec´edemment.
Exercice 2 R´egression circulaire.
1. Ecrivez une fonction´ RegressionCirculaire qui prend comme param`etres d’entr´ee une matrice n × 2 contenant les points (xi,yi) et retourne les param`etres a, b et r d´ecrivant le cercle.
2. Faˆıtes en sorte que RegressionCirculaire dessine sur un mˆeme graphe le cercle de r´egression ainsi que les donn´ees de d´epart.
3. Testez RegressionCirculaire dans votre script sur les donn´ees contenues dans les fichiers ”data4.dat” et ”data5.dat”.
Lorsque l’on demande que la fonction f soit un polynˆome en x et y de degr´e donn´e, on peut proc´eder de mani`ere similaire aux deux exemples pr´ec´edents. Par exemple, au lieu d’utiliser des cercles comme dans l’exemple simple ci-dessus, il est parfois judicieux d’utiliser des coniques
(ellipse, parabole, hyperbole), dont l’´equation est donn´ee par un polynˆome g´en´eral de degr´e 2 :
Ax2 + Bxy + Cy2 + Dx + Ey + F = 0 ,A,B,C,D,E,F ∈ R . (12)
3 R´egression non-lin´eaire : nlinfit
Pour les cas ou` la fonction f dans (1) est quelconque, Matlab fournit un outil de r´egression puissant : nlinfit. Il s’agit de trouver parmi une famille de fonctions {fa}, ou` a est un param`etre, celle approchant le mieux les donn´ees. On va impl´ementer le cas ou` la famille {fa} est une famille de courbes gaussiennes caract´eris´ees par leur moyenne, ´ecart-type et amplitude.
3.1 R´egression gaussiennes sur des donn´ees
Si l’on effectue un grand nombre de mesures d’une certaine quantit´e et si les mesures sont ind´ependantes les unes des autres, alors ces valeurs se distribuent autour de la valeur moyenne selon la loi normale ou gaussienne, dont la fonction de r´epartition est donn´ee par
. (13)
µ est la moyenne de la gaussienne, σ est l’´ecart-type (correspondant `a l’”´etalement” de la gaussienne), et A := P(µ) est l’amplitude (voir figure 3.1).
Supposons que l’on ait notre ensemble de mesures {mk}k=1,…,N, donn´e par un ensemble de N points distribu´es le long de l’axe des x. On va commencer par construire un histogramme de la fa¸con suivante. On choisit d’abord un intervalle I assez grand pour contenir toutes nos donn´ees, c’est-`a-dire tous les points {mk}. On d´ecoupe I en n intervalles (”bins”) {Ii}i=1,…,n. On note xi pour le centre de l’intervalle Ii et on compte le nombre yi de points dans chaque bin Ii. On obtient ainsi un ensemble de points {(xi,yi)}i=1,…,n dans le plan, sur lequel on peut appliquer une r´egression.
Comme on s’attend `a ce que la distribution de nos points suive une loi normale, on va appliquer une r´egression par rapport `a une famille de courbes gaussiennes
(14)
parametr´ee par A, µ et σ. C’est typiquement un cas ou` l’on peut utiliser la fonction nlinfit.

Figure 1 – Quelques courbes gaussiennes, dont les amplitudes ont ´et´e normalis´ees `a A =
.
Exercice 3 R´egression gaussienne. Le fichier ”datag.dat” contient environ 104 mesures prises au LEP de la masse de la particule W. T´el´echargez-le et importez-le dans Matlab.
1. A partir de ces donn´ees, construisez un histogramme h en utilisant la fonction histogram d´ecrite dans la Section 3.2 ci-dessous.
2. histogram nous fournit un vecteur h.BinEdges contenant les points de s´eparation des bins. Construisez un vecteur CentresBins contenant le centre de chaque bin. Plottez le vecteur de valeurs de l’histogramme, h.Values, en fonction de CentresBins et comparez `a l’histogramme pour v´erifier que tout fonctionne.
3. Cr´eez une fonction FamGauss qui admet comme variables un vecteur tridimensionel a et une variable r´eelle x. FamGauss doit retourner la valeur au point x de la gaussienne d’amplitude a(1), de moyenne a(2) et d’´ecart-type a(3). (Comme pr´ec´edemment, FamGauss doit ˆetre d´efinie dans votre script, pas dans un fichier s´epar´e.)
4. Effectuez une r´egression par rapport `a la famille de fonctions gaussiennes FamGauss sur le nuage de points produit par l’histogramme. Utilisez pour ceci la fonction nlinfit d´ecrite dans la Section 3.3. Donnez l’amplitude, la moyenne et l’´ecart-type de la courbe de r´egression.
5. Plottez cette gaussienne et v´erifiez qu’elle approxime bien l’histogramme. Pour unr´esultat esth´etique, ajoutez les trois param`etres suppl´ementaires ’k’, ’Linewidth’, 2 `a la fonction plot (pour dessiner la courbe gaussienne en noir avec une largeur de trait de deux points).
6. R´ep´etez l’op´eration avec un nombre de bins diff´erent. Est-ce que les param`etres de lacourbe de r´egression ont chang´e? Si oui, lesquels et pourquoi?
3.2 La fonction histogram
Matlab fournit la fonction histogram, qui permet de construire un histogramme `a partir de donn´ees. histogram s’utilise tr`es simplement de la mani`ere suivante :
h = histogram(m)
ou` m est le vecteur contenant les donn´ees.
La variable retourn´ee h est une variable de type histogram, qui contient toutes les informations `a propos de l’histogramme. Par exemple, le nombre de mesures dans chaque bins est contenu dans le vecteur h.Values (de taille n). Les s´eparations entre les bins sont contenues dans un vecteur h.BinEdges (de taille n+1). Le nombre, la taille et le placement des bins est effectu´e automatiquement, mais peut ˆetre ajust´e manuellement si n´ecessaire. Par exemple, on peut cr´eer un histogramme avec 50 bins par la commande h = histogram(m,50)
Pour plus d’informations, cf. help histogram.
3.3 Utilisation de nlinfit
La proc´edure pour passer la famille de fonctions {fa} comme un param`etre dans nlinfit est la suivante.
1. On cr´ee une nouvelle fonction Matlab (par exemple FuncFit, enregistr´ee dans le fichier FuncFit.m), qui va encoder notre famille de fonctions. FuncFit prend comme entr´ee le param`etre a param´etrisant la famille de fonctions fa et la variable x (dans cet ordre!). Elle retourne fa(x) en utilisant la valeur de a pass´ee en param`etre. En d’autre termes, on consid`ere la famille de fonctions param´etr´ee par a comme une unique fonction d´ependant de a et de x. a et x peuvent bien suˆr ˆetre des variables vectorielles. Par exemple, si la famille de fonctions `a minimiser est {fA,B,ω|A,B,ω ∈ R} et
fA,B,ω(x) = Asin(ωx) + B cos(ωx), (15)
on ´ecrit une fonction qui prend comme entr´ee le vecteur de param`etres (A,B,ω) ainsi que la variable x, et qui retourne fA,B,ω(x).
2. Dans le script principal, on utilise la fonction nlinfit de la fac¸on suivante :
pf = nlinfit(x,y,@FuncFit,p0).
— (x,y) sont deux vecteurs contenant les deux coordonn´ees des points sur lesquels on aimerait effectuer la r´egression.
— FuncFit est la famille de fonctions par rapport `a laquelle on aimerait effectuer la r´egression, encod´ee comme expliqu´e plus haut.
— p0 contient une valeur de d´epart des param`etres, sur lequel l’algorithme se base pour trouver la valeur optimale des param`etres. Il est n´ecessaire de fournir ces valeurs initiales. Id´ealement, le r´esultat devrait ˆetre ind´ependant de celles-ci, dans la pratique, ¸ca n’est pas toujours le cas…
— La variable retourn´ee pf contient les valeurs des param`etres qui minimisent (1).
Pour plus d’informations, voir help nlinfit.