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

Télécharger le projet
Commenter

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.
Télécharger le projet

Codes Sources

A voir également

S'inscrire maintenant

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.

Du même auteur (caco64)