IdentifiantMot de passe
Loading...
Mot de passe oublié ?Je m'inscris ! (gratuit)
Navigation

Inscrivez-vous gratuitement
pour pouvoir participer, suivre les réponses en temps réel, voter pour les messages, poser vos propres questions et recevoir la newsletter

Contribuez Discussion :

[java] Moindres carrés pour systèmes lineaires [Sources]


Sujet :

Contribuez

  1. #1
    Rédacteur
    Avatar de pseudocode
    Homme Profil pro
    Architecte système
    Inscrit en
    Décembre 2006
    Messages
    10 062
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Âge : 51
    Localisation : France, Hérault (Languedoc Roussillon)

    Informations professionnelles :
    Activité : Architecte système
    Secteur : Industrie

    Informations forums :
    Inscription : Décembre 2006
    Messages : 10 062
    Points : 16 081
    Points
    16 081
    Par défaut [java] Fitting de courbe (moindres carrés)
    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
    ALGORITHME (n.m.): Méthode complexe de résolution d'un problème simple.

  2. #2
    Rédacteur
    Avatar de pseudocode
    Homme Profil pro
    Architecte système
    Inscrit en
    Décembre 2006
    Messages
    10 062
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Âge : 51
    Localisation : France, Hérault (Languedoc Roussillon)

    Informations professionnelles :
    Activité : Architecte système
    Secteur : Industrie

    Informations forums :
    Inscription : Décembre 2006
    Messages : 10 062
    Points : 16 081
    Points
    16 081
    Par défaut Exemples d'utilisation
    Trouver la "meilleure ligne" au sens 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
     
    public static double[] bestLine(List<double[]> data) {
     
    	List<LinearEquation> system = new ArrayList<LinearEquation>();
     
    	for(double[] d : data) {
    		double x = d[0];
    		double y = d[1];
     
    		// a.1 + b.X = Y
    		//
    		// (X,Y) sont connus (données du probleme)
    		// (a,b) sont les 2 inconnues
     
    		LinearEquation equation = new LinearEquation(2);
    		equation.coefficents = new double[] { 1, x };
    		equation.residual = y;
     
    		system.add(equation);
    	}
     
    	// resolution du systeme lineaire
    	double[] C = LeastSquare.linearFit(2, system);
    	return C;
     
    	// La solution s'ecrit: Y = C[0] + C[1]*X
    }

    Trouver la "meilleure parabole" au sens 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
     
    public static double[] BestParabol(List<double[]> data) {
     
    	List<LinearEquation> system = new ArrayList<LinearEquation>();
     
    	for(double[] d : data) {
    		double x = d[0];
    		double y = d[1];
     
    		// a.1 + b.X + c.X^2 = Y
    		//
    		// (X,Y) sont connus (données du probleme)
    		// (a,b,c) sont les 3 inconnues
     
    		LinearEquation equation = new LinearEquation(3);
    		equation.coefficents = new double[] { 1, x, x*x };
    		equation.residual = y;
     
    		system.add(equation);
    	}
     
    	// resolution du systeme lineaire
    	double[] C = LeastSquare.linearFit(3, system);
    	return C;
     
    	// La solution s'ecrit: Y = C[0] + C[1]*X + C[2]*X^2
    }

    Trouver le "meilleure polynome de degré n" au sens 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
     
    public static double[] BestPolynom(int order, List<double[]> data) {
     
    	List<LinearEquation> system = new ArrayList<LinearEquation>();
     
    	for(double[] d : data) {
    		double x = d[0];
    		double y = d[1];
     
    		// a.1 + b.X + c.X^2 + ... + z.X^ordre = Y
    		//
    		// (X,Y) sont connus (données du probleme)
    		// (a,b,c,...z) sont les inconnues
     
    		LinearEquation equation = new LinearEquation(order+1);
     
    		equation.residual = y;
    		for(int n=0;n<=order;n++)
    			equation.coefficents[n] = Math.pow(x, n);
     
    		system.add(equation);
    	}
     
    	// resolution du systeme lineaire
    	double[] C = LeastSquare.linearFit(order+1, system);
    	return C;
     
    	// La solution s'ecrit: Y = C[0] + C[1]*X + C[2]*X^2 + ... + C[n]*X^n
    }

    un exemple/explication du cas "non lineaire" appliqué à un cercle:
    http://www.orbitals.com/self/least/least.htm
    ALGORITHME (n.m.): Méthode complexe de résolution d'un problème simple.

  3. #3
    Expert éminent sénior

    Profil pro
    Inscrit en
    Janvier 2007
    Messages
    10 603
    Détails du profil
    Informations personnelles :
    Âge : 66
    Localisation : France

    Informations forums :
    Inscription : Janvier 2007
    Messages : 10 603
    Points : 17 913
    Points
    17 913
    Billets dans le blog
    2
    Par défaut
    je la traduirai en C (Java moi yen a pas connaître) et je posterais le texte du papier cette semaine (peut-être aujourdhui ?) .. (je viens de rentrer de voyage..)
    "Un homme sage ne croit que la moitié de ce qu’il lit. Plus sage encore, il sait laquelle".

    Consultant indépendant.
    Architecture systèmes complexes. Programmation grosses applications critiques. Ergonomie.
    C, Fortran, XWindow/Motif, Java

    Je ne réponds pas aux MP techniques

  4. #4
    Rédacteur
    Avatar de pseudocode
    Homme Profil pro
    Architecte système
    Inscrit en
    Décembre 2006
    Messages
    10 062
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Âge : 51
    Localisation : France, Hérault (Languedoc Roussillon)

    Informations professionnelles :
    Activité : Architecte système
    Secteur : Industrie

    Informations forums :
    Inscription : Décembre 2006
    Messages : 10 062
    Points : 16 081
    Points
    16 081
    Par défaut
    Trouver la "meilleure ellipse" au sens des moindres carrés, en utilisant les formules de mathworld.wolfram.com:

    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
     
    public static double[] BestEllipsoid(List<double[]> data) {
     
    	List<LinearEquation> system = new ArrayList<LinearEquation>();
     
    	for(double[] d : data) {
    		double x = d[0];
    		double y = d[1];
     
    		// a.x^2 + 2b.xy + c.y^2 + 2d.x + 2f.y + g = 0
    		//
    		// (X,Y) sont connus (données du probleme)
    		// on pose g = -1
    		// (a,b,c,d,f) sont les 5 inconnues
     
    		LinearEquation equation = new LinearEquation(5);
    		equation.coefficents = new double[] { x*x, x*y, y*y, x, y };
    		equation.residual = 1.0;
     
    		system.add(equation);
    	}
     
    	// resolution du systeme lineaire
    	double[] C = LeastSquare.linearFit(5, system);
     
    	// evaluation des 5 inconnues
    	double a = C[0];
    	double b = C[1]/2;
    	double c = C[2];
    	double d = C[3]/2;
    	double f = C[4]/2;
    	double g = -1;
     
    	// calculs des caracteristiques de l'ellipse:
    	//
    	// d'après http://mathworld.wolfram.com/Ellipse.html
    	//
    	// formules (15), (16), ..., (23)
     
    	// verifie que l'ellipse est valide
    	double delta=a*(c*g-f*f)-b*(b*g-d*f)+d*(b*f-d*c);
    	double J=a*c-b*b;
    	double I=a+c;
    	if(delta==0) return null;
    	if(J<=0) return null;
    	if((delta/I)>=0) return null;
     
    	/* calcul du centre de l'ellipse */
    	double xcentre = (c*d-b*f)/(b*b-a*c);
    	double ycentre = (a*f-b*d)/(b*b-a*c);
     
    	/* angle de rotation (en radian) */
    	double angle = 0.5*Math.atan(2*b/(c-a));
     
    	/* longueurs des 2 demi-axes */
    	double num = 2*(a*f*f + c*d*d + g*b*b - 2*b*d*f - a*c*g);
    	double sq = Math.sqrt( 1 + 4*b*b/((a-c)*(a-c)) );
     
    	double denom1 = (b*b-a*c) * ( (c-a) * sq - (c+a) );
    	double xhalflength = Math.sqrt(num/denom1);
     
    	double denom2 = (b*b-a*c) * ( (a-c) * sq - (c+a) );
    	double yhalflength = Math.sqrt(num/denom2);
     
    	/* retourne toujours un angle positif  */
    	if (angle<0) {
    		angle+=Math.PI/2;
    		// echange les longueurs des demi-axes
    		double tmp=xhalflength;
    		xhalflength=yhalflength;
    		yhalflength=tmp;
    	}
     
    	return new double[] { xcentre, ycentre, angle, xhalflength, yhalflength};
    }

    Et un petit test pour s'assurer que ca marche...
    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
     
    // construction des points d'une ellipse
    List<double[]> data = new ArrayList<double[]>();
    for(int i=0;i<100;i++) {
    	double t = i/100.0;
    	double a = t*2*Math.PI;
     
    	// ellipse de demi-axe 17 x 10, centrée en (0,0)
    	double xe = 17*Math.cos(a);
    	double ye = 10*Math.sin(a);
     
    	// + une rotation de 45 degrés
    	double r = Math.toRadians(45);
    	double xr = xe*Math.cos(r)+ye*Math.sin(r);
    	double yr = -xe*Math.sin(r)+ye*Math.cos(r);
     
    	// + une translation de (15,8)
    	double x = xr+15;
    	double y = yr+8;
     
    	data.add( new double[] {x,y});
    }
     
    double[] E = BestEllipsoid(data);
     
    System.out.println("centre: "+E[0]+","+E[1]);
    System.out.println("angle: "+Math.toDegrees(E[2])+" (degrés)");
    System.out.println("demi-axe: "+E[3]+" x "+E[4]);
     
    // affiche:
    // centre: 14.999999999999954,8.000000000000052
    // angle: 44.99999999999916 (degrés)
    // demi-axe: 17.00000000000003 x 10.000000000000034
    ALGORITHME (n.m.): Méthode complexe de résolution d'un problème simple.

  5. #5
    Rédacteur

    Homme Profil pro
    Comme retraité, des masses
    Inscrit en
    Avril 2007
    Messages
    2 978
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Âge : 83
    Localisation : Suisse

    Informations professionnelles :
    Activité : Comme retraité, des masses
    Secteur : Industrie

    Informations forums :
    Inscription : Avril 2007
    Messages : 2 978
    Points : 5 179
    Points
    5 179
    Par défaut [java] Fitting de courbe (moindres carrés)
    Deux remarques à l'intention de pseudocode:

    1) Si la tâche à accomplir consiste à résoudre un système linéaire, ce qui est très souvent le cas, il est inutilement coûteux en temps d'exécution d'inverser la matrice (bien que ce soit ce que font la plupart des gens); le plus efficace est de la factoriser par la méthode LU dans le cas général, ou par la méthode de Cholesky, deux fois plus rapide, si elle est symétrique définie positive.

    2) Si tu cherches le vecteur qui satisfait le moins mal possible (au sens des moindres carrés) un système linéaire surdéterminé, regarde la méthode des valeurs singulières (singular value decomposition); tu trouveras toute la documentation à ce sujet dans Numerical Recipes ou dans la description de la fonction svd de MatLab.

    A ta disposition si tu as encore des questions à ce sujet.
    Calcul numérique de processus industriels
    Formation, conseil, développement

    Point n'est besoin d'espérer pour entreprendre, ni de réussir pour persévérer. (Guillaume le Taiseux)

  6. #6
    Candidat au Club
    Inscrit en
    Décembre 2006
    Messages
    3
    Détails du profil
    Informations forums :
    Inscription : Décembre 2006
    Messages : 3
    Points : 4
    Points
    4
    Par défaut
    bonjour,
    je souhaite demander si je peux utiliser cette méthode pour résoudre un système de la forme : X^t*A=o , avec A matrice non carrée(n*m)
    merci d'avance

  7. #7
    Rédacteur
    Avatar de pseudocode
    Homme Profil pro
    Architecte système
    Inscrit en
    Décembre 2006
    Messages
    10 062
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Âge : 51
    Localisation : France, Hérault (Languedoc Roussillon)

    Informations professionnelles :
    Activité : Architecte système
    Secteur : Industrie

    Informations forums :
    Inscription : Décembre 2006
    Messages : 10 062
    Points : 16 081
    Points
    16 081
    Par défaut
    Citation Envoyé par souheyeb Voir le message
    bonjour,
    je veus demander si je peut utiliser cette methode pour resoudre un systeme de la forme : X^t*A=o , avec A matrice non carrée(n*m)
    merci d'avance
    Heu... Non. La méthode des moindres carrés a pour objectif de trouver les "meilleurs" coefficients d'une courbe. La résolution du système linéaire surdéterminé est l'une des étapes de cette méthode, mais nullement le but en soi. D'ailleurs, comme l'a fait remarquer FR119492, elle n'est meme pas obligatoire.

    Pour la résolution d'un système linéaire surdéterminé, je te conseille plutôt la factorisation QR ou la décomposition SVD.
    ALGORITHME (n.m.): Méthode complexe de résolution d'un problème simple.

  8. #8
    Nouveau Candidat au Club
    Inscrit en
    Mai 2012
    Messages
    1
    Détails du profil
    Informations forums :
    Inscription : Mai 2012
    Messages : 1
    Points : 1
    Points
    1
    Par défaut
    Pour passer en 3d et chercher l'equation d'un plan , il suffit de déclarer
    double z = d[2];

    et
    de chercher l'equation lineaire z=aX+bY+c ?

    UN petit coup de main pour un informaticien non mathematicien ? (pitié!)

  9. #9
    Invité
    Invité(e)
    Par défaut Une librairie ?
    Bonjour,

    Merci pour ce sujet, pseudocode.
    On peut aussi utiliser une librairie d'algèbre linéaire pour créer le système linéaire et le résoudre... (ça fait moins de code à écrire et il y a peut-être des améliorations pour la résolution du système quand il est mal conditionné, par rapport à tA.A x = tA.b )
    Jama, par exemple ?
    http://math.nist.gov/javanumerics/jama/

    Exemple :
    Code : Sélectionner tout - Visualiser dans une fenêtre à part
    1
    2
    3
    4
    5
    6
    7
     
          double[][] array = {{1.,2.,3},{4.,5.,6.},{7.,8.,10.}};
          Matrix A = new Matrix(array);
          Matrix b = Matrix.random(3,1);
          Matrix x = A.solve(b);
          Matrix Residual = A.times(x).minus(b);
          double rnorm = Residual.normInf();

  10. #10
    Invité
    Invité(e)
    Par défaut Version avec Jama
    Citation Envoyé par FR119492 Voir le message
    1) Si la tâche à accomplir consiste à résoudre un système linéaire, ce qui est très souvent le cas, il est inutilement coûteux en temps d'exécution d'inverser la matrice (bien que ce soit ce que font la plupart des gens); le plus efficace est de la factoriser par la méthode LU dans le cas général, ou par la méthode de Cholesky, deux fois plus rapide, si elle est symétrique définie positive.
    A priori, pour une résolution "au mieux" (matrice rectangulaire, nb lignes > nb colonnes), Jama utilise une décomposition QR...
    Code Java : Sélectionner tout - Visualiser dans une fenêtre à part
    1
    2
    3
    4
    5
     
       public Matrix solve (Matrix B) {
          return (m == n ? (new LUDecomposition(this)).solve(B) :
                           (new QRDecomposition(this)).solve(B));
       }

    J'ai écrit une version dérivée du code de Xavier, adaptée à l'utilisation de Jama, voilà ce que ça donne :

    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
     
    import Jama.Matrix;
    import java.awt.geom.Point2D;
    import java.util.AbstractList;
     
    [...]
     
    	/**
             * Renvoie les paramètre de l'ellipse qui s'ajuste mieux au nuage de
             *  points 2D reçu (au sens des moindres carrés)
             * @param lstPts2D liste de points 2D (Point2D.Double) du contour de l'ellipse
             * @return un tableau qui contient les paramètres de l'ellipse
             */
    	public static double[] getParametres_ellipse_mc( AbstractList lstPts2D ) {
    		// construction du système linéaire
    		int nNb_pts = lstPts2D.size();
    		Matrix mA = new Matrix( nNb_pts, 5 );
    		Matrix vB = new Matrix( nNb_pts, 1 );
    		for( int i = 0; i < nNb_pts; ++i ) {
    			Point2D.Double pt = (Point2D.Double)lstPts2D.get( i );
    			mA.set( i, 0, pt.x * pt.x );
    			mA.set( i, 1, pt.x * pt.y );
    			mA.set( i, 2, pt.y * pt.y );
    			mA.set( i, 3, pt.x );
    			mA.set( i, 4, pt.y );
    			vB.set(i, 0, -1.0 );
    		}
    		// resolution du systeme lineaire A.x = b
    		Matrix vX = mA.solve( vB );
     
    		// evaluation des 5 inconnues
    		double a = vX.get( 0, 0 );
    		double b = vX.get( 1, 0 ) / 2.0;
    		double c = vX.get( 2, 0 );
    		double d = vX.get( 3, 0 ) / 2.0;
    		double f = vX.get( 4, 0 ) / 2.0;
    		double g = 1.0;
     
    		// calcule les caracteristiques de l'ellipse:
    		//
    		// d'après http://mathworld.wolfram.com/Ellipse.html
    		//
    		// formules (15), (16), ..., (23)
     
    		// verifie que l'ellipse est valide
    		double delta = a * ( c * g - f * f) - b * ( b * g - d * f ) + d * ( b * f - d * c );
    		double nJ = a * c - b * b;
    		double nI = a + c;
    		if ( delta == 0 ) return null;
    		if ( nJ <= 0 ) return null;
    		if ( ( delta / nI ) >= 0) return null;
     
    		// calcul du centre de l'ellipse
    		double nX_centre = ( c*d - b*f )/( b*b - a*c );
    		double nY_centre = ( a*f - b*d )/( b*b - a*c );
     
    		// angle de rotation (en radian)
    		double nAngle_rad = 0.5 * Math.atan( 2.0 * b / ( c-a ));
     
    		// longueurs des 2 demi-axes
    		double num = 2.0 * ( a * f*f + c * d*d + g* b*b - 2*b*d*f - a*c*g);
    		double sq = Math.sqrt( (a-c)*(a-c) + 4.0*b*b );
     
    		double denom1 = ( b*b - a*c ) * ( sq - ( a + c ) );
    		double nDemi_axe_x = Math.sqrt( num / denom1 );
     
    		double denom2 = ( b*b - a*c ) * ( -1.0 * sq - ( a + c ) );
    		double nDemi_axe_y = Math.sqrt( num / denom2 );
     
    		// retourne un angle entre [0 et 2\Pi[
    		if ( nAngle_rad < 0.0 ) {
    			nAngle_rad += 2.0 * Math.PI;
    		}
    		return new double[]{nX_centre, nY_centre, nAngle_rad, nDemi_axe_x, nDemi_axe_y };
    	}// getParametres_ellipse_mc

    Avec la même procédure de test que Xavier j'obtiens aussi les bons résultats :
    Code : Sélectionner tout - Visualiser dans une fenêtre à part
    1
    2
    3
    4
     
    centre: 14.999999999999991,7.999999999999996
    angle: 44.99999999999999 (degrés)
    demi-axe: 17.000000000000025 x 9.999999999999982
    Merci encore pour ce sujet, il m'a bien aidé...

  11. #11
    Invité
    Invité(e)
    Par défaut
    La version "EllipseFitting" de Michael Doube est intéressante (voir le code ici) :
    http://github.com/mdoube/BoneJ/blob/...itEllipse.java

    Voir notamment la fonction "varToDimensions"...

  12. #12
    Expert éminent sénior

    Profil pro
    Inscrit en
    Janvier 2007
    Messages
    10 603
    Détails du profil
    Informations personnelles :
    Âge : 66
    Localisation : France

    Informations forums :
    Inscription : Janvier 2007
    Messages : 10 603
    Points : 17 913
    Points
    17 913
    Billets dans le blog
    2
    Par défaut
    tu peux regarder j'en ai mise une ici-même en Fortran, facilement convertible..
    "Un homme sage ne croit que la moitié de ce qu’il lit. Plus sage encore, il sait laquelle".

    Consultant indépendant.
    Architecture systèmes complexes. Programmation grosses applications critiques. Ergonomie.
    C, Fortran, XWindow/Motif, Java

    Je ne réponds pas aux MP techniques

  13. #13
    Membre expérimenté
    Homme Profil pro
    Chercheur
    Inscrit en
    Mars 2010
    Messages
    1 218
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Paris (Île de France)

    Informations professionnelles :
    Activité : Chercheur

    Informations forums :
    Inscription : Mars 2010
    Messages : 1 218
    Points : 1 685
    Points
    1 685
    Par défaut
    Bonjour,

    inverser une matrice pour ensuite faire
    // Resolution par inversion:
    //
    // A . X = Y ==> X = A^-1 . Y

    double X[] = (A.inverse()).multiply(Y);
    c'est vraiment très maladroit!

  14. #14
    Membre confirmé

    Profil pro
    Inscrit en
    Janvier 2008
    Messages
    786
    Détails du profil
    Informations personnelles :
    Localisation : France

    Informations forums :
    Inscription : Janvier 2008
    Messages : 786
    Points : 602
    Points
    602
    Par défaut
    C'est exactement ce que je cherchais.

    (enfin je cherchais un plane fitting)

    Merci.

  15. #15
    Modérateur
    Avatar de ToTo13
    Homme Profil pro
    Chercheur en informatique
    Inscrit en
    Janvier 2006
    Messages
    5 793
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Âge : 45
    Localisation : Etats-Unis

    Informations professionnelles :
    Activité : Chercheur en informatique
    Secteur : Santé

    Informations forums :
    Inscription : Janvier 2006
    Messages : 5 793
    Points : 9 860
    Points
    9 860
    Par défaut
    En revanche, je pense que cette méthode doit pouvoir être utilisée pour trouver la meilleure gaussienne 2D passant par un ensemble de points.
    Ou ai je tort ?

    [EDIT] : J'ai définitivement tort car l'équation de la Gaussienne n'est pas linéaire à cause de l'exponentielle.
    Donc comment modifier la méthode pour résoudre dans le cas Gaussien ?
    Consignes aux jeunes padawans : une image vaut 1000 mots !
    - Dans ton message respecter tu dois : les règles de rédaction et du forum, prévisualiser, relire et corriger TOUTES les FAUTES (frappes, sms, d'aurteaugrafe, mettre les ACCENTS et les BALISES) => ECRIRE clairement et en Français tu DOIS.
    - Le côté obscur je sens dans le MP => Tous tes MPs je détruirai et la réponse tu n'auras si en privé tu veux que je t'enseigne.(Lis donc ceci)
    - ton poste tu dois marquer quand la bonne réponse tu as obtenu.

Discussions similaires

  1. Réponses: 0
    Dernier message: 20/08/2013, 15h29
  2. Réponses: 2
    Dernier message: 27/01/2009, 12h05
  3. moindre carre pour approcher une courbe en S
    Par hamska2 dans le forum Mathématiques
    Réponses: 5
    Dernier message: 19/05/2008, 10h39
  4. Moindres carres pour une droite en trois dimensions
    Par shindara dans le forum Traitement d'images
    Réponses: 4
    Dernier message: 28/06/2007, 00h08
  5. Moindres carrés pour courbe
    Par cjacquel dans le forum Mathématiques
    Réponses: 3
    Dernier message: 31/03/2007, 19h02

Partager

Partager
  • Envoyer la discussion sur Viadeo
  • Envoyer la discussion sur Twitter
  • Envoyer la discussion sur Google
  • Envoyer la discussion sur Facebook
  • Envoyer la discussion sur Digg
  • Envoyer la discussion sur Delicious
  • Envoyer la discussion sur MySpace
  • Envoyer la discussion sur Yahoo