Cet article, qui s’adresse avant tout aux plus scientifiques d’entre vous (quoique…) a pour but d’automatiser une tâche d’ajustements (ou fit par anglicisme) de courbes et d’en récupérer les informations (coefficients et compagnie). Je vais en profiter pour faire l’apologie du logiciel que je vais utiliser : gnuplot. La motivation de cet article fait suite à certaines requêtes de personnes visitant mon blog qui doivent être insatisfaites.

Position du problème

Le problème est le suivant. On dispose de deux fichiers. L’un (weight.txt) possède deux colonnes : abscisse et ordonnée d’un nuage de point. Sur ce nuage de point, je désire connaître une pente locale sur un certain intervalle. Justement, ces intervalles sont donnés dans un second fichier (range.txt) sous la forme de deux colonnes : début et fin.

On voit donc se dessiner la petite mécanique : prendre un intervalle dans le fichier range.txt, effectuer un fit sur les données de weight.txt à l’aide d’une droite, et envoyer la pente (et son incertitude au moins) dans un fichier tiers.

Gnuplot

Gnuplot est un logiciel de tracé et de tracé uniquement (mais qui n’a rien à voir avec GNU). Il ne permet pas de traiter des données, mais au moins, ce qu’il fait, il le fait bien. Enfin, quand je dis tracé seulement, nous n’en ferons aucun ici. La petite exception, c’est l’ajustement. Son usage se fait en console ce qui lui vaut les foudre de certains : « c’est d’un autre temps ». Oui, gnuplot est un vieux logiciel, mais si il a traversé les âges, c’est bien une preuve de son efficacité.
Justement, le fait qu’il soit un logiciel de type CLI (command line interface) permet d’automatiser des tâches par des « scripts » si on peut les nommer comme cela. Rapidement, les avantages. Automatisation comme je l’ai dit, adaptation rapide et conservation du « moyen » utilisé pour générer le graphique car je conserve toujours un graphique avec son fichier de commande gnuplot. Même si un jour je venais à changer de logiciel, j’ai un fichier texte interopérable qui me dit comment j’ai fait et que j’ai pu commenter et que je peux réadapter.

Bref, gnuplot c’est ce qu’il nous faut !

(A noter qu’il existe des interfaces graphiques (GUI) pour gnuplot comme cet addon pour emacs ; jamais testé par mes soins, je n’en vois pas l’intérêt ).

Automatisation

Afin d’automatiser complétement le traitement, je vais écrire un petit script perl comme j’aime le faire. Plutôt que de faire un long discours, voilà ce script que je vais commenter.

Après les ouvertures de fichiers, j’effectue les choses suivantes à chaque intervalle :

Écriture d’un fichier gnuplot.gnu avec la ligne

fit [$range[0]:$range[1]] coeffa*x+coeffb  »weight.txt » u 1:2:(0.01) via coeffa,coeffb

où $range[0] et $range[1] sont les éléments d’un tableau du script perl qui donne l’intervalle. On donne ensuite la formule d’ajustement ; ici, une droite, rien de méchant. Le fichier en question en utilisant les colonnes 1 et 2 (abscisse et ordonnée) et une erreur sur l’ordonnée de 0.01 que j’ai estimé par rapport à l’appareil de mesure. Nous verrons plus tard son intérêt. Ce « fit » se fait via les deux coefficients en question.

On referme le fichier et on lance la danse de gnuplot.

Gnuplot produit en plus de la sortie dans la console, un fichier fit.log (valeur par défaut de la variable FIT_LOG) contenant toutes les informations du processus : une précieuse ressource. En voici un extrait :

After 9 iterations the fit converged.
final sum of squares of residuals : 97.494
rel. change during last iteration : -3.45394e-07

degrees of freedom    (FIT_NDF)                        : 18
rms of residuals      (FIT_STDFIT) = sqrt(WSSR/ndf)    : 2.3273
variance of residuals (reduced chisquare) = WSSR/ndf   : 5.41633

Final set of parameters            Asymptotic Standard Error
=======================            ==========================

coeffa          = 0.0930301        +/- 0.0009025    (0.9701%)
coeffb          = -128.167         +/- 2.924        (2.281%)

correlation matrix of the fit parameters:

coeffa coeffb
coeffa          1.000
coeffb         -1.000  1.000

J’ai enlevé les informations de base (heure, fichier, fonction utilisée…) ainsi que les étapes. Je reviendrai sur certaines de ces infos…

Vous l’aurez compris, la ligne désirée est :

coeffa          = 0.0930301        +/- 0.0009025    (0.9701%)

Je vais donc utiliser l’expression rationnelle (ou régulière) suivante :

if ($_ =~ s/(w+)s+=s(.*?)s++/-s+(.*?)s+((.*?)\%)/$1 $2 $3 $4/)

qui donne dans $_ la variable, la valeur, l’erreur et le pourcentage. Ce scalaire est « splité » avec les espaces ensuite… et j’écris ce que je veux dans un fichier.

Et voilà ! c’est gagné ! Simple non ?

Précisions concernant les ajustements

L’algorithme utilisé par gnuplot est celle de Levenberg-Marquardt. Méthode connue et reconnue, robuste, implémentée dans la plupart des logiciels. Loin d’être un parfait connaisseur, je vais mettre le doigt sur quelques points clefs lorsqu’on réalise un ajustement.

Pour ceux qui ne connaissent pas du tout les ajustements, l’objectif est de réduire le carré (pour des raisons de signe) de la différence entre la fonction est le point, divisé par le poids (incertitude) du point pour toute la courbe (somme).  C’est le « khi 2 » comme on le prononce en français.

Erreurs

On considère que les incertitudes ont une distribution gaussienne (pour des raisons de théorème central limite). Ainsi si la fonction passe dans les points à « la boite d’incertitude » près, alors la racine du chi 2 divisée par le nombre de point doit être de l’ordre de 1. Dans le fit.log :

rms of residuals      (FIT_STDFIT) = sqrt(WSSR/ndf)    : 2.3273

ce qui est réconfortant. Si la valeur était bien plus grande, alors l’ajustement ne serait sans doute pas très bon. Il faut revoir la fonction choisie. Si cette valeur était très petite, alors il est fort possible que l’incertitude soit sur évaluée. La conséquence pour une droite, c’est que la pente n’est pas bien définie car différentes pentes vont passées par les barres d’erreur.
Dans le cas où aucune erreur n’est spécifiée, elle est unitaire. Si c’est la même pour tous les points, on peut l’ajouter comme dit ci-dessus dans la commande fit. Sinon, on spécifie par une troisième colonne dans le fichier la valeur pour chaque point.

Degrés de liberté

Toujours dans le log du fit :

degrees of freedom    (FIT_NDF)                        : 18

C’est la différence entre le nombre de point de l’intervalle et le nombre de coefficient. Ce nombre doit être suffisamment grand bien sûr. Cette valeur est utile pour interpréter quantitativement un chi 2 à l’aide d’une table de chi 2 afin de connaître la probabilité d’arriver à un tel chi 2.

Valeurs initiales et non linéarité

Il est important, pour ne pas dire vital, de définir les valeurs des coefficients pour l’ajustement. Par défaut, ils sont égaux à 1. Pour notre régression affine, cela ne nous affectera pas le résultat, mais dans le cas d’une fonction non linéaire, il peut exister des minima locaux que l’on veut éviter. Parfois, l’ajustement peut se retrouver piégé dans un minimum non souhaité et le résultat est erroné. Il convient d’initialiser les variables comme coeffa=10.

Matrice de corrélation

Gnuplot affiche aussi la matrice de corrélation des paramètres. Elle peut s’avérer utile lorsque des ajustements non linéaires sont réalisés. En effet, on est amené à fixer des paramètres pour ne faire varier que les autres de manière à approcher au mieux l’extremum global que j’ai expliqué ci-dessus. Dans ce cas, on a une information du style : « si je fais varier le paramètre a, de combien cela influence-t-il le paramètre b,c,d… ». Cela peut aider dans le cas de fonction compliquée à élaborer une stratégie.

Limite d’arrêt

Vous verrez dans un fit.log la ligne suivante :
limit for stopping : 1e-05
Ceci correspond au critère de convergence du fit. Si la variation relative du WSSR est inférieure à cette limite, on considère que l’on a convergé. On peut modifier cette valeur comme ceci : FIT_LIMIT=1e-7.

Deux autres paramètres sont modifiables et permettent de modifier le comportement de l’algorithme. Cependant, je ne connais pas suffisamment de choses sur l’algorithme pour m’y aventurer, mais il est bon de savoir que c’est possible.

Conclusion

Voilà un article entre science et logiciel libre. Je me devais de livrer au moins une fois des louanges de gnuplot. Malgré les logiciels que l’on m’a présenté comme miraculeux, j’y suis toujours revenu pour sa performance, sa flexibilité, son interopérabilité… je m’arrête là Par ailleurs, vous voyez que les ajustements sont en fait tout un art loin d’être trivial.