Histogrammes barres et/ou camembert

=========================================================================
Curve fit: y = b0 + b1.x + b2.x^2 + b3.x^3
-------------------------------------------------------------------------
Parameter  Est.value        Std.dev.        95% Confidence Interval
-------------------------------------------------------------------------
  b0    481.78951172      25.69367856    423.66637272;    539.91265071
  b1    -938.48356383      47.82386598  -1046.66866480;    -830.29846285
  b2    827.47418099      93.69499502    615.52137688;    1039.42698510
  b3    368.67153507      52.94578940    248.89983833;    488.44323181
-------------------------------------------------------------------------
Number of observations            : n          =    13
Residual error                    : s          =    57.4194
Coefficient of correlation        : r          =    0.9984
Coefficient of determination      : r2          =    0.9968
Adjusted coeff. of determination  : r2a        =    0.9958
Variance ratio (explained/resid.) : F(  3,  9) =  938.2036
Critical variance ratio          : F(p = 0.95) =    3.8625
-------------------------------------------------------------------------
  i        Y obs.      Y calc.      Residual      Std.dev.      Std.res.
-------------------------------------------------------------------------
  1    2615.0000    2657.3287      -42.3287      57.4194      -0.7372
  2    2474.0000    2405.8812      68.1188      57.4194        1.1863
  3    2152.0000    2188.9006      -36.9006      57.4194      -0.6426
  4    2114.0000    2024.6634      89.3366      57.4194        1.5559
  5    1862.0000    1879.0757      -17.0757      57.4194      -0.2974
  6    1364.0000    1417.7360      -53.7360      57.4194      -0.9359
  7      910.0000      965.9838      -55.9838      57.4194      -0.9750
  8      702.0000      724.3800      -22.3800      57.4194      -0.3898
  9      586.0000      583.5439        2.4561      57.4194        0.0428
 10      501.0000      481.7895      19.2105      57.4194        0.3346
 11      392.0000      341.8227      50.1773      57.4194        0.8739
 12      330.0000      284.6713      45.3287      57.4194        0.7894
 13      250.0000      296.2231      -46.2231      57.4194      -0.8050
=========================================================================

{ ******************************************************************
  Modeles de courbes de tendance (version Delphi)
  ------------------------------------------------------------------
  Necessite la bibliotheque DMath :
  http://www.unilim.fr/pages_perso/jean.debord/tpmath/dmat086.zip
  ****************************************************************** }

unit umodels;

interface

uses
  { Unites de la bibliotheque DMath }
  utypes, upolynom, ulinfit, upolfit, unlfit, ufracfit, uexpfit,
  ulogifit, upowfit, uhillfit, uevalfit, uregtest, uinvbeta;

type TModel =
  (REG_LIN,    { y = a + b.x }
  REG_POL2,  { y = a + b.x + c.x^2 }
  REG_POL3,  { y = a + b.x + c.x^2 + d.x^3 }
  REG_FRAC1,  { y = a.x / (1 + b.x) }
  REG_FRAC2,  { y = (a + b.x) / (1 + c.x) }
  REG_EXPO1,  { y = a.exp(-b.x) }
  REG_EXPO2,  { y = a + b.exp(-c.x) }
  REG_LOGIS1, { y = a / [1 + exp(-b.x + c)] }
  REG_LOGIS2, { y = a + (b - a) / [1 + exp(-c.x + d)] }
  REG_POWER1, { y = a.x^b }
  REG_POWER2, { y = a / [1 + (b / x)^c] }
  REG_USER);  { Fonction definie par l'utilisateur }

function FirstParam(Model : TModel) : Integer;
{ Retourne l'indice du premier parametre du modele }

function LastParam(Model : TModel) : Integer;
{ Retourne l'indice du dernier parametre du modele }

function FuncName(Model : TModel) : String;
{ Retourne le nom de la fonction de regression }

function ParamName(Model : TModel; I : Integer) : String;
{ Retourne le nom du parametre d'indice I }

procedure FitModel(Model      : TModel;
                  X, Y, Ycalc : TVector;
                  N          : Integer;
                  B          : TVector;
                  V          : TMatrix;
                  var Test    : TRegTest);
{ Ajustement du modele
  Entree : Model : indice du modele de regression
          X, Y  : coordonnees des points
          N    : nombre de points
  Sortie : Ycalc : valeurs calculees de Y
          B    : parametres du modele
          V    : matrice de variance-covariance
          Test  : tests statistiques }   

function RegFunc(Model : TModel; X : Float; B : TVector) : Float;
{ Calcul de la fonction de regression au point X
  (a appeler apres FitModel) }

procedure WriteResults(FileName : String;
                      Model    : TModel;
                      X, Y,
                      Ycalc    : TVector;
                      N        : Integer;
                      B        : TVector;
                      V        : TMatrix;
                      Test    : TRegTest);
{ Ecrit les resultats dans le fichier de sortie }

implementation

const
  MaxIter = 1000;    { Nb max d'iterations pour regression non lineaire }
  Tol    = 1.0E-10; { Precision des parametres de la reg. non lineaire }
  SVDTol  = 1.0E-15; { Tolerance pour la decomposition en val. singulieres }
  Alpha  = 0.05;    { Niveau de signification des tests statistiques }

function FirstParam(Model : TModel) : Integer;
begin
  case Model of
    REG_LIN    : FirstParam := 0;
    REG_POL2  : FirstParam := 0;
    REG_POL3  : FirstParam := 0;
    REG_FRAC1  : FirstParam := 1;
    REG_FRAC2  : FirstParam := 0;
    REG_EXPO1  : FirstParam := 1;
    REG_EXPO2  : FirstParam := 0;
    REG_LOGIS1 : FirstParam := 1;
    REG_LOGIS2 : FirstParam := 0;
    REG_POWER1 : FirstParam := 0;
    REG_POWER2 : FirstParam := 1;
    REG_USER  : FirstParam := 1;
  end;
end;

function LastParam(Model : TModel) : Integer;
begin
  case Model of
    REG_LIN    : LastParam := 1;
    REG_POL2  : LastParam := 2;
    REG_POL3  : LastParam := 3;
    REG_FRAC1  : LastParam := 2;
    REG_FRAC2  : LastParam := 2;
    REG_EXPO1  : LastParam := 2;
    REG_EXPO2  : LastParam := 2;
    REG_LOGIS1 : LastParam := 3;
    REG_LOGIS2 : LastParam := 3;
    REG_POWER1 : LastParam := 1;
    REG_POWER2 : LastParam := 3;
    REG_USER  : LastParam := uevalfit.LastParam;
  end;
end;

function FuncName(Model : TModel) : String;
begin
  case Model of
    REG_LIN    : FuncName := 'y = a + b.x';
    REG_POL2  : FuncName := 'y = a + b.x + c.x^2';
    REG_POL3  : FuncName := 'y = a + b.x + c.x^2 + d.x^3';
    REG_FRAC1  : FuncName := 'y = a.x / (1 + b.x)';
    REG_FRAC2  : FuncName := 'y = (a + b.x) / (1 + c.x)';
    REG_EXPO1  : FuncName := 'y = a.exp(-b.x)';
    REG_EXPO2  : FuncName := 'y = a.exp(-b.x) + c';
    REG_LOGIS1 : FuncName := 'y = a / [1 + exp(-b.x + c)]';
    REG_LOGIS2 : FuncName := 'y = a + (b - a) / [1 + exp(-c.x + d)]';
    REG_POWER1 : FuncName := 'y = a.x^b';
    REG_POWER2 : FuncName := 'y = a / [1 + (b / x)^c]';
    REG_USER  : FuncName := 'y = ' + uevalfit.FuncName;
  end;
end;

function ParamName(Model : TModel; I : Integer) : String;
begin
  if Model = REG_USER then
    ParamName := uevalfit.ParamName(I)
  else 
    ParamName := Chr(97 + I - FirstParam(Model));
end;

function RegFunc(Model : TModel; X : Float; B : TVector) : Float;
begin
  case Model of
    REG_LIN    : RegFunc := B[0] + B[1] * X;
    REG_POL2  : RegFunc := Poly(X, B, 2);
    REG_POL3  : RegFunc := Poly(X, B, 3);
    REG_FRAC1,
    REG_FRAC2  : RegFunc := FracFit_Func(X, B);
    REG_EXPO1,
    REG_EXPO2  : RegFunc := ExpFit_Func(X, B);
    REG_LOGIS1,
    REG_LOGIS2 : RegFunc := LogiFit_Func(X, B);
    REG_POWER1 : RegFunc := PowFit_Func(X, B);
    REG_POWER2 : RegFunc := HillFit_Func(X, B);
    REG_USER  : RegFunc := EvalFit_Func(X, B);
  end;
end;

procedure FitModel(Model      : TModel;
                  X, Y, Ycalc : TVector;
                  N          : Integer;
                  B          : TVector;
                  V          : TMatrix;
                  var Test    : TRegTest);
var
  I : Integer;

begin
  { Ajuster le modele }
  case Model of
    REG_LIN    : SVDLinFit(X, Y, 1, N, N * SVDTol, B, V);
    REG_POL2  : SVDPolFit(X, Y, 1, N, 2, N * SVDTol, B, V);
    REG_POL3  : SVDPolFit(X, Y, 1, N, 3, N * SVDTol, B, V);
    REG_FRAC1  : FracFit(X, Y, 1, N, 1, 1, False, MaxIter, Tol, B, V);
    REG_FRAC2  : FracFit(X, Y, 1, N, 1, 1, True, MaxIter, Tol, B, V);
    REG_EXPO1  : ExpFit(X, Y, 1, N, 1, False, MaxIter, Tol, B, V);
    REG_EXPO2  : ExpFit(X, Y, 1, N, 1, True, MaxIter, Tol, B, V);
    REG_LOGIS1 : LogiFit(X, Y, 1, N, False, False, MaxIter, Tol, B, V);
    REG_LOGIS2 : LogiFit(X, Y, 1, N, True, False, MaxIter, Tol, B, V);
    REG_POWER1 : PowFit(X, Y, 1, N, MaxIter, Tol, B, V);
    REG_POWER2 : HillFit(X, Y, 1, N, False, MaxIter, Tol, B, V);
    REG_USER  : EvalFit(X, Y, 1, N, MaxIter, Tol, B, V);
  end;

  { Calculer les Y estimes }
  for I := 1 to N do
    Ycalc[I] := RegFunc(Model, X[I], B);

  { Mettre a jour la matrice de variance-covariance
    et calculer les tests statistiques }
  RegTest(Y, Ycalc, 1, N, V, FirstParam(Model), LastParam(Model), Test);
end;

procedure WriteResults(FileName : String;
                      Model    : TModel;
                      X, Y,
                      Ycalc    : TVector;
                      N        : Integer;
                      B        : TVector;
                      V        : TMatrix;
                      Test    : TRegTest);
const
  Line1 = '-------------------------------------------------------------------------';
  Line2 = '=========================================================================';

var
  OutF  : Text;    { Fichier de sortie }
  Delta : Float;    { Residu }
  Sr    : Float;    { Ecart-type residuel }
  SB    : Float;    { Ecart-type d'un parametre }
  R    : Float;    { Coefficient de correlation }
  Tc    : Float;    { Valeur critique de t au risque Alpha }
  Fc    : Float;    { Valeur critique de F au risque Alpha }
  I    : Integer;  { Variable de Boucle }

begin
  Assign(OutF, FileName);
  Rewrite(OutF);
 
  Writeln(OutF, Line2);
  Writeln(OutF, 'Ajustement de la fonction : ', FuncName(Model));
  Writeln(OutF, Line1);

  Writeln(OutF, 'Param.    Valeur        Ecart-type          ',
          'Intervalle de confiance ', (100 * (1 - Alpha)):2:0, '%');

  Writeln(OutF, Line1);
 
  { Calculer le t de Student et le F de Snedecor }
  Tc := InvStudent(Test.Nu2, 1 - 0.5 * Alpha);
  Fc := InvSnedecor(Test.Nu1, Test.Nu2, 1 - Alpha);

  for I := FirstParam(Model) to LastParam(Model) do
    begin
      SB := Sqrt(V[I,I]);
      Writeln(OutF, ParamName(Model, I), B[I]:17:8, SB:17:8,
              (B[I] - Tc * SB):20:8, ';', (B[I] + Tc * SB):17:8);
    end;

  Writeln(OutF, Line1);

  Writeln(OutF, 'Nombre d''observations        : n          = ', N:5);

  with Test do
    begin
      Sr := Sqrt(Vr);
      Writeln(OutF, 'Erreur residuelle            : s          = ', Sr:10:4);
      if R2 <= 1.0 then
        begin
          R := Sqrt(R2);
          if (Model = REG_LIN) and (B[1] < 0.0) then R := -R;
          Writeln(OutF, 'Coefficient de correlation  : r          = ', R:10:4);
          Writeln(OutF, 'Coefficient de determination : r2          = ', R2:10:4);
        end; 
      Writeln(OutF, 'Rapport de variances        : F(', Nu1:3, ', ', Nu2:3, ') = ', F:10:4);
      Writeln(OutF, 'Valeur critique de F        : F(p = ', (1 - Alpha):4:2, ') = ', Fc:10:4);
    end;

  Writeln(OutF, Line1);
  Writeln(OutF, '  i        Y obs.      Y calc.        Residu        E.type    Res.norm.');
  Writeln(OutF, Line1);

  for I := 1 to N do
    begin
      Delta := Y[I] - Ycalc[I];
      Writeln(OutF, I:3, Y[I]:14:4, Ycalc[I]:14:4, Delta:14:4, Sr:14:4, (Delta / Sr):14:4);
    end;

  Writeln(OutF, Line2);
 
  Close(OutF);
end;

end.