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

: 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

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 <code>A = x\y</code>.

Pour cela, on peut utiliser la fonction <code>reglin()</code>, avec la syntaxe :

 

// définition du nuage de points

X = […];

// régression

[a, b, sigma] = reglin(X, Y);

 

{{note

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).

 

coefs = regress(X, Y)

 

'''Exemple'''

 

: Fichier <code>fonction.sce</code>

chdir("monrepertoire/");

 

 

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

 

: Fichier <code>generation.sce</code>

chdir("monrepertoire/");

 

 

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

 

: Fichier <code>regression.sce</code>

clf;

chdir("monrepertoire/");

+" ; 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 <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).

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

 

 

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

 

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 <code>leastsq()</code> avec la syntaxe suivante :

// fonction modèle

function y = f(A, x)

 

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

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 ;

 

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 :

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

afin d'avoir des vecteurs colonne.

 

clf;

chdir("monrepertoire/");

+" ; b = "+string(Aopt(2))+" pour "+string(ordonnee)...

+" ; variance = "+string(S))

 

[[File:Desommation gaussienne lorentzienne.svg|vignette|upright=2|Résultat de la régression.]]

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()

 

xstring(-demielargeur, hauteur/2,...

"variance : S = "+string(0.01*round(100*S))))

 

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.