Découvrir Scilab/Calcul numérique

Un livre de Wikilivres.
Sauter à la navigation Sauter à la recherche

Table des matièresIndex



3. Calcul numérique


Vocabulaire[modifier | modifier le wikicode]

Dans le vocabulaire habituel de la programmation, une procédure est un sous-programme, c'est-à-dire un bout de programme auquel on donne un nom, et que l'on exécute à la demande en « appelant » ce nom. On peut voir cela comme une macro-définition (ou simplement « macro ») : le nom de la procédure est « remplacé » par le contenu du sous-programme. Une fonction est une procédure qui renvoit une valeur.

Dans la documentation officielle de Scilab, les termes function, procedure et macro sont équivalents ; par contre, ils sont utilisés pour les sous-programme définis par l'utilisateur. Les fonctions intégrées dans Scilab, dites built-in, sont appelées « primitives ».

Lorsqu'une primitive a pour argument une fonction, cette fonction est dite « externe » (external). C'est le cas par exemple des primitives de résolution d'équations différentielle ode() ou d'optimisation optim() : on indique à ces primitives quelle est la fonction à utiliser.

Par soucis de simplicité, nous ne feront pas la différence entre les fonctions définies par l'utilisateur et les primitives, sauf lorsque cette nuance est pertinente. Les primitives sont parfois appelées « commandes » ou « instructions » dans le présent document.

Définition de fonctions extérieures[modifier | modifier le wikicode]

Une fonction extérieure est une fonction numérique définie par l'utilisateur (notez la différence entre « externe » et « extérieure »).

Pour les cas simples, la commande deff permet de définir de nouvelles fonctions pouvant s’exprimer avec les opérateurs déjà prédéfinis (les « fonctions internes » ou « primitives »). On passe deux chaînes de caractères comme paramètre :

  • la première indique le nom de la fonction et les variables utilisées en entrée et en sortie, et
  • la deuxième indique la formule.

Par exemple, la fonction

ƒ(x) = 2·x

peut se définir par

deff("[y] = f(x)", "y = 2*x")

On peut aussi utiliser l'instruction function() sous la forme

-->function [y] = f(x)
-->--> endfunction

la formule, sous la forme y = expression, étant placée à la place des trois points. Avec l'exemple ci-dessus, cela donne :

-->function [y] = f(x)
-->y=2*x
-->endfunction

(voir aussi le chapitre Programmation).

La fonction s'utilise par la suite de la même manière qu'une fonction interne (sur l'exemple ci-dessus) :

-->f(2)
 ans  =

    4.

-->feval(1,f) // calcule f(1)
 ans  =

    2.

Graphe d'une fonction[modifier | modifier le wikicode]

Il est souvent utile de pouvoir afficher des fonctions, ainsi voici un exemple de la syntaxe Scilab pour afficher la représentation graphique y = ƒ(x) d'une fonction ƒ.

t=[1:0.1:30]; // t est un vecteur des valeurs de 1 à 30 par pas de 0,1

deff("[y] = f(x)", "y = cos(x)./x")  // définition de la fonction f

plot(t, f(t)) // tracé de la fonction


Pour plus de détails voir : Découvrir Scilab/Graphiques et sons#Tracé de fonction.

Statistiques[modifier | modifier le wikicode]

Pour l'analyse statistique de données, on peut se reporter au document

(en) Anna Bassi et Giovanni Borzi, « Data Analysis and Statistics: a Scilab data mining tutorial » sur Openeering.

Tendance centrale[modifier | modifier le wikicode]

Considérons une matrice X. On peut calculer les moyennes suivantes :

  • mean(X) : calcule la moyenne de tous les éléments de X ;
  • mean(X, "r") : donne un vecteur ligne (row) contenant les moyennes colonne par colonne ;
  • mean(X, "c") : donne un vecteur colonne contenant les moyennes ligne par ligne.

On peut calculer la moyenne géométrique avec geomean(), et la moyenne harmonique avec harmean(), avec la même syntaxe.

De la même manière, on peut calculer les médianes avec median(X), median(X, "r") et median(X, "c").

La commande trimmean(X, p) permet de calculer la moyenne en éliminant les p pourcents extrêmes (les p/2 % valeurs les plus élevées, et les p/2 % valeurs les plus faibles). On peut également faire le calcul par ligne ou par colonne — trimmean(X, p, "r") ou trimmean(X, p, "c"). notons que l'on a :

  • mean(X) = trimmean(X, 0) ;
  • median(X) = trimmean(X, 100).

Dispersion[modifier | modifier le wikicode]

La commande

strange(X)

calcule l'écart (range) entre la valeur maximale et la valeur minimale de la matrice. On peut calculer les valeurs pour chaque colonne ou chaque ligne avec strange(X, "r") et strange(X, "c"). La commande

mad(X)

calcule la déviation absolue moyenne (mean absolute deviation)

on peut de même la calculer pour chaque colonne ou ligne.

La commande

iqr(X)

calcule l'écart interquartile (interquartile range).

La commande stdev calcule l'écart type (standard deviation) et la commande variance la variance des éléments de la matrice (entière, par ligne ou par colonne). Par défaut, la variance est sans biais : normalisée par m×n - 1 pour une matrice de m×n éléments, par m - 1 par colonne et par n - 1 par ligne. Si l'on calcule la variance par ligne ou par colonne, on peut ajouter l'option 1 pour calculer le moment d'ordre 2 (normalisation par m ou n), par exemple

mom2 = variance(X, "r", 1) // pour une matrice non-vecteur
mom2 = variance(X, "*", 1) // pour un vecteur

L'écart type est toujours sans biais. Si l'on veut l'estimateur normalisé par m ou n, il faut utiliser

estim_ecart_type = sqrt(variance(X, "*", 1)) // pour une matrice non-vecteur

Normalisation[modifier | modifier le wikicode]

La commande center(X) permet de centrer la matrice (équivalent de X - mean(X)). La commande wcenter(X) donne la matrice centrée réduite (équivalent de 1/stdev(X)*(X - mean(X))).

Fréquences et quartiles[modifier | modifier le wikicode]

La commande tabul(X) calcule les fréquences d'apparition des valeurs dans la matrice X. Le résultat est classé par défaut par ordre décroissant des valeurs ; pour avoir un classement croissant, on ajoute le paramètre "i" (increasing order) :

tabul(X, "i")

La commande quart(X) calcule les quartiles et les renvoie dans un vecteur de dimension 3. On peut calculer les quartiles pour chaque colonne ou chaque ligne, avec quart(X, "r") et quart(X, "c").

On peut aussi travailler avec des valeurs et leur fréquence. Si la matrice X contient les valeurs et la matrice F leur fréquence — F(i, j) est la fréquence de X(i, j) —, alors

  • meanf(X, F) donne la moyenne ;
  • stdevf(X, F) donne l'écart type ;
  • variancef(X, F) donne la variance.

Notons que F n'est pas nécessairement normalisé à 1 ou à 100 ; il peut par exemple s'agir d'un nombre d'occurrences.

Exemple

X = int(10*rand(100,1, "normal")); // génère 100 valeurs entières avec une loi normale centrée d'écart type 10
F = tabul(X); // calcule les fréquences d'apparition
x = F(:, 1); // extraction de la liste des valeurs
f = F(:, 2); // extraction de la liste des fréquences

quartiles = quart(X); // calcule les quartiles de la population

moyenne = meanf(x, f); // calcul de la moyenne pondérée par les fréquences

Considérons maintenant deux vecteurs X et Y de même dimension, et la matrice F des fréquences de X et Y : F(i, j) est la fréquence de X(i) & Y(j). Alors :

  • covar(X, Y, F) donne la covariance ;
  • correl(X, Y, F) donne le coefficient de corrélation.

Exemple

n = 10; // nombre de valeurs
X = 1:n; // abscisses
Y = 2*X + 5 + 3*rand(X); // ordonnées : droite "bruitée"
plot(X, Y); // tracé des données
F = eye(n, n); // couplage (x, y) : chaque couple est présent une fois, (X(i) & Y(i)) = 1, (X(i) & Y(j)) = 0

covariance1 = covar(X, Y, F);
covariance2 = 1/n*sum(X.*Y) - 1/n^2*sum(X)*sum(Y); // calcul "à la main" pour comparer

correlation1 = correl(X, Y, F);
correlation2 = covariance2/sqrt(variance(X, "*", 1))/sqrt(variance(Y, "*", 1)); // calcul "à la main"

Fonctions de distribution et de répartition de lois statistiques[modifier | modifier le wikicode]

Par défaut, Scilab propose une loi de distribution, pour la loi binomiale, et onze lois de répartition (probabilités cumulatives).

Pour la loi binomiale, la fonction de répartition pour n épreuves (nombre entier) avec une probabilité p (nombre réel entre 0 et 1) est donnée par

f = binomial(p, n)

La variable f est un vecteur de n + 1 composantes, avec

La fonction de répartition s'obtient avec cdfbin() :

P =  = cdfbin("PQ", k, n, p, 1-p)

donne la probabilité P d'avoir entre 0 et k succès pour n tirages avec une probabilité de réussite p. On peut d'ailleurs le comparer avec

f = binomial(p, n)
c = sum(f(1:k+1))

De manière générale, les fonctions de répartition commencent par cdf… (cumulative distribution function), par exemple :

  • P = cdfchi("PQ", x, k) : loi du χ2 à k degrés de liberté ;
  • P = cdfnor("PQ", x, mu, sigma) : loi normale d'espérance mu et d'écart type sigma ;
  • P = cdfpoi("PQ", k, lambda) : loi de Poisson de paramètre lambda ;
  • P = cdft("PQ", t, k) : loi t de Student à k degrés de liberté ;
  • … voir Scilab help >> Statistics > Cumulated Distribution Functions

Notons que pour chaque cas, on peut aussi avoir la valeur complémentaire Q = 1 - P avec la syntaxe

[P, Q] = cdf…("PQ", …)

Exemple

Ci-dessous, nous traçons la fonction de répartition de la loi normale centrée réduite (μ = 0, σ = 1). Notons que pour des calculs vectoriels, l'espérance et l'écart type sont des vecteurs de même dimension que X.

X = -3:0.1:3;
mu = zeros(X);
sigma = ones(X);
[P, Q] = cdfnor("PQ", X, mu, sigma);
plot(X', [P', Q'])

Les mêmes fonctions permettent de rechercher les quantiles ou fractiles :

  • k = cdfbin("S", n, P, 1-P, p, 1-p) : nombre de succès k sur n épreuves ayant pour paramètre p, donnant une probabilité cumulée P ;
  • x = cdfchi("X", k, P, 1-P) : valeur de x pour laquelle on a une probabilité cumulée P avec une loi du χ2 à k degrés de liberté ;
  • x = cdfnor("X", mu, sigma, P, 1-P) : valeur de x pour laquelle on a une probabilité cumulée P avec une loi normale d'espérance mu et d'écart type sigma ;
  • k = cdfpoi("S", lambda, P, 1-P) : valeur de k pour laquelle on a une probabilité cumulée P avec une loi de Poisson de paramètre lambda ;
  • t = cdft("T", k, P, 1-P) : valeur de t pour laquelle on a une probabilité cumulée P avec une loi t de Student à k degrés de liberté ;

Exemple

Ci dessous, nous déterminons le quantile pour un risque α bilatéral de 5 %, c'est-à-dire que l'on rejette 0,25 % des cas les plus hauts et les plus bas, pour une loi normale centrée réduite puis pour une loi t de Student à 30 degrés de liberté

-->x = cdfnor("X", 0, 1, 0.975, 0.025)
 x  =
 
    1.959964
 
-->t = cdft("T", 30, 0.975, 0.025)
 t  =
 
    2.0422725

Enfin, ces fonctions permettent de déterminer les paramètres des loi pour atteindre un quantile visé :

  • loi binomiale :
    • n = cdfbin("Xn", p, 1-p, P, 1-P, k) : nombre d'épreuves ayant pour paramètre p afin d'avoir un quantile P valant k,
    • [p, q] = cdfbin("PrOmpr", P, 1-P, k, n) : paramètre p de chaque épreuve afin d'avoir un quantile P pour k succès parmi n ;
  • loi du χ2 :
    k = cdfchi("Df", P, 1-P, x) : valeur du nombre de degrés de liberté (degrees of freedom) k pour avoir un quantile P valant x ;
  • loi normale :
    • mu = cdfnor("Mean", sigma, P, 1-P, x) : valeur de l'espérance (moyenne) pour avoir un quantile P valant x si l'on a un écart type sigma,
    • sigma = cdfnor("Std", P, 1-P, X, mu) : valeur de l'écart type (standard deviation) pour avoir un quantile P valant x si l'on a une espérance de mu ;
  • loi de Poisson :
    lambda = cdfpoi("Xlam", P, 1-P, k) : valeur du paramètre λ pour avoir un quantile P valant k ;
  • loi t de Student :
    k = cdft("Df", P, 1-P, t) : nombre de degrés de liberté k pour avoir un quantile P valant t.

Générateurs pseudo-aléatoires[modifier | modifier le wikicode]

Les générateurs de nombres aléatoires permettent de tester des programmes ou bien de faire des simulations par la méthode de Monte Carlo.

Scilab possède une fonction rand() qui par défaut génère un nombre entre 0 et 1 avec une loi de distribution uniforme. Par exemple, pour simuler un jet de dé à six faces, on peut utiliser :

de = 1 + int(6*rand())

La fonction peut générer une matrice m×n avec les syntaxes

rand(m, n)
rand(X)

où X est une matrice m×n. Notons que si l'on veut générer n valeurs, il faut bien taper rand(n, 1) ou rand(1, n).

On peut lui faire générer une valeur selon une loi normale centrée réduite avec la clef "normal" :

rand(m, n, "normal")

Exemple

Supposons par exemple que l'on a une poutre faite dans un acier dont la résistance intrinsèque, la limite d'élasticité Re, est garantie au minimum à 235 MPa. En réalité, Re suit une loi normale d'espérance μRe = 269 MPa et d'écart type σRe = 17 MPa. Dans les conditions d'utilisation, l'acier de cette poutre doit résister à une contrainte S (stress) de 214 MPa ; en fait, les conditions de chargement sont variables et suivent une loi normale de paramètres μS = 214 MPa et σS = 15 MPa.

On veut déterminer statistiquement dans combien de cas la contrainte S est supérieure à la résistance Re.

n = 1e6; // nombre d'événements
R = 269 + 17*rand(n, 1, "normal"); // résistance
S = 214 + 14*rand(n, 1, "normal"); // contrainte
Mbool = S>R; // comparaison des valeurs (matrice de booléens)
moyenne = sum(Mbool)/n; // les "vrai" comptent pour 1 et les "faux" pour 0
disp(string(100*moyenne)+" %") // affichage du résultat

La commande grand() (get random number) permet d'avoir des nombres aléatoires selon de nombreuses lois, par exemple :

  • Y = grand(m, n, "bin", N, p) : matrice m×n de nombres aléatoires tirés selon une loi binomiale de N essais avec une probabilité p chacun ;
  • Y = grand(m, n, "nor", mu, sigma) : pour une loi normale d'espérance mu et d'écart type sigma ;
  • Y = grand(m, n, "unf", a, b) : pour une loi uniforme sur [a ; b[ (b exclus) ;
  • Y = grand(m, n, "uin", a, b) : pour une loi uniforme sur [a ; b] (b inclus) ;

Dans l'exemple précédent, on aurait donc :

R = grand(n, 1, "nor", 369, 17); // résistance
S = grand(n, 1, "nor", 314, 14); // contrainte

Autres lois statistiques[modifier | modifier le wikicode]

Des modules complémentaires Atoms permettent d'avoir accès à d'autres loi.

Par exemple le module CASCI donne accès à de nombreuses fonctions de distribution et de répartition, notamment des lois log-normale (pdflognormal(), cdflognormal()), exponentielle (pdfexponential(), cdfexponential()) ou de Weibull (pdfweibull(), cdfweibull()). Ainsi :

p = pdfnormal(x, mu, sigma); // densité d'une loi normale en x
y = cdflognormal(x, lambda); // répartition d'une loi log-normale d'espérance mu et d'écart type sigma
y = cdfexponential(x, mu, sigma); // répartition d'une loi exponentielle de paramètre lambda
y = cdfweibull(x, k, lambda, theta); // répartition d'une loi de Weibull
   // de paramètre de forme k, de paramètre d'échelle lambda et de paramètre de position thêta

Le module fournit également le calcul des quantiles avec les fonctions idf… (inverse cumulative distribution functions), et des générateurs de nombres aléatoires avec les fonctions rnd…. Par exemple :

x = idfweibull(y, k, lambda, theta); // quantile y d'une loi de Weibull
   // de paramètre de forme k, de paramètre d'échelle lambda et de paramètre de position thêta
X = rndweibull(n, k, lambda, theta); // génère n nombres aléatoires selon la loi de Weibull

Le module Stixbox définit également un certan nombre de fonctions statistiques utiles.

Interpolation, extrapolation et lissage[modifier | modifier le wikicode]

Interpolation dans le plan[modifier | modifier le wikicode]

Présentation du problème[modifier | modifier le wikicode]

On dispose de points expérimentaux (xi, yi), qui sont censés suivre une loi y = ƒ(x) inconnue, et l'on désire avoir une valeur y correspondant à une valeur x arbitraire, mais contenue dans l'intervalle [min{xi}, max{xi}]. Une manière de faire consiste à utiliser une interpolation ; Scilab propose deux démarches :

  • interpolation linéaire : si xixxi + 1, alors on prend y sur le segment de droite [(xi, yi) ; ((xi + 1, yi + 1)] :
     ;
  • interpolation par un polynôme de degré trois, en utilisant des splines cubiques d'Hermite :
    y = p(x)
    la fonction p étant un polynôme de degré trois par parties, les portions de polynôme étant raccordées aux points expérimentaux (même valeur, même pente).


Pour plus de détails voir : wikipedia:fr:Interpolation linéaire et wikipedia:fr:Spline cubique d'Hermite.

Interpolation linéaire[modifier | modifier le wikicode]

Le nuage de n points (xi, yi)1 ≤ in est décrit par deux vecteurs : x, qui contient les xi classés par ordre croissant, et y, qui contient les yi correspondant.

Pour un interpolation linéaire, on considère une valeur xp arbitraire, ou un vecteur de valeurs, les valeurs de yp correspondantes sont obtenues par

yp = interp1(x, y, xp)

Plutôt que d'avoir les valeurs interpolées, on peut préférer les valeurs les plus proches ; pour cela, on utilise l'option "nearest" :

yp = interp1(x, y, xp, "nearest")

On peut aussi utiliser la commande d'interpolation à n dimensions linear_interpn() :

yp = linear_interpn(xp, x, y)

Si l'on veut faire une extrapolation, c'est-à-dire avoir une valeur hors de l'intervalle [x1, xn], on peut ajouter une option :

  • "by_zero" : y = 0 hors intervalle ;
  • "C0" : y = y1 si x < x1, y = yn si x > xn ;
  • "natural" : si la fonction est supposée monotone hors de l'intervalle de définition, on utilise le segment de droite le plus proche ;
  • "periodic" : si la fonction est supposée périodique, on répète le motif obtenu.

Par exemple :

yp = linear_interpn(xp, x, y, "natural")

Interpolation par une spline cubique[modifier | modifier le wikicode]

Pour une interpolation par spline cubique d'Hermite, la fonction d'interpolation est décrite par les coordonnées des points et un vecteur de coefficients d (les dérivées), soit trois vecteurs ( x, y, d). Le vecteur d est obtenu par la commande splin() :

d = splin(x, y)

si l'on a maintenant une valeur xp arbitraire, ou un vecteur de valeurs, les valeurs de yp correspondantes sont obtenues par

yp = interp(xp, x, y, d)

Lissage[modifier | modifier le wikicode]

L'interpolation considère que les données sont peu bruitées ; les points mesurés sont « fiables », les valeurs des x et y sont « exactes ». Si le signal est bruité, il faut alors effectuer un lissage.

Lissage par une spline cubique.

Scilab propose la commande lsq_splin() qui utilise des splines cubique comme modèle, déterminées par la méthode des moindres carrés (least square). Si l'on a deux vecteurs x et y de m valeurs chacun, on choisit n abscisses xs pour déterminer la spline, avec dans l'idéal n << m. Typiquement :

// ***** Génération des données *****

m = 200; // nombre de points mesurés
bruit = 0.2; // amplitude du bruit
x = linspace(0, 2*%pi, m);
y = sin(x) + 0.2*rand(x, "normal"); // valeurs bruitées

// ***** Détermination de la spline *****

n = 6; // nombre de points de contrôle de la spline
xs = linspace(x(1), x($), n); // abscisse des points de contrôle
// répartition uniforme
[ys, d] = lsq_splin(x, y, xs); // détermination de la spline

// ***** Calcul des valeurs lissées de y *****

yliss = interp(x, xs, ys, d); // calcul des données lissées

// ***** Tracé *****

plot(x, y) // données originales
plot(x, yliss, "r") // données lissées

Ainsi, Scilab a construit une spline cubique passant par les abscisses xs (en l'occurrence 0, 2π/5, 4π/5, 6π/5, 8π/5 et 2π), et a déterminé les dérivées d en ces n points.

Si l'on veut éviter des effets de bord, on peut choisir une répartition différente des points de contrôle. On utilise fréquemment une répartition sinusoïdale pour avoir plus de points aux bords qu'au centre (comme pour les polynômes de Tchebychev).

theta = linspace(-%pi, 0, n); // angles uniformément répartis
xs = mean([x(1), x($)]) + 0.5*(x($) - x(1))*cos(theta);
// abscisse des points de contrôle

Ceci est en particulier important si l'on veut extrapoler, c'est-à-dire calculer des valeurs de y pour des x hors des points d emesure.

Cela ne présente d'intérêt que pour n ≥ 4.

Lissage et détermination des dérivée et dérivée seconde par l'algorithme de Savitzky-Golay.

L'utilisateur peut évidemment coder un autre algorithme. L'image ci-contre montre le calcul de la dérivée et de la dérivée seconde par l'algorithme de Savitzky-Golay ; cliquer sur l'image pour avoir accès à la page de description indiquant le code utilisé.

On peut aussi effectuer un simple lissage par la méthode des moyennes glissantes en utilisant la fonction convol() (voir ci-après). Si l'on appelle p le nombre de points de l'intervalle glissant :

p = 10; // largeur de l'intervalle

h = 1/p*ones(1, p); // fonction porte pour le lissage

yliss = conv(h, y);

plot(x, y) // données originales
plot(x, yliss(1+p/2:$-p/2+1), "r") // données lissées

Interpolation d'une surface[modifier | modifier le wikicode]

Présentation du problème[modifier | modifier le wikicode]

On dispose de valeurs zi à des positions déterminées (xi, yi), qui sont censés suivre une loi z = ƒ(x, y) inconnue, ce qui définit donc une surface dans l'espace, et l'on désire avoir une valeur z correspondant à un point (x, y) arbitraire.

Interpolation linéaire sur un quadrillage[modifier | modifier le wikicode]

Les données sont décrites par :

  • un vecteur x de dimension m et un vecteur y de dimension n, strictement croissants, qui définissent un quadrillage du plan, donc un ensemble de m×n points {(x1, y1) ; (x1, y2) ; … ; (x1, yn) ; (x2, y1) ; … ; (xm, yn) ;
  • une matrice z de dimension m×n, contenant les valeurs correspondante : zi, j = ƒ(xi, yj).

L'interpolation linéaire se fait de la même manière que précédemment : on définit un nouveau couple de valeurs xp et yp (ou un couple de vecteurs défnissant un nouveau quadrillage), et l'on a :

zp = linear_interp(xp, yp, x, y, z)

c'est-à-dire qu'à partir des points (x, y, z), on définit les points interpolés (xp, yp, zp).

Interpolation cubique sur un quadrillage[modifier | modifier le wikicode]

Comme précédemment, les données sont décrites par :

  • un vecteur x de dimension m et un vecteur y de dimension n, strictement croissants, qui définissent un quadrillage du plan, donc un ensemble de m×n points {(x1, y1) ; (x1, y2) ; … ; (x1, yn) ; (x2, y1) ; … ; (xm, yn) ;
  • une matrice z de dimension m×n, contenant les valeurs correspondante : zi, j = ƒ(xi, yj).

Et comme dans le plan, la fonction d'interpolation est décrite par (x, y, C), les valeurs coordonnées des points du quadrillage et une matrice de coefficients C obtenue par

C = splin2d(x, y, z)

Si l'on considère un point du plan (xp, yp) (ou un couple de vecteurs défnissant un nouveau quadrillage), les valeurs de z sur ces nouveaux points sont :

zp = interp2d(xp, yp, x, y, C)

Interpolation cubique sur un nuage de points dispersés[modifier | modifier le wikicode]

On ne dispose pas toujours de relevés de valeurs sur des coordonnées régulières formant un quadrillage. Dans ce cas-là, on utilise la méthode de Shepard (pondération des plus proches voisins). Les n points expérimentaux (xi, yi, zi) sont décrits par une matrice M de trois colonnes et de n lignes

La fonction d'interpolation est définie par une liste typée T obtenue par

T = cshep2d(M)

Si l'on considère un point du plan (xp, yp) (ou un couple de vecteurs défnissant un nouveau quadrillage), les valeurs de z sur ces nouveaux points sont :

zp = eval_cshep2d(xp, yp, T)

Régression[modifier | modifier le wikicode]

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 :

  • 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 = [X1 ; X2 ; … ; Xn] ;
  • Y = [Y1 ; Y2 ; … ; Yn].

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))2/A(3)) : une fonction gaussienne ;
  • ƒ(A, x) = A(1)·x3 + A(2)·x2 + 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 = ∑i r2(A, Xi, Yi) ;

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[modifier | modifier le wikicode]

Régression linéaire simple[modifier | modifier le wikicode]

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

ce qui peut s'écrire sous forme matricielle

A×X = Y

avec

A = (a b)
Y = (y1 y2yn )

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ésidus.

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 point 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[modifier | modifier le wikicode]

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

ƒ(a1, a2, …, am, b, x1, x2, …xn) = a1·x1 + a2·x2 + … + an·xn + b

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

X = [X11, X12, … X1n ;
     X21, X22, … X2n ;
     …
     Xm1, Xm2, … Xmn]

où Xij 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, … , Yn]

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

Régression polynomiale[modifier | modifier le wikicode]

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 x3 - 2·x2 + 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[modifier | modifier le wikicode]

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

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.

Le problème de la précision et de l'exactitude[modifier | modifier le wikicode]

Scilab est destiné à faire du calcul. La précision, l'exactitude des résultats est donc une préoccupation primordiale.

La première chose est bien évidemment d'utiliser des méthodes, des algorithmes corrects, validés, documentés. Mais la mise en code de ces algorithmes doit également être correcte. Bien évidemment, la méthode choisie doit être retranscrite fidèlement et sans erreur mais il faut aussi prendre en compte certains « effets de bord » dus à la représentation des nombres.

En effet, Scilab calcule par défaut avec des nombre décimaux codés en virgule flottante à double précision selon la norme IEEE754. Cela permet de représenter des nombres dont la valeur absolue est comprise entre 10-307 et 10308 ; si un calcul donne une valeur inférieure, cela retourne la valeur 0 (erreur de soupassement), et s'il donne une valeur supérieure, cela retourne la valeur Inf ou –Inf (erreur de dépassement). Le nombre de chiffres significatifs est de l'ordre de 16 ; ainsi, si deux nombres diffèrent au 17e chiffre significatif, ils seront considérés comme égaux (erreur d'arrondi), leur différence sera nulle, leur rapport sera de 1.

Il faut donc s'assurer d'une part que le valeur finale est représentable mais aussi que les résultats des calculs intermédiaires sont eux aussi représentables. Par exemple le calcul de abs(-10^-200) donne bien le résultat 10-200, alors que sqrt((-10^-200)^2) donne 0, puisque le résultat intermédiaire (-10200)2 n'est pas représentable (soupassement). Dans certains cas, un algorithme itératif, cherchant une valeur approchée, peut donner un meilleur résultat qu'un calcul direct utilisant une formule exacte… Par exemple, pour calculer l’hypoténuse d'un triangle rectangle (norme euclidienne), on peut utiliser la méthode de Moler et Morrison (1983)[1], ce que fait d'ailleurs la commande Scilab norm().

Le problème de la précision peut parfois se manifester sur des calculs « banals ». Par exemple, le vecteur [-5.1:0.2:5.1] devrait contenir les 51 éléments {–5,1 ; –4,9 ; … ; 4,9 ; 5,1}. Or, les cinquante additions successives de la quantité 0,2, du fait des erreurs d'arrondis, font que le dernier élément calculé est légèrement supérieur à 5,1 (il vaut environ 5,1 + 2⋅10-15) et donc la liste s'arrête à 4,9. Ce genre de problème peut se gérer avec la commande nearfloat() :

-5.1:0.2:nearfloat("succ", 5.1)

ou bien, pour la génération de vecteurs, avec la commande linspace().

Si l'on doit faire des calculs avec des cas « extrêmes » — nombres très proches, ou bien ayant une valeur absolue élevée ou faible —, on s'intéressera à des algorithmes « robustes ». Ceux-ci font fréquemment intervenir une normalisation des données ou une factorisation pertinente. Dans la mesure du possible, si l'on a le choix entre une fonction « fait maison » et une fonction Scilab faisant le même travail, on favorisera la fonction Scilab qui est censée être optimisée.

Si des calculs intermédiaires font intervenir des nombres positifs très grands ou très petits, on peut utiliser la fonction logarithme qui a initialement été inventée pour cela. Par exemple, pour calculer

avec

on peut utiliser la fonction logarithme de la fonction gamma, gammaln()[2] :

1-exp(gammaln(365 + 1)-(gammaln(365 - 23 + 1) + 23*log(365)))

On peut également retravailler la formule en

et utiliser

1 - prod( (343/365):(1/365):(364/365) )

Pour effectuer des tests, Scilab dispose de la constante spéciale %eps qui est égale à la plus petite valeur absolue représentable.

Quoi qu'il en soit, la mise en œuvre d'un algorithme doit toujours s'accompagner de tests avec des valeurs pour lesquelles le résultat est connu. Les résultats de référence peuvent être calculés à la main, obtenus avec un logiciel de référence, ou bien les données peuvent être générées en fonction du résultat attendu.

Par ailleurs, un résultat de qualité fait figurer la précision du calcul.

On pourra se reporter à Baudin 2010a.

Voir aussi[modifier | modifier le wikicode]

Dans Wikipédia
Liens externes
Sur Scilab.org
Autres

Notes[modifier | modifier le wikicode]

  1. Baudin 2010b, p. 44–45
  2. voir (en) Computing with very large numbers in Scilab sur liste de discussion Scilab-users et gammaln sur Aide de Scilab

Calculs élémentaires < > Calcul différentiel et intégral