Régression polynomiale

Soyez le premier à donner votre avis sur cette source.

Snippet vu 2 345 fois - Téléchargée 9 fois

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.
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.

A voir également

Ajouter un commentaire

Commentaires

cs_pseudo3
Messages postés
270
Date d'inscription
mardi 24 juillet 2007
Statut
Membre
Dernière intervention
7 juin 2018
-
Bonjour,

Juste une remarque au sujet des paramètres à transmettre à la procedure Regression2(var X,Y,Fit: array of Real; IndiceMax,DegRePolyn:Integer);
Inutile de transmettre IndiceMax car IndiceMax:= High(X) ou High(Y) ou High(Fit)

Ce serait génial si en plus on pouvait ne pas être obligé de transmettre DegrePolyn de sorte que la procédure détermine elle-même le degré optimal du polygone qui passe le plus près des points expérimentaux.

A+.
cs_Gerard
Messages postés
123
Date d'inscription
jeudi 10 janvier 2002
Statut
Membre
Dernière intervention
7 août 2018
> cs_pseudo3
Messages postés
270
Date d'inscription
mardi 24 juillet 2007
Statut
Membre
Dernière intervention
7 juin 2018
-
La remarque est très juste! Merci de l'avoir suggérée!
Par contre pour le degré du polynôme la réponse est évidente: il s'agit du nombre de points expérimentaux moins 1 me semble-t-il...
> cs_pseudo3
Messages postés
270
Date d'inscription
mardi 24 juillet 2007
Statut
Membre
Dernière intervention
7 juin 2018
-
Bonjour,

Ayant testé le code j'ai constaté un comportement bizarre : Tant que le degré du polynôme est inférieur à 11 la courbe se trace bien entre les points expérimentaux,mais à partir de 11 elle est presque entièrement en dehors des clous : quelle en est la cause ???

A+.

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.