Moindres carrés surface tendance ordre 2

Version imprimable

Bonjour à tous,

j'ai un nuage de points que je veux modéliser par une surface. J'ai développé un algo qui fait une régression d'ordre 1(plan) et qui fonctionne. J'ai ensuite trouvé sur internet les formules pour la modélisation d'une surface d'ordre 2

http://support.articque.com/samples/doc/surf_t2.htm

Voici l'algo mais il ne marche pas

où est la faute?

Code:

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
114
115
116
117
118
119
120
121
function [a,b,c,d,e,f,g,h,i] = reg_ordre2(N,X,Y,Z)
 
Xi = 0;
Yi = 0;
Zi = 0;
Xi2 = 0;
Yi2 = 0;
Xi3 = 0;
Yi3 = 0;
Xi4 = 0;
Yi4 = 0;
 
XiYi = 0;
XiYi2 = 0;
XiYi3 = 0;
XiYi4 = 0;
 
 
Xi2Yi = 0;
Xi2Yi2 = 0;
Xi2Yi3 = 0;
Xi2Yi4 = 0;
 
Xi3Yi = 0;
Xi3Yi2 = 0;
Xi3Yi3 = 0;
Xi3Yi4 = 0;
 
Xi4Yi = 0;
Xi4Yi2 = 0;
Xi4Yi3 = 0;
Xi4Yi4 = 0;
 
XiZi = 0;
YiZi = 0;
XiYiZi = 0;
Xi2Zi = 0;
Yi2Zi = 0;
Xi2YiZi = 0;
XiYi2Zi = 0;
Xi2Yi2Zi = 0;
 
 
for i = 1:N
    Xi = Xi + N*X(i);
    Yi = Yi + N*Y(i);
    Xi2 = Xi2 + N*X(i)^2;
    Yi2 = Yi2 + N*Y(i)^2;
    Xi3 = Xi3 + N*X(i)^3;
    Yi3 = Yi3 + N*Y(i)^3;
    Xi4 = Xi4 + N*X(i)^4;
    Yi4 = Yi4 + N*Y(i)^4;
end
 
for i = 1:N
    for t = 1:N
        XiYi = XiYi + X(i)*Y(t);
        XiYi2 = XiYi2 + X(i)*Y(t)^2;
        XiYi3 = XiYi3 + X(i)*Y(t)^3;
        XiYi4 = XiYi4 + X(i)*Y(t)^4;
 
        Xi2Yi = Xi2Yi + X(i)^2*Y(t);
        Xi2Yi2 = Xi2Yi2 + X(i)^2*Y(t)^2;
        Xi2Yi3 = Xi2Yi3 + X(i)^2*Y(t)^3;
        Xi2Yi4 = Xi2Yi4 + X(i)^2*Y(t)^4;
 
        Xi3Yi = Xi3Yi + X(i)^3*Y(t);
        Xi3Yi2 = Xi3Yi2 + X(i)^3*Y(t)^2;
        Xi3Yi3 = Xi3Yi3 + X(i)^3*Y(t)^3;
        Xi3Yi4 = Xi3Yi4 + X(i)^3*Y(t)^4;
 
        Xi4Yi = Xi4Yi + X(i)^4*Y(t);
        Xi4Yi2 = Xi4Yi2 + X(i)^4*Y(t)^2;
        Xi4Yi3 = Xi4Yi3 + X(i)^4*Y(t)^3;
        Xi4Yi4 = Xi4Yi4 + X(i)^4*Y(t)^4;
 
        XiZi = XiZi + X(i)*Z(i,t);
        YiZi = YiZi + Y(t)*Z(i,t);
        XiYiZi = XiYiZi + Y(t)*X(i)*Z(i,t);
        Xi2YiZi = Xi2YiZi + Y(t)*X(i)^2*Z(i,t);
        XiYi2Zi = XiYi2Zi + Y(t)^2*X(i)*Z(i,t);
        Xi2Yi2Zi = Xi2Yi2Zi + Y(t)^2*X(i)^2*Z(i,t);
        Xi2Zi = Xi2Zi + X(i)^2*Z(i,t);
        Yi2Zi = Yi2Zi + Y(t)^2*Z(i,t);
        Zi = Zi + Z(i,t);
    end
end
 
P = [Xi2    XiYi   Xi2Yi  Xi3    XiYi2  Xi3Yi  Xi2Yi2 Xi3Yi2 Xi;
     XiYi   Yi2    XiYi2  Xi2Yi  Yi3    Xi2Yi2 XiYi3  Xi2Yi3 Yi;
     XiYi   XiYi2  Xi2Yi2 Xi3Yi  XiYi3  Xi3Yi2 Xi2Yi3 Xi3Yi3 XiYi
     Xi3    Xi2Yi  Xi3Yi  Xi4    Xi2Yi2 Xi4Yi  Xi3Yi2 Xi4Yi2 Xi2;
     XiYi2  Yi3    XiYi3  Xi2Yi2 Yi4    Xi2Yi3 XiYi4  Xi2Yi4 Yi2;
     Xi3Yi  Xi2Yi2 Xi3Yi2 Xi4    Xi2Yi3 Xi4Yi2 Xi3Yi3 Xi4Yi3 Xi2Yi;
     Xi2Yi2 XiYi3  Xi2Yi3 Xi3Yi2 XiYi4  Xi3Yi3 Xi2Yi4 Xi3Yi4 XiYi2;
     Xi3Yi2 Xi2Yi3 Xi3Yi3 Xi4Yi2 Xi2Yi4 Xi4Yi3 Xi3Yi4 Xi4Yi4 Xi2Yi2;
     Xi     Yi     XiYi   Xi2    Yi2    Xi2Yi  XiYi2  Xi2Yi2 N^2];
 
V = [XiZi;
     YiZi;
     XiYiZi;
     Xi2Zi;
     Yi2Zi;
     Xi2YiZi;
     XiYi2Zi;
     Xi2Yi2Zi;
     Zi];
 
R = inv(P)*V;
 
a = R(1);
b = R(2);
c = R(3);
d = R(4);
e = R(5);
f = R(6);
g = R(7);
h = R(8);
i = R(9);
 
return

18/04/2007, 17h29
Jerome Briot

Citation:

Envoyé par tony57

Voici l'algo mais il ne marche pas

Pour faire court... c'est quoi qui marche pas ? Matlab retourne une erreur ou bien le résultat retourné n'est pas le bon ?

Note: pour la résolution du système, il est préférable d'utiliser l'opérateur backslash

il y a plus simple :

ta surface est en X^i*Y^j, i et j allant chacun de 0 à 2 sans que leur somme dépasse 2.

Tu construit la matrice M
Code:

1 2 3 4 5 6 M = [] ; X=X(:) ; Y=Y(:) ; % au besoin for i=0:2, for j=0:2, if (i+j)>2, continue, end M = [M (X.^i) .* (Y.^j)] ; end, end
Maintenant, tu veux résoudre aux moindres carrés l'équation en k :
M * k = Z
en Matlab, la solution est directe grâce à l'anti-division

Code:

k = M \ Z ;

maintenant tu n'as plus qu'à te servir de tes coefficients k pour calculer une nouvelle grille
Code:

1 2 3 4 5 6 7 8 [mesX, mesY] = meshgrid(x0:dx:x1, y0:dy:y1) ; s = size(mesX) ; % pour le remettre comme il était après mesX=mesX(:) ; mesY=mesY(:) ; monM = [] ; for i=0:2, for j=0:2, if (i+j)>2, continue, end monM = [monM (mesX.^i) .* (mesY.^j)] ; end, end mesZ = monM * k ;
et profiter du résultat

Code:

surf (reshape(mesX,s), reshape(mesY,s), reshape(mesZ,s)) ;

En résumé :
toi, tu as construit la matrice M' * M et tu l'as inversée.
c'est correct (à un bug près ;) ) mais :
1/ il y a une écriture plus simple de ce que tu as fait : k = inv (M' * M) * M' * Z et,
2/ Matlab fait mieux avec l'antidivision : k = M\Z ;

Edit : autres commentaires transférés dans le post suivant

OL

Bjr,

Ton problème est le deuxième en deux jours qui se résouds facilement en connaissant l'anti-division de Matlab (le back slash \).

C'est pas très étonnant parceque dierrière l'opérateur \, il y a vraiment du monde, et toute une théorie qui s'appelle la modélisation linéaire. celui qui n'a pas de notions à-dessus e peut pas deviner que son problème se résoud en une seule ligne : k = X\Z.

Je t'ai fait un petit tuto d'introduction à la modélisation linéaire. Il présente d'une manière un peu plus complète et générale ce que je t'ai déja expliqué au post précédent. J'espère qu'il te sera utile. L'ensemble du programme est recopié en bloc dans la dernière boite.

% 1/ préparation des données à simuler
Code:

1 2 3 4 5 6 7 8 9 10 % 1.a) X et Y sont tirés dans [0 1[ taille=600; X=rand(taille,1) ; Y=rand(taille,1) ; % 1.b) on simule Z = X^a * Y^b après avoir ajouté du bruit à X et Y bruit = 0.1 ; a=1.75 ; b=-0.9 ; Z=(X+bruit*rand(taille,1)).^a .* (Y+bruit*rand(taille,1)).^b ;
% 2/ définition du modèle à résoudre
PS. Tu verra que j'ai viré les dernières boucles for du code. En Matlab, quand on a un programme avec des for ... end imbriqués, c'est une bonne idée d'aller étudier l'aide ou poser une petite question sur un forum.
Code:

1 2 3 4 5 6 7 px=2 ; % degré max en X py=2 ; % degré max en Y ptot=2 ; % degré max total [i,j]=meshgrid(0:px,0:py) ; % génération de toutes les possibilités u=find(i+j<=ptot) ; % écarter les degrés trop élevés i=i(u) ; j=j(u) ; nk=length(u) ; % nombre de paramètres du modèle
% 3/ Résolution du problème
J'ai mis deux méthodes. la méthode classique de cours, et la méthode Matlab. Commençons par la classique :
Le problème posé s'écrit sous forme matricielle M*k = Z. (voir construction de M dans le bloc de code suivant).
L'esprit général est de multiplier à gauche les deux membres de l'égalité par inv(M), ce qui donnerait la solution k = inv(M) * Z
Mais il y a un os : M n'est pas carrée ! donc clairement pas inversible.
on contourne ce problème en multipliant M par sa transposée M'. Le produit (M' * M) est non seulemnt une matrice carrée, mais c'est en plus une matrice dite "symétrique définie positive". (symétrique : elle est égale à sa transposée. définie positives : toutes ses valeurs propres sont >0. Ces matrices ont des propriétés remaquables, en particulier pour leur décomposition.)

M' * M * k = M' * Z

En multipliant, encore à gauche, par l'inverse de (M' * M) :

inv(M' * M) * (M' * M) * k = inv (M'*M) * M' * Z

ce qui donne la solution puisque tout le tremblement dans le membre de gauche s'annule :

k = inv (M' * M) * M' * Z

Ouf !!!

Et Matlab ?

Matlab te fait tout ça en un clin d'oeuil et plus proprement. Sans boucles for imbriquées, et, surtout, sans inversion de matrice (ni demandée, ni cachée), grâce à ces fameuses propriétés de décomposition des matrices symétriques définies positives qui sont pleinement exploitées sous le capot du puissant opérateur \.

solution Matlab : k = M\Z ;
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 % 3.a) on calcule 4 matrices (taille x nk) nx = length(X) ; % nombre de données ? X=repmat(X,1,nk) ; Y=repmat(Y,1,nk) ; i=repmat(i',nx,1) ; j=repmat(j',nx,1) ; % 3.b) on les assemble dans la matrice M M=X.^i.*Y.^j ; % 3.c1) résolution régulière en Matlab k_matlab = M\Z ; % 3.c2) la méthode que tu as employée (dans une version plus compacte) % C'est celle qu'on apprend en cours : le truc qui faut comprendre mais que personne n'utilise. k_manuel = inv(M' * M) * M' * Z ; % 3.d) comparaison des coefficients du polynome par les deux méthodes % C'est quasiment la meme chose, mais pas tout à fait. % Dans certains cas dégénérés, l'inversion peut etre fausse et l'anti-divison contnuer à trouver le bon résultat. disp('Réultats de chaque calcul') ; disp ([k_matlab k_manuel]) ; disp('différence relative entre les deux méthodes') ; disp(k_matlab-k_manuel) ;
% 4. Calcul des Z simulés
Code:

1 2 3 4 5 6 Z_matlab = M*k_matlab ; Z_manuel = M*k_manuel ; differences = (Z_matlab-Z_manuel) ; disp ('l''ecart type des differences en Z entre chaque version est') ; disp (std(differences)) ;
% 5. calcul du coeficient de correlation
Code:

1 2 3 residus = Z-Z_matlab ; r2 = (var(Z)-var(residus))/var(Z) ; disp(sprintf('coefficient de correlation r2 = %f', r2)) ;
Finalement tout le code pour copier/coller et jouer avec
Code:

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 % 1/ préparation des données à simuler % 1.a) X et Y sont tirés dans [0 1[ taille=6000; X=rand(taille,1) ; Y=rand(taille,1) ; % 1.b) on simule Z = X^a * Y^b après avoir ajouté du bruit à X et Y bruit = 0.1 ; a=1.75 ; b=-0.9 ; Z=(X+bruit*rand(taille,1)).^a .* (Y+bruit*rand(taille,1)).^b ; % 2/ définition du modèle à résoudre px=13 ; % degré max en X py=13 ; % degré max en Y ptot=13 ; % degré max total [i,j]=meshgrid(0:px,0:py) ; % génération de toutes les possibilités u=find(i+j<=ptot) ; % écarter les degrés trop élevés i=i(u) ; j=j(u) ; nk=length(u) ; % finalement, combien de paramètres à calculer ? % 3/ Résolution du problème % 3.a) on calcule 4 matrices (taille x nk) nx = length(X) ; % combien de données ? X=repmat(X,1,nk) ; Y=repmat(Y,1,nk) ; i=repmat(i',nx,1) ; j=repmat(j',nx,1) ; % 3.b) on les assemble dans la matrice M M=X.^i.*Y.^j ; % 3.c1) résolution normale par la formule que je t'ai indiquée k_matlab = M\Z ; % 3.c2) la méthode que tu as employée (dans une version plus compacte) k_manuel = inv(M' * M) * M' * Z ; % 3.d) comparaison des coefficients du polynome par les deux méthodes % C'est quasiment la meme chose, mais pas tout à fait. % Dans certains cas dégénérés, l'inversion peut etre fausse et l'anti-divison contnuer à trouver le bon résultat. disp('Réultats de chaque calcul') ; disp ([k_matlab k_manuel]) ; disp('différence relative entre les deux méthodes') ; disp(k_matlab-k_manuel) ; % 4. Calcul des Z simulés Z_matlab = M*k_matlab ; Z_manuel = M*k_manuel ; differences = (Z_matlab-Z_manuel) ; disp ('l''ecart type des differences en Z entre chaque version est') ; disp (std(differences)) ; % 5. calcul du coeficient de correlation residus = Z-Z_matlab ; r2 = (var(Z)-var(residus))/var(Z) ; disp(sprintf('coefficient de correlation r2 = %f', r2)) ;
Et pour aller plus loin ?

ca devient tellement facile d'ajuster n'umporte quoi sur autre chose qu'à un moment, on a besoin de limites. "Je suis déja au degré 13 avec 105 paramètres. Je n'ai qu'une ligne à modifier pour passer au degré 14. Jusqu'ou s'arrêtera-t-il ?

Il y a pour cela une méthode rigoureuse qui psse par le calcul d'une grandeur appelée F. F représente la quantité de variance expliquée par un modèle meilleur, par rapport à la variance déja expliquée par un modèle moins bon.

J'expliquerait ça une autre fois dans un autre mini tuto...