Régression non-linéaire polynomiale par la méthode des moindres carrés

Soyez le premier à donner votre avis sur cette source.

Vue 22 638 fois - Téléchargée 1 088 fois

Description

N'ayant pas trouvé de code simple traitant la régression non-linéaire par des polynômes, je me suis décidé à la programmer moi-même.

Ce code me paraît simple à reprendre.
Il suffit d'indiquer les valeurs x et y sur lesquelles porte la régression dans un fichier nommé "valeurs_init.txt". Ce fichier doit être placé dans le même répertoire que le fichier .bas

Ensuite, il fait lancer la procédure main.
Par défaut le polynôme de tendance est de degré 6, il est possible d'en changer en indiquant le degré voulu lors de l'appel de la procédure régression.

Le code génère alors les valeurs x et y = a0 + a1*x + ... + a6*x^6

Source / Exemple :


Option Explicit

Public x() As Double                'Données d'entrée (xi et yi)
Public y() As Double
Public ai_p() As Double             'Valeurs des coefficients du polynôme de tendance

Sub main()
Dim i As Integer, j As Integer
Dim ligne As String
Dim virg_deb As Long
Dim virg_fin As Long

Dim a As String
    
    i = 0
    'Boucle de récupération des données xi et yi à traiter
    Open "valeurs_init.txt" For Input As #1
    While Not EOF(1)
        i = i + 1
        ReDim Preserve x(1 To i)
        ReDim Preserve y(1 To i)
        
        Line Input #1, ligne
        
        'Remplacement des blancs par des séparateurs (chr(47) = "/")
        ligne = Replace(ligne, Chr(9), Chr(47))
        ligne = Replace(ligne, Chr(32), Chr(47))
        ligne = Chr(47) & ligne & Chr(47)
        
        'Récupération des données x(i) et y(i)
        'Recherche des séparateurs (/) et stockage des valeurs indiquées entre les séparateurs
        virg_deb = InStr(1, ligne, Chr(47))
        virg_fin = InStr((virg_deb + 1), ligne, Chr(47))
        x(i) = CDbl(Mid(ligne, (virg_deb + 1), (virg_fin - virg_deb - 1)))
        
        virg_deb = InStr((virg_fin), ligne, Chr(47))
        virg_fin = InStr((virg_deb + 1), ligne, Chr(47))
        y(i) = CDbl(Mid(ligne, (virg_deb + 1), (virg_fin - virg_deb - 1)))
    Wend
    
    Close #1
    
    'Appel de la procédure de régression polynômiale :
    'le premier argument est le degré du polynôme de tendance chosi, x() et y() sont les données à traiter
    'ai_p() est un argument de retour renvoyant les coef ai du polynôme
    Call regression_polynomiale(6, x(), y(), ai_p())
    
    
    Dim xf() As Double
    Dim yf() As Double
    ReDim xf(1 To UBound(x))
    ReDim yf(1 To UBound(y))
    
    'Boucle générant un fichier texte contenant les valeurs xf(i) et yf(i) approchées
    Open "valeurs_fin.txt" For Output As #1
    For i = 1 To UBound(x)
        xf(i) = x(i)
        yf(i) = ai_p(0)
        For j = 1 To (UBound(ai_p()))
            yf(i) = yf(i) + ai_p(j) * xf(i) ^ j
        Next j
        
        Print #1, xf(i) & Chr(9) & Format(yf(i), "####0.000")
    Next i
    
End Sub

Sub regression_polynomiale(deg_pol As Integer, xi() As Double, yi() As Double, coef_pol() As Double)
Dim i As Integer, j As Integer
Dim k As Integer, l As Integer
Dim nb_pts As Integer
Dim mat_A() As Double, mat_B() As Double
Dim mat_At() As Double, mat_Bt() As Double
Dim dim_mat As Integer

dim_mat = deg_pol + 1

'Nombre de points
nb_pts = UBound(xi)
If nb_pts <> UBound(yi) Then MsgBox "Problème, il n'y a pas le même nombre de valeurs pour les x et les y"

ReDim mat_A(1 To dim_mat, 1 To dim_mat)
ReDim mat_B(1 To dim_mat)
ReDim mat_At(1 To dim_mat, 1 To dim_mat)
ReDim mat_Bt(1 To dim_mat)
ReDim coef_pol(0 To deg_pol)

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
''           Etape 1 :
'' Mise en équation du problème de minimisation par la méthode des moindres carrés
'' on cherche à minimiser : E(a0,..., ap) = somme(yi-somme(aj*xi^j)²) où p est le degré du pôlynome
'' les dérivées partielles de E par rapport à chacun des coef du polynôme de tendance doivent être nulles,
'' c'est une condition nécessaire !
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' les xi et yi sont transmis en paramètres
' la procédure renvoie alors en sortie 2 matrices mat_A et mat_B tq :
' mat_A * vect_des_aj = mat_B
Call mise_en_éq(nb_pts, deg_pol, xi(), yi(), mat_A(), mat_B())

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
''           Etape 2 :
'' On réalise alors la triangulation de la mat_A et les opérations équivalentes sur mat_B
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' mat_A et mat_B sont données en arguments d'entrée
' mat_At et mat_Bt sont renvoyées en arguments de sortie, après triangulation de mat_A
Call triangulation(dim_mat, mat_A(), mat_B(), mat_At(), mat_Bt())

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
''           Etape 3 :
'' Cette procédure résout alors le système et renvoie les coefficients aj du pôlynome de tendance
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'les matrices mat_At et mat_Bt sont données en arguments d'entrée
'les coef aj sont renvoyées en arguments de sortie sous le nom de variable coef_pol
Call resolution(dim_mat, mat_At(), mat_Bt(), coef_pol())

End Sub

Sub mise_en_éq(n As Integer, p As Integer, vect1() As Double, vect2() As Double, mat1() As Double, mat2() As Double)

Dim i As Integer, j As Integer, k As Integer

For k = 0 To p          ' on balaie les lignes
    For j = 0 To p      ' on balaie les colonnes
        mat1(k + 1, j + 1) = 0
        For i = 1 To n  ' on somme tous les xi^(k+j)
            mat1(k + 1, j + 1) = mat1(k + 1, j + 1) + vect1(i) ^ (k + j)
        Next i
    Next j
Next k

For k = 0 To p          ' on balaie les lignes
        mat2(k + 1) = 0
    For i = 1 To n      ' on somme tous les (yi*xi^(k))
        mat2(k + 1) = mat2(k + 1) + (vect2(i) * vect1(i) ^ (k))
    Next i
Next k

End Sub

Sub triangulation(n As Integer, mat1() As Double, mat2() As Double, mat3() As Double, mat4() As Double)

Dim i As Integer, j As Integer, k As Integer
Dim mat3k() As Double, mat4k() As Double
ReDim mat3k(1 To n, 1 To n, 1 To n) As Double
ReDim mat4k(1 To n, 1 To n) As Double

For i = 1 To n
    mat4k(i, 1) = mat2(i)
    For j = 1 To n
        mat3k(i, j, 1) = mat1(i, j)
    Next j
Next i

For k = 1 To (n - 1)
    For i = 1 To n
        For j = 1 To n
            If i <= k Then
                mat3k(i, j, (k + 1)) = mat3k(i, j, k) ' on conserve les mêmes valeurs
                mat4k(i, (k + 1)) = mat4k(i, k)
            Else
                mat4k(i, (k + 1)) = mat4k(i, k) - (mat3k(i, k, k) * mat4k(k, k) / mat3k(k, k, k))
                If j <= k Then
                    mat3k(i, j, (k + 1)) = 0
                Else
                    mat3k(i, j, (k + 1)) = mat3k(i, j, k) - (mat3k(i, k, k) * mat3k(k, j, k) / mat3k(k, k, k))
                End If
            End If
        Next j
    Next i
Next k

For i = 1 To n
    mat4(i) = mat4k(i, n)
    For j = 1 To n
        mat3(i, j) = mat3k(i, j, n)
    Next j
Next i
End Sub

Sub resolution(n As Integer, mat3() As Double, mat4() As Double, coef_pol() As Double)

Dim i As Integer, j As Integer, k As Integer
Dim xsol() As Double
Dim som_aijxj As Double
ReDim xsol(1 To n)

xsol(n) = mat4(n) / mat3(n, n)
coef_pol(n - 1) = xsol(n)

For i = (n - 1) To 1 Step (-1)
    som_aijxj = 0
    For j = i + 1 To n
        som_aijxj = som_aijxj + mat3(i, j) * xsol(j)
    Next j
    xsol(i) = 1 / mat3(i, i) * (mat4(i) - som_aijxj)
    coef_pol(i - 1) = xsol(i)
Next i

End Sub

Conclusion :


J'espère que ce code s'avèrera utile à certains. N'hésitez pas à me faire part de vos remarques.

Codes Sources

A voir également

Ajouter un commentaire Commentaires
Messages postés
3
Date d'inscription
dimanche 3 novembre 2013
Statut
Membre
Dernière intervention
3 novembre 2013

Rien à dire, viens de verifier avec EXcel, les coefs sont justes, un est certes à 0 et n'apparait pas dans le cheminement du code, mais bien, simplement, génial, vais pouvoir travailler avec et approfondir ...excellent job, félicitations !!!
Messages postés
3
Date d'inscription
dimanche 3 novembre 2013
Statut
Membre
Dernière intervention
3 novembre 2013

Bonjour. Super cool comme soft. Quelqu'un aurait-il vérifié si les calculs sont fiables ? Cela serait génial et m'éviterait de le faire...sourire. Merci à tous
Messages postés
69
Date d'inscription
dimanche 16 mars 2003
Statut
Membre
Dernière intervention
14 mars 2015
2
Bonjour,

Ce code m'est utile, je l'utilise pour une régression quadratique, les coordonnes en x,y du pic me seront très utiles.
merci.
Messages postés
1
Date d'inscription
samedi 26 mars 2011
Statut
Membre
Dernière intervention
13 avril 2011

Bonjour, je suis actuellement en train de lancer ce code sous VB 2010 sans résultats. Si quelqu'un voudrait bien m'aider à l'adapter pour cette dernière version je lui en serait grandement reconnaissant.
Messages postés
132
Date d'inscription
mercredi 6 mars 2002
Statut
Membre
Dernière intervention
27 novembre 2012
1
Code Très intéressant,

il y a moyen de calculer l’écart Type ? et la quadratique ?
Afficher les 14 commentaires

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.