Outils d'utilisateurs

Outils du Site


mc:moindre.pl

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 (<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);
#------------------------------------------------------------------------------
#
mc/moindre.pl.txt · Dernière modification: 2010/02/12 02:22 par newacct