Régression polynomiale
Contenu du snippet
En m'inspirant largement du travail de cs_Gimli j'ai transformé son code de façon à ce que ce soit une procédure qui renvoie directement la courbe.
C'est ainis plus facile à appeler dans un projet.
Le tableau dynamique Fit doit être dimensionné avant l'appel.
Si l'un des tableaux est d'un type différent (integer, TDateTime...) il suffit de modifier la séquence d'appel.
C'est ainis plus facile à appeler dans un projet.
Le tableau dynamique Fit doit être dimensionné avant l'appel.
procedure Regression2(var X,Y,Fit: array of Real; IndiceMax,DegRePolyn:Integer);
var
//Le tableau Fit doit être dimensionné avant l'appel à cette fonction
a : Array of Double; // Coefficients polynomiaux cherchés
c : Array of Array of Double; // Coefficients de Cramer
b : Array of Double; // Vecteur constant du système de Cramer
i, j, k : Integer; // variables de boucle
d : Double;
t : Double;
Pt : TPoint;
{ Calcul du déterminant par la méthode de Gauss }
function Determinant : Double; // Calcul le déterminant de la matrice c
var
c2 : Array of Array of Double;
l, m, q, r : Integer;
s : Double;
begin
SetLength(c2, DegrePolyn + 1, DegrePolyn + 1);
for l := 0 to DegrePolyn do for m := 0 to DegrePolyn do c2[l, m] := c[l, m];
for l := 0 to DegrePolyn do
begin
m := l;
while (c2[l, m] = 0) and (m < DegrePolyn) do inc(m);
if c2[DegrePolyn, l] = 0 then
begin
Result := 0;
Exit;
end
else begin
{ On échange la colonne l et la colonne m }
for q := 0 to DegrePolyn do
begin
t := c2[q, l];
c2[q, l] := c2[q, m];
c2[q, m] := t;
end;
{ Pour les colonnes suivantes : Cq <- Cq - (aq / am) * Cm }
for q := l + 1 to DegrePolyn do
begin
s := c2[l, q] / c2[l, l];
for r := 0 to DegrePolyn do c2[r, q] := c2[r, q] - s * c2[r, l]
end;
end;
end;
Result := 1;
for l := 0 to DegrePolyn do Result := Result * c2[l, l];
end;
begin
{ La régression polynomiale consiste à minimiser la distance d'un polynome P
à des points expérimentaux }
{ Si on note (x1, y1),..., (xn, yn) ces points, x = (x1,...,xn), y = (y1,...,yn)
alors cela revient à minimiser ||y - P(x)|| (norme 2 euclidienne) donc P(x)
doit être le projeté orthogonal de y sur le sous-espace engendré par les
puissances successives de x.
La nullité des produits scalaires (y - P(x) | power(x, i) ) = 0 nous donne
système de Cramer cherché }
{ Initialisation }
SetLength(a, DegrePolyn + 1);
{ La résolution de la régression polynomiale conduit à 1 système de Cramer }
{ Calcul des coefficients du système de Cramer }
SetLength(c, DegrePolyn + 1, DegrePolyn + 1);
SetLength(b, DegrePolyn + 1);
for i := 0 to DegrePolyn do
begin
for j := 0 to DegrePolyn do
begin
c[i, j] := 0;
for k := 0 to IndiceMax do c[i, j] := c[i, j] + Power(x[k], i + j);
end;
end;
for i := 0 to DegrePolyn do
begin
b[i] := 0;
for j := 0 to IndiceMax do b[i] := b[i] + y[j] * Power(x[j], i);
end;
{ Résolution du système de Cramer }
d := Determinant;
for j := 0 to DegrePolyn do
begin
{ On permute la colonne j et la colonne b }
for i := 0 to DegrePolyn do
begin
t := b[i];
b[i] := c[i, j];
c[i, j] := t;
end;
a[j] := Determinant / d;
{ On re-permute la colonne j et la colonne b }
for i := 0 to DegrePolyn do
begin
t := b[i];
b[i] := c[i, j];
c[i, j] := t;
end;
end;
{ On calcule la fonction polynomiale trouvée }
for i := 0 to IndiceMax do
begin
Fit[I] := 0;
for j := 0 to DegrePolyn do Fit[I] := Fit[I] + a[j] * Power(X[I], j);
end;
end;
end.
Si l'un des tableaux est d'un type différent (integer, TDateTime...) il suffit de modifier la séquence d'appel.
Vous n'êtes pas encore membre ?
inscrivez-vous, c'est gratuit et ça prend moins d'une minute !
Les membres obtiennent plus de réponses que les utilisateurs anonymes.
Le fait d'être membre vous permet d'avoir un suivi détaillé de vos demandes et codes sources.
Le fait d'être membre vous permet d'avoir des options supplémentaires.