Ajustements d'équations non linéaires avec des exemples MATLAB

Bien qu'il y ait de nombreuses ressources sur la théorie derrière les ajustements (fits) avec équations non linéaires, les exemples, souvent critiques pour la bonne compréhension du modèle, font très souvent défaut. C'est la majeure motivation derrière la rédaction de cet article : offrir des exemples de codes illustrant la théorie.

Les méthodes de Gauss-Newton et de Levenberg-Marquardt sont présentées. La qualité des ajustements est également abordée par le calcul de quantités statistiques (RMSE, déviation standard, erreurs sur les paramètres).

Votre avis et vos suggestions sur cet article m'intéressent !

Alors après votre lecture, n'hésitez pas : 2 commentaires Donner une note à l'article (5)

Article lu   fois.

L'auteur

Site personnel

Liens sociaux

Viadeo Twitter Facebook Share on Google+   

I. Introduction

Ce tutoriel va vous aider à comprendre comment ajuster une équation non linéaire à des données. Je ne suis ni mathématicien ni informaticien. Donc, excusez-moi par avance des inexactitudes que je vais écrire, ou des oublis. Par exemple, je ne vais pas définir mes variables dans un certain groupe donné (pardon, les mathématiciens) et mes exemples codés ne seront pas les plus rapides et certainement pas les plus malins (pardon, les programmeurs). Mais, ça marche ! Et c'est le principal. Ce sujet a été traité maintes fois dans plusieurs livres comme cités dans la bibliographie et sur Internet comme sur la page Méthode des moindres carrés de Wikipédia. Je ne prétends pas réinventer la roue, mais juste donner quelques exemples implémentant les sujets discutés ici.

Tout d'abord, je vais supposer que vous savez comment résoudre un système d'équations linéaires. Ceci est très bien expliqué par Jean-Marc Blanc dans son tutoriel intitulé Résolution des systèmes linéaires.

De manière générale, un système d'équations linéaires a un ensemble de solutions exactes, qui peut être trouvé en une étape, par une inversion de matrice ou par l'utilisation de l'opérateur \ dans MATLAB. Afin de rester général, je vais utiliser l'inversion de matrice et la fonction inv. Par exemple, pour trouver l'ensemble des solutions à une droite de la forme :

formula (1)

C'est-à-dire, trouver les solutions Image non disponible et formula connaissant une série de formula mesurés à une « position » formula (où le mot « position » doit être utilisé de manière générale, car pouvant correspondre à une vitesse, un temps, une température…). De manière générale, l'utilisation de la forme matricielle est plus compacte et va donc être utilisée par la suite. Dans ce cas, l'équation (1) devient :

formula (2)

Dans cette équation, formula correspond aux solutions que l'on veut trouver. Donc ici, formula est un vecteur contenant deux éléments : formula et formula. formula est une matrice dont le nombre de colonnes est égal au nombre de variables/solutions (i.e. deux dans notre cas) et dont le nombre de lignes est égal au nombre de points sur la droite (i.e. la taille de formula). Immédiatement, on peut voir que la première colonne va multiplier le premier élément de formula et la deuxième colonne va multiplier le deuxième élément de formula. Dans ce cas, si l'on considère que le premier élément de formula est formula et le deuxième est formula, alors la première colonne de formula sera le vecteur unité et la deuxième colonne sera formula. Le but est de trouverformula :

formula (3)

Cependant, l'inversion d'une matrice ne peut se faire que si cette matrice est carrée. Ceci est contourné par l'introduction de la transposée de formula qui va être multipliée puis divisée :

formula (4)

On peut écrire dans MATLAB :

 
TéléchargerSélectionnez
function Linear_test

x = [82 ; 271 ; 340.4 ; 431 ; 443.3 ; 516.5 ; 533.6 ; 539.2];
y = [358.5 ; 453.3 ; 488.4 ; 534.1 ; 540.1 ; 576.4 ; 585.2 ; 588.2];

totsize = length(x);
J(1:totsize,1) = ones(totsize,1);
J(1:totsize,2) = x;

solution = inv(J.'*J)*(J.'*y);

D'autres méthodes peuvent être utilisées dans MATLAB pour obtenir le même résultat :

 
Sélectionnez
solution = J\y;

ou :

 
Sélectionnez
solution = lscov(J,y);

ou encore :

 
Sélectionnez
solution = polyfit(x,y,1);

Si on avait voulu un ajustement d'ordre 2, il aurait simplement fallu ajouter une troisième colonne à formula (Image non disponible dont les éléments auraient été formula et la solution formula aurait été un vecteur de trois éléments.

II. Ajustement non linéaire par la méthode des moindres carrés

II-A. Méthode de Gauss-Newton

II-A-1. Explication simpliste

La méthode de Gauss-Newton consiste à transformer un système d'équations non linéaires en un système d'équations linéaires. Prenons une équation non linéaire formula. Elle dépend de deux paramètres formula et formula et est évaluée à divers moments formula (où formula peut être vu comme une « position » au sens général) : formula. Cette fonction devra être minimisée pour n'importe quel formula; c'est-à-dire que formula est implicitement définie comme étant :

formula (5)

formula est la fonction s'ajustant aux données expérimentales formula. Nous supposons que cette équation a un ensemble de solutions, de telle manière que l'on ait le couple de solutions formula pour lequel :

formula (6)

Puis nous développons la fonction en une série de Taylor :

formula (7)

formula et formula correspondent aux dérivées partielles formula et formula respectivement.

L'équation précédente devient :

formula (8)

Avec Image non disponible et formula.

De manière générale, formula et formula représentent les paramètres à ajuster de la fonction. formula et formula sont des valeurs particulières de formula et formula pour lesquelles la fonction formula s'annule.

À ce point, il est important de remarquer plusieurs choses.

Premièrement, la série de Taylor a été tronquée, ce qui signifie que la solution ne peut pas être exacte, mais seulement approximée. À ce titre, la solution devra être trouvée par itérations.

Deuxièmement, nous pouvons alors réécrire formula et formula, où formula et formula sont substitués par formula et formula respectivement, formula désignant la formula itération. formula et formula (formula et formula dans l'équation) sont déjà connus et utilisés. Les seules inconnues sont formula et formula, qui convergent vers les vraies solutions formula et formula après plusieurs itérations, en incrémentant à chaque fois formula et formula : formulaet formula.

Finalement, l'équation est linéaire et peut donc être résolue très simplement pour formula et formulacomme expliqué précédemment.

Il est maintenant nécessaire de généraliser ce que nous avons vu et d'utiliser une formulation matricielle pour le problème. La fonction formula est une fonction de deux variables à optimiser évaluée en n points et dépendant du temps formula. Je rappelle ici que les variables sont les coefficients d'une équation et non les points sur l'axe des ordonnées (ici le temps) qui sont connus. On peut alors définir le vecteur :

formula (9)

Le jacobien du système est alors défini par :

formula (10)

Le vecteur formula permettant l'incrémentation et le vecteur contenant les paramètres :

formula (11)

Donc on a :

formula (12)
   

Il suffit ensuite de répéter ces calculs (calcul de formula et de ses dérivées après chaque incrémentation de formula).

Enfin, il faut que l'utilisateur donne les paramètres initiaux que le programme va utiliser comme point de départ pour la minimisation.

II-A-2. Exemple

Un exemple de problème d'ajustement non linéaire est donné sur cette page : Algorithme de Gauss-Newton. Le code MATLAB permettant de réaliser l'ajustement est donné ci-dessous.

Les données expérimentales formula sont modelées par la fonction :

formula (13)

Les données expérimentales sont connues, de même que formula, et il nous faut déterminer les deux paramètres formula et formula pour ajuster au mieux la fonction aux données. La figure suivante montre la fonction à ajuster pour différentes valeurs de ses paramètres.

Image non disponible
Représentation de la fonction à ajuster pour différentes valeurs de ses paramètres.
 
TéléchargerSélectionnez
function NewtonMethod

%p est un vecteur de dimension 2x1 contenant les parametres

p(1,1) = 0.8;
p(2,1) = 0.4;

%F = y - (p(1).*t)./(p(2)+t); -> forme de l'equation a minimiser

%dF1 = -t./(p(2)+t);          -> forme de la premiere derivee en p(1)
%dF2 = p(1)*x./((p(2)+t).^2); -> forme de la seconde derivee en p(2)

%donnees experimentales
t = [0.038 ; 0.194 ; 0.425 ; 0.626 ; 1.253 ; 2.500 ; 3.740];
D = [0.050 ; 0.127 ; 0.094 ; 0.2122 ; 0.2729 ; 0.2665 ; 0.3317];

%il y aura 5 iterations
for inc = 1:5
    
    %calcul du jacobien
    J(1:length(t),1) = -t./(p(2)+t);
    J(1:length(t),2) = p(1)*t./((p(2)+t).^2);
    
    %calcul des residus
    F = D-(p(1).*t)./(p(2)+t);
    
    %il y a plusieurs manieres de calculer
    %l'incrementation Delta
    Delta = -inv(J.'*J)*(J.'*F);
    %celle-ci est la plus simple sous MATLAB :
    %Delta = -J\F;
    
    %incrementation des parametres
    p = p+Delta;
    
end

%calcul de la RMSE
rmse = (sum(F.^2)/length(F))^0.5;

%calcul de l'ecart type (standard deviation)
stdv = (sum(F.^2)/(length(F)-length(p)))^0.5;

%evalutation des donnees ajustees
ta = 0:0.01:4;
Da = (p(1).*ta)./(p(2)+ta);

%affichage des resultats
fprintf('\nB1 = %f\nB2 = %f\n\n', p);
fprintf('rmse = %f\n', rmse);
fprintf('stdv = %f\n\n', stdv);

%affichage du graphique
figure

h = plot(t, D, 'or', ta, Da, 'b');

legend(h, {'Donnees experimentales' 'Meilleur ajustement'}, 'location', 'southeast')

title('Ajustement avec la methode de Gauss-Newton');
xlabel('t');
ylabel('Donnees a ajuster');

Les valeurs calculées sont les suivantes :

 
Sélectionnez
B1 = 0.361835
B2 = 0.556256

rmse = 0.033475
stdv = 0.039608

Voici un aperçu du graphique avec les données expérimentales et le meilleur ajustement :

Ajustement Gauss-Newton avec MATLAB
Données expérimentales et meilleur ajustement (méthode de Gauss-Newton)

II-B. Méthode de Levenberg-Marquardt

II-B-1. Explication détaillée

Tout d'abord, nous définissons formula, un vecteur qui contient tous les paramètres à ajuster. Le formula paramètre, après formula itérations est défini par :

formula (14)

Si l'on veut ajuster des données expérimentales formula à une équation formula qui dépend de formula (ordonnées), formula étant la formula donnée, qui sont connues et de paramètres formula, alors on peut définir la fonction formula par une série de Taylor tronquée à la première dérivée :

formula (15)

La quantité à minimiser est le résidu :

formula (16)

formula est un vecteur contenant des résidus, mais obtenu avec les « vieux » paramètres (les seuls auxquels on ait accès à une itération formula). Nous voulons minimiser cette quantité :

formula (17)

Celle-ci représente la somme de toutes les « distances » (qui dépend de l'unité des abscisses) entre les valeurs expérimentales et le modèle. Plus cette valeur est petite et plus le modèle représente au mieux les données. La seule variable dans l'équation (17) est le vecteur formula. On dérive donc l'équation (17) par rapport à chaque élément formula et on essaye de trouver formula qui annule la dérivée. Tout d'abord, on peut noter qu'au premier ordre :

formula (18)

Car les résidus sont supposés représenter le minimum et leurs dérivées s'annulent (bien que ce ne soit pas exactement vrai, surtout lors des premières itérations). Deuxièmement, en notant que Image non disponible ne dépend pas explicitement de Image non disponible, il reste :

formula (19)

Si Image non disponible, alors formula. Il nous reste donc :

formula (20)

Ce qui est équivalent à écrire en termes matriciels :

formula (21)

On voit ici que donner de bons paramètres initiaux est essentiel pour éviter que l'ajustement ne diverge, car la dérivée qui s'annule veut dire que l'on a soit un minimum soit un maximum et il vaut mieux bien sûr que notre algorithme cherche le minimum. Ceci dépend en partie de l'initialisation des paramètres.

L'algorithme de Levenberg-Marquardt consiste en l'introduction d'un terme d'amortissement, formula, dans l'équation (21) :

formula (22)

Ce terme permet d'ajuster la fonction en utilisant les caractéristiques de la méthode « steepest-descent » (lorsque les valeurs des paramètres initiaux formula sont loin des valeurs finales, avec formula large) ou les caractéristiques de la méthode de Gauss-Newton (lorsque formula est proche des valeurs finales, avec formula petit). Une valeur typique de formula pour commencer l'ajustement est formula. Ensuite, la RMSE, pour « root mean square error » (ce terme sera défini avec précision plus loin dans ce tutoriel) doit être calculée plusieurs fois (et donc la fonction elle-même) afin de déterminer la nouvelle valeur de formula. À chaque itération, la RMSE est calculée pour formula et pour formula. Ensuite, à chaque itération, on a les conditions :

  • si Image non disponible conduit à une moins bonne RMSE que le Image non disponible original, alors Image non disponible jusqu'à ce qu'une meilleure RMSE soit calculée et le « vieux » formula n'est pas remplacé par le nouveau ;
  • si Image non disponible donne la meilleure RMSE, alors Image non disponible, Image non disponible est remplacé et l'itération continue ;
  • si Image non disponible donne la meilleure RMSE, alors on garde Image non disponible, formula est remplacé et l'itération continue.

Bien que cette méthode fasse appel à plus de calculs, elle est stable, en ce sens que les paramètres formula, au pire ne convergeront pas, mais dans ce cas, ils ne divergeront pas, car formula n'est pas recalculé pendant que formula augmente.

II-B-2. Exemple

L'exemple précédent est repris, en utilisant la méthode de Levenberg-Marquardt :

 
TéléchargerSélectionnez
function LevenbergMarquardtMethod

%p est un vecteur de dimension 2x1 contenant les parametres

p(1,1) = 0.4;
p(2,1) = 0.2;

%F = y - (p(1).*t)./(p(2)+t); -> forme de l'equation a minimiser

%dF1 = -t./(p(2)+t);          -> forme de la premiere derivee en p(1)
%dF2 = p(1)*x./((p(2)+t).^2); -> forme de la seconde derivee en p(2)

%donnees experimentales
t = [0.038 ; 0.194 ; 0.425 ; 0.626 ; 1.253 ; 2.500 ; 3.740];
D = [0.050 ; 0.127 ; 0.094 ; 0.2122 ; 0.2729 ; 0.2665 ; 0.3317];

%initialisation de lambda
lambda0 = 0.01;
rmsea = 15;
rmse = 10;
inc = 0;

%for inc = 1:10
while abs(rmsea-rmse(end))>1e-5 && inc<500 && lambda0 < 1e7
    
    inc = inc+1;
    
    %calcul du jacobien
    J(1:length(t),1) = -t./(p(2)+t);
    J(1:length(t),2) = p(1)*t./((p(2)+t).^2);
    
    %calcul des residus
    F = D-(p(1).*t)./(p(2)+t);
    
    %calul de la RMSE lors de la premiere iteration
    if inc==1
        rmse = (sum(F.^2)/length(F))^0.5;
    else
        %mise a jour de lambda
        lambda = lambda0;
        %generation de la matrice unite
        tailleJJ = length(J.'*J);
        %calul des delta en fonction de lambda
        Delta1 = -inv((J.'*J)+lambda*eye(tailleJJ))*(J.')*F;
        lambda=lambda0/10;
        Delta2 = -inv((J.'*J)+lambda*eye(tailleJJ))*(J.')*F;
        p1 = p;
        p2 = p;
        %mise a jour des parametres dans un nouveau vecteur
        p1 = p1+Delta1;
        %calcul des nouveaux residus et de la nouvelle
        %RMSE pour lambda=lambda0
        F = D-(p1(1).*t)./(p1(2)+t);
        rmse1 = (sum(F.^2)/length(F))^0.5;
        %calcul des nouveaux residus et de la nouvelle
        %RMSE pour lambda=lambda0/10
        p2 = p2+Delta2;
        F = D-(p2(1).*t)./(p2(2)+t);
        rmse2 = (sum(F.^2)/length(F))^0.5;
        %condition sur la mise a jour des parametres
        %voir texte pour plus de details
        if rmse1>rmse && rmse2>rmse
            lambda0 = lambda0*10;
        elseif rmse2<rmse1
            lambda0 = lambda0/10;
            p = p2;
            rmsea = rmse;
            rmse = rmse2;
        elseif rmse2>=rmse1
            p = p1;
            rmsea = rmse;
            rmse = rmse1;
        end
        
    end
    
end

%calcul de la RMSE
rmse = (sum(F.^2)/length(F))^0.5;

%calcul de l'ecart type (standard deviation)
stdv = (sum(F.^2)/(length(F)-length(p)))^0.5;

%evalutation des donnees ajustees
ta = 0:0.01:4;
Da = (p(1).*ta)./(p(2)+ta);

%affichage des resultats
fprintf('\np1 = %f\ntp2 = %f\n\n',p);
fprintf('rmse = %f\n',rmse);
fprintf('stdv = %f\n\n',stdv);

%affichage du graphique
figure

h = plot(t, D, 'or', ta, Da, 'b');

legend(h, {'Donnees experimentales' 'Meilleur ajustement'}, 'location', 'southeast')

title('Ajustement avec la methode de Levenberg-Marquardt');
xlabel('t');
ylabel('Donnees a ajuster');

Les valeurs calculées sont les suivantes :

 
Sélectionnez
p1 = 0.361697
tp2 = 0.555446

rmse = 0.033475
stdv = 0.039608

Voici un aperçu du graphique avec les données expérimentales et le meilleur ajustement :

Ajustement Levenberg-Marquardt avec MATLAB
Données expérimentales et meilleur ajustement (méthode de Levenberg-Marquardt)

III. Calcul de quantités statistiques

III-A. Qualité d'un ajustement

La qualité d'un ajustement peut être évaluée à l'aide de deux critères :

  • la racine de la somme des carrés des résidus ou « root mean square error » (RMSE) en anglais ;
  • l'écart type ou déviation standard ou encore « standard deviation » en anglais.

III-A-1. Racine de la somme des carrés des résidus ou « root mean square error » (RMSE)

La RMSE est la quantité la plus basique, car elle ne dépend que de la distance verticale séparant la donnée expérimentale de la valeur calculée :

formula (23)

Où N est le nombre total de points expérimentaux, formula est une donnée expérimentale et formula est la valeur correspondante calculée. C'est cette valeur qui est minimisée lors du fit.

Si l'on prend l'exemple donné au chapitre II.A, le vecteur formula (résidu) prend les valeurs [0.0269 0.0334 -0.0627 0.0206 0.0223 -0.0295] duquel on calcule que RMSE=0.033.

III-A-2. Écart type ou « standard deviation »

L'écart type est différent de la RMSE en ce sens que l'on prend en compte le nombre de variables formula pour modeler les données. L'écart type dépend alors du nombre de degrés de liberté Image non disponible :

formula (24)

Où N, formula et formula réfèrent aux mêmes quantités présentées précédemment.

Si l'on prend l'exemple donné au chapitre II.A, le vecteur formula (résidu) prend les valeurs [0.0269 0.0334 -0.0627 0.0206 0.0223 -0.0295] duquel on calcule que stdv=0.040 avec formula variables.

III-B. Erreurs sur les données et les paramètres

III-B-1. Pondération des données

Dans certains cas, une donnée expérimentale peut être déterminée avec une précision moindre que les autres données. Il est alors naturel de vouloir donner moins d'importance à cette donnée dans l'ajustement. Le poids est lié à formula représentant l'inverse au carré de l'erreur associée à un point formula; et à la variance formula :

formula (25)

formula si formula et 1 si formula, formula et formula représentant l'index des données. Donc formula est une matrice carrée diagonale de taille formula ou formula est le nombre total de données expérimentales. Il y a plusieurs manières de déterminer l'erreur dans une donnée expérimentale, mais pour simplifier ici, on suppose qu'elle est proche de la RMSE et nous prenons la valeur 0.035. Habituellement, lorsque je ne connais pas l'erreur avec certitude, je prends formula=RMSE;

La variance est définie par l'équation :

formula (26)

formula est le nombre de variables utilisées dans le modèle. Si le modèle choisi ajuste les données et que le bruit qui reste est donc gaussien, formula peut être facilement déterminé et on s'attend à une variance très proche de 1.

La version pondérée des équations (21) (Gauss-Newton) et (22) (Levenberg-Marquardt) est :

formula (27)

et :

formula (28)

Dans ce cas, la valeur minimisée devient la RMSE pondérée :

formula (29)

L'écart type pondéré devient :

formula (30)

Il est bon de noter ici que seul le poids relatif des données est important pour la minimisation. Si les poids sont tous les mêmes, alors avoir tous les poids à 5000 ou 1 ne changera rien. Ce qui compte, c'est le poids relatif, c'est-à-dire à quel point les poids sur une ou plusieurs mesures sont différents des autres.

Un exemple d'application de la méthode de Gauss-Newton utilisant la pondération est donné ci-dessous. Ici, c'est la RMSE pondérée qui est minimisée. Pour tester l'algorithme, vous pouvez changer une valeur expérimentale et augmenter fortement l'erreur associée et vérifier que les valeurs finales ne changent que très peu, que la RMSE pondérée reste quasiment inchangée, alors que la RMSE explose.

 
TéléchargerSélectionnez
function NewtonMethod_ponderation

%p est un vecteur de dimension 2x1 contenant les parametres

p(1,1) = 0.8;
p(2,1) = 0.4;


%F = y - (p(1).*t)./(p(2)+t); -> forme de l'equation a minimiser

%dF1 = -t./(p(2)+t);          -> forme de la premiere derivee en p(1)
%dF2 = p(1)*x./((p(2)+t).^2); -> forme de la seconde derivee en p(2)

%donnees experimentales
t = [0.038 ; 0.194 ; 0.425 ; 0.626 ; 1.253 ; 2.500 ; 3.740];
D = [0.050 ; 0.127 ; 0.094 ; 0.2122 ; 0.2729 ; 0.2665 ; 0.3317];

%ponderation
weight = [0.03 0.05 0.03 0.04 0.03 0.03 0.06];
w = eye(length(t));
w(w~=0) = 1./(weight.^2);

%il y aura 5 iterations
for inc = 1:5
    
    %calcul du jacobien
    J(1:length(t),1) = -t./(p(2)+t);
    J(1:length(t),2) = p(1)*t./((p(2)+t).^2);
    
    %calcul des residus
    F = D-(p(1).*t)./(p(2)+t);
    %calcul de la variance
    variance2 = sum(((F.^2).')./(weight.^2))./(length(t)-length(p));
    %calul de la matrice de ponderation
    w = variance2*w;
    %il y a plusieurs manieres de calculer l'incrementation Delta
    Delta = -inv(J.'*w*J)*(J.'*w*F);
    
    %incrementation des parametres
    p = p+Delta;
    %retour à la matrice de ponderation
    w = w/variance2;
    
end

%calcul de la RMSE
rmse = (sum(F.^2)/length(F))^0.5;

%calcul de la RMSE ponderee
rmsew = (sum(diag(w).*(F.^2))/sum(diag(w)))^0.5;

%calcul de l'ecart type (standard deviation)
stdv = (sum(F.^2)/(length(F)-length(p)))^0.5;

%evalutation des donnees ajustees
ta = 0:0.01:4;
Da = (p(1).*ta)./(p(2)+ta);

%affichage des resultats
fprintf('\nB1 = %f\nB2 = %f\n\n',p);
fprintf('rmse = %f\n',rmse);
fprintf('rmsew = %f\n',rmsew);
fprintf('stdv = %f\n\n',stdv);

%affichage du graphique
figure

h = plot(t, D, 'or', ta, Da, 'b');

legend(h, {'Donnees experimentales' 'Meilleur ajustement'}, 'location', 'southeast')

title('Ajustement avec la methode de Gauss-Newton et ponderation');
xlabel('t');
ylabel('Donnees a ajuster');

Les valeurs calculées sont les suivantes :

 
Sélectionnez
B1 = 0.361939
B2 = 0.625569

rmse = 0.034376
rmsew = 0.035185
stdv = 0.040675

Voici un aperçu du graphique avec les données expérimentales et le meilleur ajustement :

Ajustement Gauss-Newton avec pondération des données avec MATLAB
Données expérimentales et meilleur ajustement (méthode de Gauss-Newton avec pondération des données)

III-B-2. Erreurs sur les paramètres - Matrice de corrélation

Avant de rentrer dans le vif du sujet, il est important de prendre conscience que les calculs qui suivent et qui permettent d'évaluer l'erreur sur les paramètres ne sont valides que si et seulement si la fonction est appropriée pour décrire les données. C'est-à-dire que si la fonction n'est pas appropriée, alors les paramètres n'auront pas de sens et encore moins l'erreur sur les paramètres. De plus, les erreurs systématiques, souvent très dures à trouver (car cachées et donc non détectées), ne sont pas prises en compte. En gros, ces erreurs n'ont de sens que si le bruit est supposé gaussien et donc aléatoire. Finalement, l'autre point à noter est que pour que les paramètres aient du sens, il faut généralement avoir assez de points par paramètre. Généralement, une valeur de dix points par paramètre (ou plus) est suffisante pour que l'erreur donnée corresponde à 1formula ou une déviation standard, ou bien il faut que le nombre de degrés de liberté soit large (Image non disponible, en rappelant que le nombre de degrés de liberté est défini par Image non disponible). Je vais continuer à utiliser le même modèle avec les mêmes données, pour assurer une continuité entre les divers chapitres de cet article, mais sachez que je n'ai pas assez de points par paramètre pour m'assurer que le vrai paramètre est contenu dans formula à 68,3 %.

L'équation permettant de calculer la déviation standard pour chaque paramètre est :

formula (31)

Image non disponible est un vecteur de même dimension que le nombre de paramètres, contenant les déviations standards de chaque paramètre, rangées dans le même ordre que dans le jacobien. Image non disponible est définie dans l'équation (26). On retrouve ici la matrice de pondération formula. On peut voir ici toute son importance, car les erreurs sur les paramètres dépendent directement des incertitudes liées aux mesures. Alors que dans les équations (27) et (28), seul le poids relatif de chaque mesure par rapport aux autres est important, le poids absolu devient critique pour bien calculer l'erreur associée à chaque paramètre. Personnellement, j'affiche toujours la valeur formula, qui me dit combien de fois mon paramètre est plus grand que sa déviation standard. Dans le cas où ce nombre est inférieur à 1, cela veut dire que l'erreur sur mon paramètre est plus grande que le paramètre lui-même, une indication que ce paramètre n'est sûrement pas utile ou approprié dans le modèle utilisé, ou encore qu'il n'y a pas assez de données pour le caractériser.

La matrice de covariance, de taille formula est calculée de manière très similaire, en utilisant l'équation :

formula (32)

De laquelle on tire la matrice de corrélation :

formula (33)

On peut voir que toutes ces matrices et valeurs dépendent du jacobien, qui représente les dérivées de la fonction s'ajustant aux données. La matrice de covariance permet de déterminer les relations entre paramètres, c'est-à-dire, de quelle manière un paramètre formula va être affecté si je change la valeur d'un paramètre formula (en le fixant), tout en continuant d'essayer de minimiser mes données. Est-ce que le paramètre formula est très dépendant de la valeur que le paramètre formula va prendre ? Par exemple, est-ce qu'un terme constant introduit pour corriger un éventuel décalage (offset) des données est important pour déterminer les paramètres déterminant la courbe des données ? Intuitivement, la réponse est non. La matrice de corrélation est la matrice de covariance normalisée, pour pouvoir discuter des dépendances entre les divers paramètres. Dans l'idéal, on veut des dépendances minimales entre paramètres, pour que les solutions soient le plus orthogonales possible et éviter d'avoir des paramètres qui vont être en compétition pour représenter une même caractéristique de la courbe à ajuster. Par exemple, essayer de modéliser une gaussienne avec un terme quadratique d'ajustement de la ligne de base est dangereux, car le terme quadratique peut essayer de partiellement s'ajuster à la gaussienne (si celle-ci est bruitée par exemple), ce qui aura une incidence sur la valeur finale des paramètres décrivant la gaussienne. Une forte corrélation entre un ou plusieurs paramètres de la gaussienne et ce terme quadratique est alors attendue.

Les paramètres et les erreurs sont généralement affichés en donnant la valeur du paramètre puis la déviation standard, entre parenthèses, là où s'arrête le dernier chiffre significatif. Voici un exemple provenant directement du code MATLAB donné ci-après. Image non disponible et formula. On peut alors écrire formula ou formula ce qui signifie formula et formula respectivement. Pour les deux autres paramètres (formula, formula, formula et formula), on aurait formula et formula. Certaines personnes donnent les paramètres avec Image non disponible (la confiance dans l'intervalle donné étant alors de 95,5 %).

Les diverses quantités introduites dans ce chapitre sont calculées dans l'exemple de code ci-dessous. Un décalage a été ajouté sur les données expérimentales formula et le modèle a été ajusté pour prendre en compte ce changement. Le modèle est donc maintenant : Image non disponible.

 
TéléchargerSélectionnez
function NewtonMethod_ponderation_error_parameters

%p est un vecteur de dimension 3x1 contenant les parametres

p(1,1) = 0.8;
p(2,1) = 0.4;
p(3,1) = 0.5;

%F = y - (p(1).*t)./(p(2)+t); -> forme de l'equation a minimiser

%dF1 = -t./(p(2)+t);          -> forme de la premiere derivee en p(1)
%dF2 = p(1)*x./((p(2)+t).^2); -> forme de la seconde derivee en p(2)

%donnees experimentales
t = [0.038 ; 0.194 ; 0.425 ; 0.626 ; 1.253 ; 2.500 ; 3.740];
D = [0.050 ; 0.127 ; 0.094 ; 0.2122 ; 0.2729 ; 0.2665 ; 0.3317];
D = D + 0.57;
 
weight = [0.03 0.02 0.035 0.05 0.04 0.035 0.07];
w = eye(length(t));
w(w~=0) = 1./(weight.^2);

%il y aura 10 iterations
for inc = 1:5
    
    %calcul du jacobien
    J(1:length(t),1) = -t./(p(2)+t);
    J(1:length(t),2) = p(1)*t./((p(2)+t).^2);
    J(1:length(t),3) = -ones(1,length(t));
    
    %calcul des residus
    F = D-((p(1).*t)./(p(2)+t)+p(3));
    %calcul de la variance
    variance2 = sum(((F.^2).')./(weight.^2))./(length(t)-length(p));
    %calul de la matrice de ponderation
    w = variance2*w;
    %il y a plusieurs manieres de calculer
    %l'incrementation Delta
    Delta = -inv(J.'*w*J)*(J.'*w*F);
    
    %incrementation des parametres
    p = p+Delta;
    w = w/variance2;
    
end

%Deviation standard des parametres
stdp = sqrt(variance2)*diag(sqrt((inv((J.')*w*J))));

%Parametres divises par stdp
p_stdp = p./stdp;

%Calcul de la matrice variance-covariance pour calculer
%la matrice de correlation qui va etre affichee a l'ecran.
%La matrice est symetrique et seulement une moitie doit etre calculee.
covmat = sqrt(variance2)*inv(J.'*w*J);                
for inc = 1:3
    for tt = 1:inc
        CorrelationMatrix(inc,tt) = covmat(inc,tt)/sqrt(covmat(inc,inc)*covmat(tt,tt));
    end
end

%calcul de la RMSE
rmse = (sum(F.^2)/length(F))^(0.5);

%calcul de la RMSE ponderee
rmsew = (sum(diag(w).*(F.^2))/sum(diag(w)))^0.5;

%calcul de l'ecart type (standard deviation)
stdv = (sum(F.^2)/(length(F)-length(p)))^(0.5);

%evalutation des donnees ajustees
ta = 0:0.01:4;
Da = (p(1).*ta)./(p(2)+ta) + p(3);

%affichage des resultats
fprintf('\np1 = %f\np2 = %f\np3 = %f\n\n',p);
fprintf('error_p1 = %f\nerror_p2 = %f\nerror_p3 = %f\n\n',stdp);
fprintf('p1_stdp1 = %f\np2_stdp2 = %f\np3_stdp3 = %f\n\n',p_stdp);
  
fprintf('          p(1)    p(2)    p(3)');
for inc = 1:3
    fprintf('\np(%d)\t', inc);
    for tt = 1:inc
        fprintf(' %4.3f\t',CorrelationMatrix(inc,tt));
    end
end
fprintf('\n\n');

fprintf('rmse = %f\n',rmse);
fprintf('rmsew = %f\n',rmsew);
fprintf('stdv = %f\n\n',stdv);

figure
h = plot(t, D, 'or', ta, Da, 'b');

legend(h, {'Donnees experimentales' 'Meilleur ajustement'}, 'location', 'southeast')

title({'Ajustement avec la methode de Gauss-Newton' 'avec calcul des erreurs sur les parametres'});
xlabel('t');
ylabel('Donnees a ajuster');

Les valeurs calculées sont les suivantes :

 
Sélectionnez
p1 = 0.311792
p2 = 0.774445
p3 = 0.617076

error_p1 = 0.074810
error_p2 = 0.726572
error_p3 = 0.041097

p1_stdp1 = 4.167804
p2_stdp2 = 1.065890
p3_stdp3 = 15.015118

          p(1)    p(2)    p(3)
p(1)     1.000
p(2)     0.604     1.000
p(3)     0.085     0.793     1.000

rmse = 0.032253
rmsew = 0.029378
stdv = 0.042667

Voici un aperçu du graphique avec les données expérimentales et le meilleur ajustement :

Ajustement Gauss-Newton avec pondération des données et erreurs sur les paramètres avec MATLAB
Données expérimentales et meilleur ajustement (méthode de Gauss-Newton avec calcul d'erreurs sur les paramètres)

IV. Conclusion

J'espère que cet article pourra vous être utile pour commencer à faire des ajustements non linéaires par vous-même, ou au moins vous aider à y voir plus clair sur ce que qui se passe lors d'une minimisation. Les méthodes proposées sont parmi les plus communément utilisées pour faire des ajustements linéaires : méthode de Gauss-Newton et de Levenberg-Marquardt. L'utilisation de la pondération permet de mieux traiter les données en minimisant l'importance de celles qui sont les plus bruitées. Enfin, les calculs de la RMSE (pondérée ou non), de la déviation standard (pondérée ou non) ainsi que l'erreur sur les paramètres, permettent de déterminer la qualité de l'ajustement. Tout ceci mis bout à bout permet d'effectuer une analyse complète de données par un ajustement non linéaire.

Bien sûr, je ne suis pas à l'abri des erreurs et si vous en voyez. N'hésitez pas à me contacter pour me dire ce qui ne va pas !

Vos remarques et corrections sont les bienvenues dans cette discussion 2 commentaires Donner une note à l'article (5)

Il reste encore bien des sujets à couvrir, comme coder des simplexes (qui ne font pas appel à des dérivées), ou la minimisation orthogonale, qui pourront s'ajouter à ce tutoriel, mais plus tard.

Voilà, bon courage dans le monde de la minimisation !

V. Bibliographie

Voici quelques liens qui vous permettront d'approfondir le sujet développé dans ce tutoriel :

VI. Remerciements

Je remercie Dut pour l'aide apportée pendant la rédaction de ce tutoriel ainsi que ClaudeLELOUP pour ses corrections orthographiques.

En complément sur Developpez.com

  

Copyright © 2014 Damien Forthomme. Aucune reproduction, même partielle, ne peut être faite de ce site et de l'ensemble de son contenu : textes, documents, images, etc. sans l'autorisation expresse de l'auteur. Sinon vous encourez selon la loi jusqu'à trois ans de prison et jusqu'à 300 000 € de dommages et intérêts.