====== Régression à l'aide de Perl et MetaPost ======
J'ai écrit (avec l'aide de Jean-Michel Sarlat) un script Perl permettant la création d'un fichier MetaPost où figure le tracé des points des données et le tracé de la régression. Le fichier MetaPost utilise //graph//.
===== La régression =====
L'avantage (pour moi) de ce script, est qu'il peut réaliser beaucoup de régressions. En effet, on précise l'ordre de régression (voir //Syntax//), c'est-à-dire le degré du polynôme qui est issu de la régression.
Par exemple l'ordre 3 donne une fonction du type : a+bx+cx^2+dx^3.
Il est bien entendu possible de réaliser une régression //exponentielle// (voir //Syntax//)
===== Utilisation =====
On doit disposer d'un fichier //monfichier.dat// avec les données dont on souhaite réaliser une régression. Le script crée un fichier MetaPost //MetaRgmonfichier.mp// "décoré".
J'entends par "décoré" qu'il figure des légendes et l'équation de la régression. Après, rien n'empêche de modifier le fichier MetaPost.
===== Syntax =====
Dans le cas général :
moindre.pl monfichier.dat [-d affichage des décimales] [légende ordonnées] [légende abscisses] [degrés de régression]
Le ''-d'' indique que c'est une option (non obligatoire).
==== Légendes ====
On rentre les légendes comme chaîne de caractère, i.e. entre guillemets :
moindre.pl monfichier.dat "Ordonnées" "Abscisses" 1
Ici, entre les guillemets, on utilise du code LaTeX, on peut donc utiliser des formules mathématiques!!!
==== Degrés de régression ====
Le degré de régression est, comme indiqué plus haut, le degré du polynôme "fonction régression".
=== Régression exponentielle ===
Pour réaliser une régression exponentielle, on tape ''exp'' à la place du nombre.
==== Option "nombre de décimale" ====
L'affichage de l'équation de la régression est par défaut à 3 décimales, si l'on veut changer cela on indique ''-d'' puis le nombre de décimale souhaité, ceci, //juste// après l'appel moindre.pl.
Exemple :
moindre.pl -d 6 monfichier.dat 2
==== Fichier coeff.dat ====
Un fichier coeff.dat est créé dans lequel sont inscrits les coefficients du polynôme.
==== Exemple ====
Voici un exemple présenté très joliement par Jean-Michel Sarlat :
[[http://melusine.eu.org/syracuse/journal/2008/07/21/perl/]]
==== Le code ====
#!/usr/bin/perl
# -----------------------------------------------------------------------------
# MOINDRE
# Maxime Chupin
# -----------------------------------------------------------------------------
# version 1.0
use Math::Matrix;
use Cwd;
use Getopt::Std; # Usage du module standard d'acquisition des options
# sur la ligne de commande.
#---------------- récupération du chemin --------------------------------------
my $chemin = cwd();
#------------------------------------------------------------------------------
#------------------- les options ----------------------------------------------
getopts("d:");
my $defaut_d = 3;
$opt_d and $defaut_d = $opt_d;
print ("$defaut_d\n");
#------------------------------------------------------------------------------
#---------------- subroutine d'un matrice unicolonne --------------------------
sub MatCol {
my @list = @_;
my $l = [];
foreach $v (@list){
push @$l, $v;
}
my $mat = new Math::Matrix($l);
return $mat->transpose;
}
#------------------------------------------------------------------------------
#---------------------- listes des points mesures -----------------------------
my @xi = ();
my @yi = ();
#------------------------------------------------------------------------------
#---------- récupération des listes de points d'un fichiers externe -----------
# Jean-Michel Sarlat, cette opération trie les listes par ordre croissant de x
@p = (); # Table de références à des points.
# Ouverture en lecture seule du fichier xy.dat
open(DAT,"$chemin/$ARGV[0]") or die "Impossible d'ouvrir $ARGV[0]... : $!";
while () { # Lecture ligne par ligne du fichier (handle DAT)
chomp; # Suppresion des blancs finaux, y compris le retour à la ligne
s/^\s*//; # Suppression des blancs initiaux (remplacement des blancs par rien...)
if($_){ # Traite la suite si la ligne est non vide
my $r = {}; # Création d'une référence vide...
($r->{x},$r->{y}) = split /\s+/; # Découpage de la ligne, le séparateur est
# une succession de blancs.
push @p, $r;} # Ajout de $r à la table des points.
}
# Tri de la liste des points selon x.
@p = sort { $a->{x} <=> $b->{x} } @p;
# Construction des tables des x et y.
foreach (@p) {
push @xi, $_->{x};
push @yi, $_->{y};
}
#------------------------------------------------------------------------------
#---------------------- min Max @xi et pas ------------------------------------
my $minxi;
my $maxxi;
$minxi = $xi[0];
$maxxi = $xi[-1];
$pas = abs($minxi-$maxxi)/50;
#------------------------------------------------------------------------------
#----------------------------- Les définitions des "outils" -------------------
my $taille = @xi; #comme son nom l'indique
my @S=(); # tableau des sommes des xi^k, k de 0 à 2*$ordre
my @W=(); # tableau des sommes des yi$xi^k, k de 0 à $ordre
my $ordre;
($prefixe = $ARGV[0]) =~ s/\.dat$//; #récupération du préfixe
if("$ARGV[1]" eq "exp"){
$ordre = 1;
}
else{
$ordre = $ARGV[1]; #ordre de regression
}
my $matrice = new Math::Matrix([1]) ;# définition de la matrice des Sk
#------------------------------------------------------------------------------
#-------------------- regression exponentielle --------------------------------
my @saveyi = @yi;
if("$ARGV[1]" eq "exp"){
for my $i(0..($taille-1)){
$yi[$i] = log($yi[$i]);
}
}
#------------------------ Calculs des S_k et W_k ------------------------------
for my $k ( 0 .. (2*$ordre) )
{
for my $i ( 0 .. ($taille-1) )
{
if($k==0){
$xik=1;
}
else{
$xik=1;
$pro=$xi[$i];
for my $l (1..$k){
$xik*=$pro;
}
}
$S[$k]+=$xik;
if($k<=$ordre){
$W[$k]+=$yi[$i]*$xik;
}
}
}
#------------------------------------------------------------------------------
$matW = MatCol(@W);
for my $i (0 .. $ordre){
for my $j (0 .. $ordre){
$matrice->[$i]->[$j] = $S[$i+$j];
}
}
$sys = $matrice->concat($matW);
$sol=$sys->solve;
#------------------------------------------------------------------------------
#------------------- liste des coordonnées de la régression -------------------
my @xr=();
my @yr=();
my $incr=0;
open(DATTH, ">$chemin/regr$ARGV[0]");
for(my $i=$minxi; $i<=$maxxi+$pas; $i+=$pas){
my $sum=0;
$xr[$incr]=$i;
#si regression exponentielle y=exp(ax+b)
if("$ARGV[1]" eq "exp"){
$yr[$incr] = exp(($sol->[0]->[0]))*(exp(($sol->[1]->[0])*$xr[$incr]));
}
#regression "normal"
else{
for my $j (0 .. $ordre+1){
$sum+=($sol->[$j]->[0])*(($xr[$incr]**($j)));
}
$yr[$incr]=$sum;
}
print (DATTH "$xr[$incr] $yr[$incr]\n");
$incr++;
}
close(DATTH);
#------------------------------------------------------------------------------
#------------------------ on imprime les coefficients -------------------------
open(FIC,">$chemin/coef$ARGV[0]");
print(FIC "Régression d'ordre $ARGV[1]\n $sol\n");
close(FIC);
#------------------------------------------------------------------------------
#----------------------- On crée le fichier metapost --------------------------
my $legende;
my $coeffexp;
open(FICMP, ">$chemin/MetaRg$prefixe.mp");
# Définition de la légende "regression"
# Définition de la légende "regression"
# JMS - Retouche pour montrer l'utilisation l'instruction sprintf. Elle vient du
# C et est très utile pour le formatage des résultats !
if("$ARGV[1]" eq "exp"){
$legende = sprintf("%0.${defaut_d}f \\text{e}^{%0.${defaut_d}f x}",
exp($sol->[0]->[0]),exp($sol->[1]->[0]));
}
else {
$legende = sprintf("%0.${defaut_d}f",$sol->[0]->[0]);
for my $i (1 .. $ordre){
$legende .= sprintf("%+0.${defaut_d}f x",$sol->[$i]->[0]);
$i > 1 and $legende .= "^{$i}";
}
}
# Calcul des coordonnées de "fixation" de la légende
my $xleg = $minxi+1*($maxxi-$minxi)/2;
# on trouve le $yi max
my $yleg;
if("$ARGV[1]" eq "exp"){
@yi = @saveyi;
}
for(my $i=0; $i<=$taille-1; $i+=1){
if($i==0){
$yleg=$yi[$i];
}
else{
if($yi[$i]>$yleg){
$yleg=$yi[$i];
}
}
}
print ("$xi[-1] $xleg $yleg \n");
print FICMP << "ENDOFMP";
verbatimtex
%&latex
\\documentclass{article}
\\usepackage[latin1]{inputenc}
\\usepackage{amsmath}
\\usepackage[frenchb]{babel}
\\begin{document}
etex
picture bullet;
bullet = image(draw origin withpen pencircle scaled 1.5mm withcolor green);
input graph;
beginfig(1);
draw begingraph(17cm,12cm);
glabel.lft(btex \\Large $ARGV[2] etex rotated 90,OUT);
glabel.bot(btex \\Large $ARGV[3] etex,OUT);
gdraw "$ARGV[0]" plot bullet ;
gdraw "regr$ARGV[0]" withcolor red withpen pencircle scaled 1.3pt;
autogrid(grid.bot,grid.lft) withcolor .85white ;
setrange(whatever,whatever,whatever,whatever);
picture legende;
legende = btex Régression : \$\\boxed{$legende}\$ etex;
glabel(image(unfill bbox legende; draw legende) , ($xleg,$yleg));
frame;
endgraph;
endfig;
end.
ENDOFMP
close(FICMP);
#------------------------------------------------------------------------------
#