RÉGRESSION NON-LINÉAIRE POLYNOMIALE PAR LA MÉTHODE DES MOINDRES CARRÉS

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