Discussion :

[ford_escort]

Salut a tous,
je souhaiiterai implementer un algorithme pour faire du traitements de donnees.
J aurai besoin de fitter une courbe de la forme a*x(t)^3 +b*x(t)^2 = t
Afin d avoir les valeurs de a et b pour chaque courbes.
Merci
Ford

[méphistopheles]

fitter

c'est un mot qui existe ou c'est une faute de frappe?

[ford_escort]

C est un anglicisme, je ne sais pas trop si il y a un ou deux "t"......

[FrancisSourd]

Il y a par exemple la fonction fit de gnuplot.

[ford_escort]

Je sais que gnuplot le fait (ou xmgrace voire meme excel on ne sait jamais)
Mais je doit traiter des tonnes de fichiers de data, donc tous les entrer a la main serait une vraie perte de temps.
Vu que pas mal de soft le font, je me suis dit que c etait peut etre un truc classique et que l algo devait etre connu

[Matthieu Brucher]

Tout dépend de la méthode que tu veux utiliser, mais les moindres carrés de base semblent bien adaptés.

[j.p.mignot]

un petit code en pascal pour répondre au problème
cela est un moindre carrés non itératif car pour ce fit on peut formuler le résidu à minimiser se façon à ce que les coefs recherchés soient solution d'un système de kramer.
Bien entendu un fit moindres carrés itératif général fonctionne. J'en ai aussi donné un source code en pascal et en C dans un de mes précédents posts ( novembre 2005 )

ici on prend p points PX(a),PY(a) a=1..p et un polynôme de degré n
y = somme(i=0..n) Ci X^i

on nomme Ra (a=1..p) l'écart somme(i=0..n) Ci * PX(a)^i - Y(a)

et on minimise S = somme (a=1..p) (Ra^2) en résolvant le système de n+1 équations à n+1 inconnues défini par @S/@Ck = 0 pour tout k de 0 à n soit n+1 équations
L'implémentation proposée propose une solution par triangularisation
du système

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
 
 
{$n+,g+,a+}
unit math_rx;
interface
const     deg_mx     = 20;               { degre max du fit }
          data_mx    = 16383;             { nombre max de data }
 
type
          a0         = array[0..deg_mx,0..deg_mx] of extended;
          b0         = array[0..deg_mx]      of extended;
          xy         = array[1..data_mx]     of single;
 
function  polynome_fits(
                        px,py : pointer;  { pointe les tableaux de x,y de type single dont
                                            le nombre d'element ne doit pas
                                            exceder data_mx }
                        pf : pointer;     { pointe 1 tableau 0..dg of extended pour retour coef }
                        dg,ndat : byte;   { degre ( doit etre <= dg_mx) et nombre de data (doit etre <= ndata_mx )}
                        var ki2 : single; { pour retour du residu }
                        _weight : pointer) :   { pointeur sur tableau de poids (NIL si non souhait‚)}
                        boolean;
implementation
 
const     epsilon1   = 1e-25; { arbitrairement petit }
 
var
    a_mat     : a0;
    arot,
    piv,
    b_mat     : b0;
    brot      : extended;
    function puiss0(s : extended; v : byte) : extended;
    var se : extended;
       begin
       if v=0 then
          puiss0:=1
       else
       if abs(s) < 1e-40 then
          puiss0:=0
       else
          begin
          se:=exp(v*ln(abs(s)));
          if (se < 0) then if ((v and 1) = 1) then se:=-se;
          puiss0:=se;
          end;
       end;
function polynome_fits( px,py : pointer; pf : pointer ; dg,ndat : byte;
                      var ki2 : single ; _weight : pointer) : boolean;
   procedure make_a_b;
   var i,j,k : integer;
      begin
      if _weight=nil then
         begin
         for i:=0 to dg do
            begin
            b_mat[i]:=0;
            for k:=1 to ndat do
               b_mat[i]:=b_mat[i] + xy(py^)[k] *
                         puiss0( xy(px^)[k],i);
            for j:=0 to dg do
               begin
               a_mat[i,j]:=0;
               for k:=1 to ndat do
                  a_mat[i,j]:=a_mat[i,j] +
                         puiss0( xy(px^)[k],i+j)
               end
            end
         end
      else
         begin
         for i:=0 to dg do
            begin
            b_mat[i]:=0;
            for k:=1 to ndat do
               begin
               b_mat[i]:=b_mat[i] +
                      xy(py^)[k] * puiss0( xy(px^)[k],i) *
                      xy(_weight^)[k];
               end;
            for j:=0 to dg do
               begin
               a_mat[i,j]:=0;
               for k:=1 to ndat do
                  begin
                  a_mat[i,j]:=a_mat[i,j] +
                       puiss0(xy(px^)[k],i+j)  *
                       xy(_weight^)[k];
                  end;
               end
            end
         end
      end;
   var toto : boolean;
       ktest : integer;
   procedure permut (i,ns : integer);
   var j,k : integer;
      begin
      inc(ktest);
      brot:=b_mat[i];
      for j:=i to dg do arot[j]:=a_mat[i,j];
      for j:=i to dg-1 do
         begin
         for k:=i to dg do a_mat[j,k]:=a_mat[j+1,k];
         b_mat[j]:=b_mat[j+1]
         end;
       for k:=i to dg do
          a_mat[dg,k]:=arot[k];
       b_mat[dg]:=brot;
       toto:=false
       end;
   function triangle : boolean;
   var i,j,k : integer;
      begin
      triangle:=false;
      for i:=0 to dg-1 do
         begin
         ktest:=0;
            repeat
            toto:=true;
            if ktest > dg then exit;
            if a_mat[i,i]=0 then permut(i,dg)
            until toto;
         for j:=i+1 to dg do piv[j]:=a_mat[i,j]/a_mat[i,i];
         for j:=i+1 to dg do
            begin
            for k:=i to dg do
               a_mat[j,k]:=a_mat[j,k]-piv[j]*a_mat[i,k];
            b_mat[j]:=b_mat[j]-piv[j]*b_mat[i]
            end
         end;
      triangle:=true
      end;
   procedure solu;
   var u   : extended;
       i,j : integer;
      begin
      for i:=dg downto 0 do
         begin
         u:=b_mat[i];
         if i < dg then
            for j:=i+1 to dg do u:=u-a_mat[i,j] * b0(pf^)[j];
         u:=u/a_mat[i,i];
         b0(pf^)[i]:=u;
         end
      end;
   var s,n : extended;
       i,j : integer;
   begin
   polynome_fits:=false;
   if dg > deg_mx  then exit;
   if ndat <= dg   then exit;
   make_a_b;
   if not triangle then exit;
   solu;
   ki2:=0;
   n:=0;
   for i:=1 to ndat do
      begin
      s:=0;
      for j:=0 to dg do
         s := s + puiss0(xy(px^)[i],j) * b0(pf^)[j];
      n:=n + sqr(xy(py^)[i]);
      ki2:=ki2 + sqr(xy(py^)[i]-s)
      end;
   if n > 1e-100 then
      ki2:=sqrt(ki2/n)*100
   else
      ki2:=1e35;
   polynome_fits:=true
   end;
end.