Le "fitting" de courbe (ou "ajustement" en francais) a pour objectif de trouver la meilleure courbe passant par un ensemble de points.
Par "meilleure courbe" on entend celle qui minimise une erreur. Dans notre cas, cette erreur est l'ecart quadratique entre la courbe et les points.
Lorsque la courbe recherchée est une forme linéaire, la methode des "moindres carrés" permet de calculer les coefficients de la courbe en résolvant un système linéraire.
Le code:
La classe Matrix, qui s'occupe de stocker et inverser une matrice:
Code java : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113 /** * Structure et Operations sur une matrice * * @author Xavier PHILIPPEAU * */ public class Matrix { // les données de la matrice ordonnnée en [x=colonne][y=ligne] public double value[][]; // les dimensions de la matrice public final int rows, cols; // constructeur de matrice quelconque (vide) public Matrix(int cols, int rows) { this.cols = cols; this.rows = rows; this.value = new double[this.cols][this.rows]; } // constructeur de matrice carré (pré-remplies) public Matrix(int size, double[][] squarearray) { this.cols = size; this.rows = size; this.value = squarearray; } /** * inversion par pivot de Gauss * * @return l'inverse de la matrice (ou null si non inversible) */ public Matrix inverse() { //Une matrice ne peut etre inversée que si elle est carrée. if (this.rows != this.cols) { System.err.println("Inversion : La matrice n'est pas carré !!"); return null; } // Création de la matrice de travail T = [ this | Identité ] Matrix T = new Matrix(this.cols * 2, this.cols); for (int y = 0; y < this.rows; y++) { for (int x = 0; x < this.cols; x++) { T.value[x][y] = this.value[x][y]; if (x == y) T.value[this.cols + x][y] = 1; } } // Pour chaque ligne de la matrice T for (int x = 0; x < T.rows; x++) { // Cherche la ligne avec le pivot max (en valeur absolue) int bestline = x; double pivot = T.value[x][bestline]; for (int y = x + 1; y < T.rows; y++) { if (Math.abs(T.value[x][y]) > Math.abs(pivot)) { bestline = y; pivot = T.value[x][bestline]; } } if (pivot == 0) { System.err.println("Inversion : Le pivot est nul,inversion impossible !!"); return null; } // Echange des lignes (si necessaire) if (bestline != x) { double tmp; for (int t = 0; t < T.cols; t++) { tmp = T.value[t][x]; T.value[t][x] = T.value[t][bestline]; T.value[t][bestline] = tmp; } } // Normalisation de la ligne du pivot for (int t = 0; t < T.cols; t++) T.value[t][x] /= pivot; // elimination des autres lignes for (int y = 0; y < T.rows; y++) { if (y==x) continue; double coef = T.value[x][y]; for (int t = 0; t < T.cols; t++) T.value[t][y] -= coef * T.value[t][x]; } } // recupere la partie droite de T qui contient l'inverse de la matrice // (la partie gauche devrait contenir l'identité) Matrix inverse = new Matrix(this.cols,this.rows); for (int y = 0; y < this.rows; y++) for (int x = 0; x < this.cols; x++) inverse.value[x][y] = T.value[this.cols+x][y]; return inverse; } /** * Multiplication de la matrice par un vecteur * * @param X le vecteur * @return le vecteur résultat: this[][] * X[] */ public double[] multiply(double[] X) { double[] Y = new double[this.rows]; for (int y = 0; y < this.rows; y++) for (int x = 0; x < this.cols; x++) Y[y] += this.value[x][y]*X[x]; return Y; } }
La classe LinearEquation qui s'occupe de modeliser une equation lineaire (coefs + valeur)
Code java : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44 /** * Structure d'une equation lineaire * * coef[0]*X0 + coef[1]*X1 + ... + coef[n]*Xn = residu * * ou X=(X0,X1,...,Xn) est le vecteur "inconnu" * * @author Xavier PHILIPPEAU */ public class LinearEquation { /** * Les coefficients de l'equation */ public double[] coefficents; /** * Le résidu (= le résultat) de l'equation */ public double residual; /** * Constructeur * * @param unknownsCount Le nombre d'inconnues */ public LinearEquation(int unknownsCount) { this.coefficents = new double[unknownsCount]; } /** * Constructeur * * @param coefficents Les coefficients de l'equation * @param residual Le résidu (= le résultat) de l'equation */ public LinearEquation(double[] coefficents, double residual) { this.coefficents = coefficents; this.residual = residual; } }
La classe LeastSquare qui s'occupe de construire et résoudre le système d'equations lineaires par la méthode des "moindres carrés":
Code java : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66 import java.util.List; /** * résolution du système d'equations par la methode des moindres carrés * * @author Xavier PHILIPPEAU */ public class LeastSquare { /** * Resolution d'un systeme lineaire * * @param systemsize Nombre d'inconnues dans le systeme * @param system Le systeme d'equations lineaires * @return le vecteur solution du système */ public static double[] linearFit(int systemsize, List<LinearEquation> system) { // taille du systeme a resoudre (= nombre d'inconnues) int n = systemsize; // On cherche a trouver X qui resoud le systeme: // // J.X = R // // on multiplie par Jt, la transposée de J: // // Jt.J.X = Jt.R // // on pose: // n // A[c][l] = Jt.J = Somme[ coefficient[c][k] * coefficient[l][k] ] // k=0 // // n // Y[l] = Jt.R = Somme[ residual[k] * coefficient[l][k] ] // k=0 // // et "n", le nombre d'equation dans le systeme. // // Le systeme a resoudre devient ainsi: // // A.X = Y Matrix A = new Matrix(n,n); double Y[] = new double[n]; for(LinearEquation equation:system) { for (int l = 0; l < n; l++) { for (int c = 0; c < n; c++) { A.value[c][l] += equation.coefficents[c] * equation.coefficents[l]; } Y[l] += equation.residual * equation.coefficents[l]; } } // Resolution par inversion: // // A . X = Y ==> X = A^-1 . Y double X[] = (A.inverse()).multiply(Y); return X; } }
PRomu@ld : sources
Partager