Ajustement de courbe

Ajustement par itérations d'une courbe bruitée par un modèle de pic asymétrique (méthode de Gauss-Newton avec facteur d'amortissement variable).

L'ajustement de courbe est une technique d'analyse d'une courbe expérimentale, consistant à construire une courbe à partir de fonctions mathématiques et d'ajuster les paramètres de ces fonctions pour se rapprocher de la courbe mesurée — on parle donc aussi d'ajustement de paramètres. On utilise souvent le terme anglais curve fitting, profile fitting ou simplement fitting, pour désigner cette méthode ; on utilise souvent le franglais « fitter une courbe » pour dire « ajuster une courbe ».

On utilise des méthodes de régression. Dans les cas simples, il s'agit de régression multilinéaire si la loi est linéaire pour tous les paramètres, ou de régression polynomiale lorsque l'on utilise un polynôme pour simuler le phénomène (les paramètres physiques pouvant être déduits des coefficients du polynôme).

Les méthodes de régression classiques permettent de déterminer les paramètres à partir de calculs sur les données, mais sont inapplicables si la fonction est trop complexe. Il faut alors travailler par essai-erreur pour se rapprocher d'une solution, au sens de la méthode des moindres carrés. La solution n'est pas nécessairement unique.

Ajustement algébrique ou géométrique

On distingue deux types d'ajustement :

• les ajustements algébriques consistent à minimiser l'écart vertical (en y) entre la courbe modèle et les points expérimentaux ;
• les ajustements géométriques consistent à minimiser la distance perpendiculairement à la courbe modèle ; c'est le cas par exemple de la régression circulaire ou elliptique, qui consiste à trouver le cercle, resp. l'ellipse, la plus proche des points

Dans le cas d'un ajustement géométrique, on parle de la méthode des moindres carrés totaux (total least square, TLS) : en effet, on prend en compte les deux coordonnées x et y pour la détermination de l'écart quadratique.

Fonction modèle utilisée

Régression sur un nuage de points par un polynôme de degré croissant.

Dans certains cas, on a un modèle théorique permettant de prévoir la forme de la courbe ; la méthode d'ajustement permet de déterminer les paramètres de l'échantillon. Dans d'autres cas, on utilise une fonction empirique ; on s'intéresse alors en général à la surface, la largeur ou à la position du maximum de la fonction.

On utilise souvent des polynômes : ce sont des modèles qui se calculent facilement. Pour les formes de type « pic » (distributions unimodales), on utilise fréquemment des fonctions gaussiennes, lorentziennes ou bien des combinaisons (fonctions ou pseudo-fonctions de Voigt), ou encore des fonctions de Pearson. Les courbes présentant un amortissement comportent fréquemment une composante exponentielle négative (fonction en e^-x) ; les courbes en S peuvent être modélisées par une fonction sigmoïde.

De manière générale, on a une fonction f modèle ayant n paramètres p¹, p², …, pⁿ qui relie l'abscisse x à l'ordonnée y :

$y=f\left(p^{1},p^{2},\ldots ,p^{n},x\right)$

et on compare cette fonction avec les m points expérimentaux

$\left[\left(x_{1},y_{1}^{\text{exp}}\right),\left(x_{1},y_{1}^{\text{exp}}\right),\ldots ,\left(x_{m},y_{m}^{\text{exp}}\right)\right].$

La fonction f peut être parfois décomposée en plusieurs fonctions f₁, f₂… soit qu'elle en est la somme, le produit, le produit de convolution…

$f=f_{1}+f_{2}+\ldots \qquad {\text{ ou }}\qquad f=f_{1}\times f_{2}\times \ldots \qquad {\text{ ou }}\qquad f=f_{1}*f_{2}*\ldots$

Démarche de régression

Considérations générales

On calcule les points y_i^cal de la courbe simulée :

$y_{i}^{\text{cal}}=f\left(x_{i}\right)$.

Habituellement, on utilise un algorithme de régression visant à minimiser l'écart quadratique entre la courbe simulée et la courbe expérimentale ; on parle d'algorithme de minimisation de l'erreur.

On définit donc un facteur de fiabilité (reliability factor) R :

$\mathrm {R} ={\sqrt {\frac {\sum _{i}\left(y_{i}^{\mathrm {exp} }-y_{i}^{\mathrm {cal} }\right)^{2}}{\sum _{i}{y_{i}^{\mathrm {exp} }}^{2}}}}$

où y_i^exp est le i-ème point mesuré et y_i^cal est le i-ème point calculé. Le facteur R est similaire dans son expression au coefficient de corrélation multiple (R² étant alors le coefficient de détermination). On utilise plus couramment le facteur de fiabilité pondérée (weighted reliability factor) R_wp :

$\mathrm {R_{wp}} ={\sqrt {\frac {\sum _{i}w_{i}\cdot \left(y_{i}^{\mathrm {exp} }-y_{i}^{\mathrm {cal} }\right)^{2}}{\sum _{i}w_{i}\cdot {y_{i}^{\mathrm {exp} }}^{2}}}}$

où w_i est le poids attribué au point i ; ce poids représente l'incertitude associée au point i.

La régression consiste à trouver le jeu de paramètres (p¹, p², …, pⁿ) tel que R_wp est minimal. On a donc notamment

$\forall i,{\frac {\partial \mathrm {R_{wp}} }{\partial p^{i}}}=0$.

Dans les cas simples, le système d'équations que l'on obtient peut se résoudre « simplement » ; en particulier pour les régressions multilinéaires, on peut le résoudre par inversion de matrice.

Régression non linéaire

Lorsqu'il n'est pas possible de résoudre simplement le système d'équations, on peut utiliser la méthode de Newton-Raphson.

C'est un algorithme itératif : à chaque étape, on regarde de quelle manière une petite variation d'un paramètre fait varier la fonction de fiabilité R. Le jeu de paramètre retenu pour l'étape suivante est obtenue par extrapolation de ces petites variations, par linéarisation : on cherche les valeurs des paramètres qui rendraient R nul si les variations de R étaient proportionnelles aux variations des paramètres. Le calcul s'arrête lorsque l'on n'arrive plus à diminuer le facteur de fiabilité, on parle de convergence.

Si l'on note (p^k) l'ensemble des paramètres, l'intensité calculée en chaque point i à l'étape j s'exprime par

$y_{i}^{\mathrm {cal} }=f\left(p_{j}^{1},p_{j}^{2},\ldots ,p_{j}^{n},x_{i}\right)$.

Pour simplifier les calculs, on peut faire un développement limité du premier ordre de cette fonction $f$, alors en appliquant des incréments Δp^k aux paramètres, on peut écrire :

$f\left(p_{j}^{1}+\Delta p_{j}^{1},\ldots ,p_{j}^{n}+\Delta p^{n},x_{i}\right)\simeq f\left(p_{j}^{1},\ldots ,p_{j}^{n},x_{i}\right)+\sum _{k}{\frac {\partial f}{\partial p^{k}}}\left(x_{i}\right)\cdot \Delta p^{k}$

en imposant y_i^cal = y_i^exp, on obtient ainsi un système d'équations linéaires (l'indice j est omis pour plus de clarté) :

${\begin{pmatrix}{\frac {\partial f}{\partial p^{1}}}(x_{1})&{\frac {\partial f}{\partial p^{2}}}(x_{1})&\cdots &{\frac {\partial f}{\partial p^{n}}}(x_{1})\\{\frac {\partial f}{\partial p^{1}}}(x_{2})&{\frac {\partial f}{\partial p^{2}}}(x_{2})&\cdots &{\frac {\partial f}{\partial p^{n}}}(x_{2})\\\vdots &\vdots &&\vdots \\{\frac {\partial f}{\partial p^{1}}}(x_{m})&{\frac {\partial f}{\partial p^{2}}}(x_{m})&\cdots &{\frac {\partial f}{\partial p^{n}}}(x_{m})\\\end{pmatrix}}{\begin{pmatrix}\Delta p^{1}\\\Delta p^{2}\\\vdots \\\Delta p^{n}\end{pmatrix}}+{\begin{pmatrix}f\left(p^{1},\ldots ,p^{n},x_{1}\right)\\f\left(p^{1},\ldots ,p^{n},x_{2}\right)\\\vdots \\f\left(p^{1},\ldots ,p^{n},x_{m}\right)\\\end{pmatrix}}={\begin{pmatrix}y_{1}^{\mathrm {exp} }\\y_{2}^{\mathrm {exp} }\\\vdots \\y_{m}^{\mathrm {exp} }\\\end{pmatrix}}$

qui peut se résoudre par une pseudo-inversion de matrice, ce qui permet de calculer les Δp^k et donc les valeurs des paramètres à l'étape suivante j+1 :

si l'on appelle A la matrice $\left({\frac {\partial f}{\partial p^{k}}}(x_{i})\right)_{i,k}$, le système d'équations s'écrit :

$A\left(\Delta p_{j}^{k}\right)_{k}+\left(y_{i}^{\text{cal}}\right)_{i}=\left(y_{i}^{\text{exp}}\right)_{i}$

La matrice pseudo-inverse B a les mêmes dimensions que ^tA (la matrice transposée de A) et vérifie :

• A × B × A = A ;
• B × A × B = B ;
• A × B et B × A sont hermitiennes.

La matrice pseudo-inverse peut se calculer par décomposition en valeurs singulières.

On prend alors

$\left(\Delta p_{j}^{k}\right)_{k}=B\left[\left(y_{i}^{\text{exp}}\right)_{i}-\left(y_{i}^{\text{cal}}\right)_{i}\right]$

et

$\left(p_{j+1}^{k}\right)_{k}=\left(p_{j}^{k}\right)_{k}+\left(\Delta p_{j}^{k}\right)_{k}$.

La pseudo-inversion est justifiée par le fait que le premier membre du système d'équations s'écrit alors :

$A\times \left(\Delta p_{j}^{k}\right)_{k}=A\times B\times \left[\left(y_{i}^{\text{exp}}\right)_{i}-\left(y_{i}^{\text{cal}}\right)_{i}\right]$ ;

si on le multiplie par B, il devient :

$B\times A\times B\times \left[\left(y_{i}^{\text{exp}}\right)_{i}+\left(y_{i}^{\text{cal}}\right)_{i}\right]=B\times \left(y_{i}^{\text{exp}}\right)_{i}+B\times \left(y_{i}^{\text{cal}}\right)_{i}$

et donc on a bien

$B\times \left[A\times \left(\Delta p_{j}^{k}\right)_{k}-\left(y_{i}^{\text{cal}}\right)_{i}\right]=B\times \left(y_{i}^{\text{exp}}\right)_{i}$.

Notons qu'il s'agit d'une condition suffisante, mais pas nécessaire.

Il s'agit bien d'une méthode de Newton-Raphson, puisque l'on linéarise la fonction et que l'on cherche le point qui annule R (R = 0 si y_i^cal = y_i^exp pour tout i).

En remarquant que A est en fait la matrice jacobienne de $f$ (A = J), on voit que le système d'équations est constitué des équations normales. Les incréments des paramètres peuvent alors se calculer par :

$\left(\Delta p^{k}\right)_{k}=-\left({^{\text{t}}}\!J\times J\right)^{-1}\times {^{\text{t}}}\!J\times \left[\left(y_{i}^{\text{exp}}\right)_{i}-\left(y_{i}^{\text{cal}}\right)_{i}\right]$.

On retrouve alors l'algorithme de Gauss-Newton.

Parfois, la linéarisation est « trop violente » : si la fonction présente une « courbure » importante lorsque les paramètres varient, l'application de la « loi tangente » peut éloigner de la solution plutôt que rapprocher. On peut donc utiliser un coefficient d'atténuation d ≤ 1, on a alors :

$\left(p_{j+1}^{k}\right)_{k}=\left(p_{j}^{k}\right)_{k}+d\cdot \left(\Delta p_{j}^{k}\right)_{k}$.

On peut aussi utiliser un algorithme plus complexe utilisant les dérivées secondes.

Il peut arriver que le programme converge vers un minimum local de R qui n'est pas le minimum absolu. Pour éviter cette situation, on utilise l'algorithme de Metropolis-Hastings : lorsque l'on a convergé, on fait varier les paramètres d'une certaine valeur et on recommence le processus pour voir si l'on arrive à converger vers un autre jeu de paramètres ayant un facteur de fiabilité plus faible.

Si la simulation était parfaite, le facteur de fiabilité aurait une valeur dépendant du rapport signal sur bruit. Si l'on sait calculer ce rapport signal sur bruit, c'est-à-dire si l'on connaît la loi de probabilité régissant les fluctuations du signal, on peut alors déterminer un facteur de fiabilité incompressible R₀. Si la simulation est parfaite, on a alors R = R₀ de fait, la qualité de la simulation est souvent exprimée par le rapport R/R₀, qui doit tendre vers 1 au fur et à mesure des étapes itératives.

Régression circulaire

La régression circulaire consiste à trouver le « meilleur cercle », au sens des moindres carrés, décrivant un ensemble de points. C'est un problème de régression géométrique non linéaire. Cependant, un choix astucieux de la fonction d'erreur permet de se ramener à un problème linéaire.

Régression elliptique

Régression elliptique.

La régression elliptique consiste à trouver la « meilleure ellipse », au sens des moindres carrés, décrivant un ensemble de points.

Applications

• Déconvolution
• Méthode de Rietveld
• Compression de courbe

Voir aussi

Articles connexes

• Approximation de fonction
• Interpolation (mathématiques)
• Régression multilinéaire
• Algorithme de Levenberg-Marquardt