Le wiki de Syracuse

====== 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 : <code bash> moindre.pl monfichier.dat [-d affichage des décimales] [légende ordonnées] [légende abscisses] [degrés de régression] </code> 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 : <code bash> moindre.pl monfichier.dat "Ordonnées" "Abscisses" 1 </code> 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 : <code bash> moindre.pl -d 6 monfichier.dat 2 </code> ==== 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 ==== <code perl> #!/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 (<DAT>) { # 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); #------------------------------------------------------------------------------ # </code>