Discussion :

[membreComplexe12]

Salut tous,

je cherche à faire des régressions avec des modèles simples mais je ne sais pas trop comment choisir le modèle et l'identifier (j'ai fais plusieurs tests qui me paraissait correctes mais la courbe obtenue n'est pas bonne...)

Tout d'abord, voici mes deux courbes que je cherche à identifier :
(cf. PJ image plus données pour la première courbe)

Je vous explique ce que je dois respecter :
1°) je veux utiliser polyfit pour identifier mes modeles car fonction très rapide
2°) il faut que le modèle est une évolution monotone. Donc pour la premiere courbe il faut que la dérivée du modèle soit toujours positive est sur le second modele toujours négative.

Voici ce que j'ai fais :
1°) j'ai choisi un modèle Y=A.X^B
ce qui me donne : ln(Y)=ln(A)+B.ln(X)
=> du coup, j'ai fais

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
p=polyfit(log(X),log(Y),1)
A=exp(p(2));
B=p(1);
modele=A*X.^B;
figure(1)
plot(X,modele)

2°) j'ai choisi un modèle Y=A.ln(B.X)
ce qui me donne : exp(Y)=exp(A)+B.X
=> du coup, j'ai fais

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
p=polyfit(X,exp(Y),1)
A=log(p(2));
B=p(1);
modele=A.log(B.X);
figure(1)
plot(X,modele)

3°) j'ai choisi un modèle Y=A.exp(B.X)
ce qui me donne : ln(Y)=ln(A)+B.X
=> du coup, j'ai fais

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
p=polyfit(X,ln(Y),1)
A=exp(p(2));
B=p(1);
modele=A.exp(B.X);
figure(1)
plot(X,modele)

Avec un de ces trois modèles je m'attendais à avoir un résultat correct pour ma première courbe mais "que nenni"

=> du coup, je ne sais plus trop quoi faire... je ne sais pas si le problème vient de mes modèles où de ma programmation...

[le fab]

pas de réponse à proprement parler sur ton soucis de fitting (faut peut être essayé des fonction plus complexes), par contre il y a plus rapide que polyfit, surtout pour un polynôme du 1er degré, en utilisant \ :
remplace

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

p=polyfit(log(X),log(Y),1)

par

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

p=log(X)\exp(Y);

[membreComplexe12]

Bonjour le Fab et merci pour cette premiere réponse

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

p=polyfit(log(X),log(Y),1)

par

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

p=log(X)\exp(Y);

en effet, je vois le principe (il faut que je fasse sur papier pour m'en persuader). Merci pour l'astuce.

=> ma courbe est particulièrement simple donc je suis un peu étonné que des modèles simples ne fonctionnent pas...

Si quelqu'un peux me donner son avis là dessus...
merci

[Jean Dumoncel]

Bonjour,

Avec un de ces trois modèles je m'attendais à avoir un résultat correct pour ma première courbe mais "que nenni"

pas très explicite, tu aurais pu développer un peu les résultats que tu obtenais...

Ta première valeur est (0,0) donc si tu prends le log tu obtiens (-Inf,-Inf) d'ou les NaN que tu dois recevoir de polyfit.

Comme Le Fab ton modèle me semble un peu simpliste, je te conseille de montrer tes courbes dans le forum algorithme pour savoir quelle fonction pourrait s'ajuster sur ce type de courbe.

[membreComplexe12]

Merci tous d'avoir pris le temps de répondre

Ta première valeur est (0,0) donc si tu prends le log tu obtiens (-Inf,-Inf) d'ou les NaN que tu dois recevoir de polyfit.

as oui, en effet. Je viens de modifier mon code pour prendre en compte ceci

Bonjour,
pas très explicite, tu aurais pu développer un peu les résultats que tu obtenais...

ah oui, en effet. Désolé.
Ci dessous voici le code que j'utilise (il est utilisable directement si vous voulez voir ce que ça donne il suffit juste d'avoir le fichier data.txt dans le dossier de travail)

Code :

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
%% fit
data=load('data.txt');
X=data(:,1);
Y=data(:,2);
%% modele puissance
p=polyfit(log(X(2:end)),log(Y(2:end)),1)
A=exp(p(2));
B=p(1);
modele=A*X.^B;
figure(1)
plot(X,Y,'b')
hold on
plot(X,modele,'*r')
%% modele log
p=polyfit(X,exp(Y),1)
A=log(p(2));
B=p(1);
modele=zeros(1,length(X));
modele(2:end)=A*log(B*X(2:end));
figure(2)
plot(X,Y,'b')
hold on
plot(X,modele,'*r')
%% modele exponentiel
p=polyfit(X(2:end),log(Y(2:end)),1)
A=exp(p(2));
B=p(1);
modele=A*exp(B*X);
figure(3)
plot(X,Y,'b')
hold on
plot(X,modele,'*r')

ce code me donne trois courbes (1,2 et 3) doit voici en PJ l'image

1°) la première courbe me donne un fit pas trop loin de la réalité mais il ne colle pas parfaitement aux données. Du coup, ça ne me va pas.

2°) le deuxieme me donne des réssultats énormes : très très loin de la réalité. De plus j'ai un warning qui me dit qu'il y a des nombres complexes dans mon truc

3°) la dernière courbe donne un truc qui ressemble à quelque chose mais la courbure est inversée

Comme Le Fab ton modèle me semble un peu simpliste, je te conseille de montrer tes courbes dans le forum algorithme pour savoir quelle fonction pourrait s'ajuster sur ce type de courbe.

je viens de poser des questions sur leur forum. Merci du conseil

A+

[le fab]

au vue de tes point pourquoi tu ne fite pas un simple polynôme ?
une droite collerait mieux que tes 3 tentatives, et si tu veux t'en rapprocher encore plus un polynôme du 2ème degré ou plus

ta deuxième tentative pourrait être pas mal avec un déplacement vertical :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

modele=A.log(B.X)+Y(1)

[membreComplexe12]

merci d'avoir pris le temps de repondre

au vue de tes point pourquoi tu ne fite pas un simple polynôme ?
une droite collerait mieux que tes 3 tentatives, et si tu veux t'en rapprocher encore plus un polynôme du 2ème degré ou plus

- en fait avec un polynome de degrès 2 je n'arrive pas à coller parfaitement à la courbe.
- Pour un polynome de degrès plus important par contre je colle parfaitement à la courbe.

Par contre le soucis avec le polynome c'est que localement il peut y avoir de petites oscillations et du coup ça ne colle pas vraiment avec ce que je veux faire :
=> je veux que la dérivée en chaque point soit positive....

ta deuxième tentative pourrait être pas mal avec un déplacement vertical :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

modele=A.log(B.X)+Y(1)

je testerai, merci.

=> par contre pourquoi me parle t il de nombres complexes ?

au vue de tes point pourquoi tu ne fite pas un simple polynôme ?
une droite collerait mieux que tes 3 tentatives, et si tu veux t'en rapprocher encore plus un polynôme du 2ème degré ou plus

ta deuxième tentative pourrait être pas mal avec un déplacement vertical :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

modele=A.log(B.X)+Y(1)

[ol9245]

Ta courbe n'est pas un polynôme de X car elle est asymptotique. Les polynomes ne le sont pas.
Si tu tiens à la faire rentrer au chausse pied dans un polynome, voila ce qu'il faut faire :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
k0 = min(X) ;
k1 = max(X) ;
Xk = (X-k0)/k1 ; % pour éviter des overfow
 
% on cherche maintenant une régression entre Xk et Y
powers = [0 2 4 6 8 10 12 14 16 18 20] ;
U = bsxfun(@power, Xk, powers) ;
K = U\Y ;
 
Y2 = (U*K) ;
figure(1)
hold off
plot(Xk, Y2, '-')
hold on 
plot(Xk, Y, '.')

Le problème c'est que même en prenant des puissances paires, tu n'a aucune garantie que les coefficients soient tous positifs. la preuve :

Code :

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
K>> K
 
K =
 
   1.0e+07 *
 
    0.0000
    0.0003
   -0.0038
    0.0337
   -0.1779
    0.5853
   -1.2254
    1.6291
   -1.3292
    0.6063
   -0.1183

Conclusion : ta fonction n'est vraiment pas un polynome de X.
Si tu tiens à conserver les polynomes poour des raisons algorithmiques (il est joli le petit code hein ), alors il faut que tu transforme tes y (chnagement de variable) de manière à ce qu'ils aillent à l'infini (pour singer un vrai polynome, que tu vas ajuster ensuite). Selon les familles de courbes sur lesquelles tu as à travailler, il y a pas mal de solutions pour le changement de variable.

en voila une au pif (changement de variable par tangente. tout est dedans même les données tu peux copier le tout dans un fichier vierge, le sauver et le jouer.

Code :

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
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
function my_fit2
M = data ; % lecture de tes données dans une matrice M à deux colonnes
X = M(:,1) ;
Y = M(:,2) ;
k0 = min(X) ;
k1 = max(X) ;
k2 = max(Y) ;
Xk = X2K(X, k0,k1) ;
Yk = Y2K(Y, k2) ;
% on cherche maintenant une régression entre Xk et Yk
powers = 0:2:10 ;
 
U = bsxfun(@power, Xk, powers) ;
K = U\Yk ; % le modèle
Y2 = (U*K) ; % les valeurs simulées
 
figure(1)
hold off
plot(X, K2Y(Y2, k2), '-')
hold on 
plot(X, Y, '.') 
end
 
% les fonctions de changement de variable en X et en Y
 
function X = X2K(X, k0, k1)
X = (X-k0)/k1 ;
end
 
function X = K2X(X, k0, k1)
X = X.*k1 + k0 ;
end
 
function Y = Y2K(Y, y0)
Y = tan(Y * (pi/2*0.9 / y0)) ;
end
 
function Y = K2Y(Y, y0)
Y = atan(Y) ./  (pi/2*0.9 / y0);
end
 
% Les données
function M = data
M=[0.000000  0.000000
0.000054  3.651103
0.000104  7.686931
0.000155  12.064052
0.000218  15.651276
0.000273  20.024326
0.000333  23.893334
0.000393  27.667778
0.000451  31.897439
0.000506  35.934300
0.000556  39.924951
0.000611  43.862777
0.000643  46.967480
0.000711  51.964196
0.000771  56.134343
0.000834  60.010673
0.000890  64.257914
0.000943  68.258370
0.001000  71.762958
0.001061  76.102560
0.001115  80.074136
0.001172  83.800245
0.001221  88.377166
0.001273  92.318190
0.001338  95.957944
0.001385  99.886904
0.001443  103.607786
0.001517  107.931473
0.001565  111.906755
0.001614  115.816621
0.001670  119.670741
0.001716  123.884518
0.001772  127.976461
0.001829  131.715339
0.001892  135.723343
0.001956  139.241143
0.002011  143.270729
0.002063  147.554154
0.002120  151.224619
0.002168  155.449189
0.002209  159.103544
0.002290  163.281087
0.002332  167.176223
0.002383  171.311151
0.002460  174.313960
0.002515  178.763084
0.002565  182.445918
0.002617  186.285859
0.002661  190.659420
0.002726  194.440866
0.002785  198.482577
0.002834  202.389011
0.002900  206.279309
0.002964  210.308885
0.003010  213.543672
0.003081  217.667689
0.003109  221.520643
0.003167  225.786957
0.003227  229.381494
0.003282  233.410481
0.003353  237.254472
0.003395  240.912387
0.003450  244.819457
0.003512  248.422991
0.003554  251.914270
0.003618  255.743586
0.003673  259.421120
0.003740  263.104048
0.003791  266.179981
0.003846  269.455373
0.003897  272.625880
0.003940  275.424635
0.004010  277.720568
0.004066  280.165750
0.004138  281.885816
0.004178  283.458683
0.004235  284.696281
0.004283  285.961878
0.004335  286.827765
0.004389  287.512115
0.004465  287.919197
0.004508  288.581952
0.004576  289.222779
0.004617  289.621007
0.004685  289.662900
0.004732  290.623260
0.004805  290.913975
0.004848  290.890326
0.004912  291.345561
0.004949  291.419170
0.005018  291.932728
0.005066  292.157341
0.005137  292.180274
0.005189  292.676951
0.005239  292.616698
0.005297  293.153295
0.005344  293.136046
0.005399  293.540013
0.005469  293.656930
0.005527  293.712079
0.005576  293.756993
0.005630  293.955046
0.005675  294.092305
0.005742  294.487700
0.005799  294.488345
0.005861  294.869446
0.005911  294.580046] ;
end

Ca reste du chasse pied.
Si ça te suffit pas, il faut un vrai modèle intelligent de ta variable et l'utilisation d'un ajustement non linéaire du type fminfun

[membreComplexe12]

Merci beaucoup d'avoir pris le temps de te pencher sur mon problème c'est très gentil

Je vais étudier ce que tu viens me donner afin de bien saisir tous les éléments.

Merci

ps: je crois qu'à présent j'ai tous les éléments pour arriver à faire quelque chose

[le fab]

très joli code ol !!
je connaissais pas bsxfun (plus sympa que de le faire à la main ) ...

sinon tu pourrais expliquer un peu

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

Xk = (X-k0)/k1 ; % pour éviter des overfow

je vois bien la différence en le faisant ou pas, mais j'ai du mal à comprendre

j'ai déjà fais qqch dans le genre, mais juste

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

Xk = X-k0

et ça me paraissait suffisant
(je précise que je faisais une régression linéaire sur un polynôme d'ordre 1 : dérivée glissante, et c'est surement pour ça que ça marche)

question subsidiaire :
comment tu récupères tes coefficients du polynôme fonction de X (plutot que Xk)

[ol9245]

très joli code ol !!
je connaissais pas bsxfun (plus sympa que de le faire à la main ) ...

sinon tu pourrais expliquer un peu

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

Xk = (X-k0)/k1 ; % pour éviter des overfow

je vois bien la différence en le faisant ou pas, mais j'ai du mal à comprendre

j'ai déjà fais qqch dans le genre, mais juste

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

Xk = X-k0

et ça me paraissait suffisant
(je précise que je faisais une régression linéaire sur un polynôme d'ordre 1 : dérivée glissante, et c'est surement pour ça que ça marche)

question subsidiaire :
comment tu récupères tes coefficients du polynôme fonction de X (plutot que Xk)

Merci du compliment. J'apprécie
En fait mon attrait pour le beau code me prend plus de temps que je devrais normalement y consacrer. Donc à la fin je suis à la bourre et le code qui va en production est de la bouillie . Le pire ennemi du beau code, c'est le beau code !!

pour les polynomes, je met tous les X dans une boite 0-1, comme ça toutes les puiisances nièmes sont dans la même fourchette, donc toutes les erreurs d'arrondi ont la même magnitude. Cela n'est utile qu'à partir du moment où il y a plus de 14 ordres de grandeur entre la plus petite puissance (en général le chiffre 1) et la plus grande. Ca va très vite. Par exemple 1000^5 est déjà hors des clous.

Si tu as besoin des coefficinet du polynôme en X et non pas en Xk (pour quoi ?), il y a des règles de math qui s'appliquent. Si j'en avais besoin, je chercherai dans l'ordre :
* aide mémoire de mathématiques, Bronstein & Semendiaev, Eyrolles, 934pp
* calcul formel sous Mapple ou Matlab,
* Google
* inversion du problème avec Matlab ?

mais bon, je n'en vois aps l'usage, à part de faire des erreurs d'arrondi dans les calculs. Mon habitude est de passer par Xk puis d'inverser le changement de variable initial.

[le fab]

pourquoi
bonne question

en fait la seule régression linéaire que j'ai codé jusqu'à aujourd'hui est sur un polynôme d'ordre 1 (dérivée par régression linéaire sur n points glissants)
on est loin des 14, passer par une boite 0-1 n'est donc pas utile

cependant comme j'applique ma fonction sur un vecteur de plus de 300 000 points (issu d'une simulation sur un an dans simulink), j'ai été contraint de faire la moitié du travail, cad travailler sur un Xk = X - min(X) pour ramener les X à l'"origine" du repère temporelle et éviter des erreurs d'arrondis sur des nombres très grands

quand j'ai vu ton (beau) code, je me suis dis que j'allais faire de même et travailler sur un Xk = (X - min(X))/max(x) (même si au vue de tes dernières explications c'est loin d'être utile sur mon polynôme d'ordre 1)

or ce qui m’intéresse dans ma dérivée, c'est la pente, donc le coefficient a du polynôme a.x+b
d’où ma question "comment récupérer les coefficients ?"

mais bon, vu la complexité de mon probleme, ma solution initiale de ramener mes n points glissants à une origine est suffisante pour éviter les pb d'arrondis et plus simple pour récupérer a sans rajouter des risques d'erreur

m'enfin, c'était surtout pour ma culture générale

[membreComplexe12]

régression linéaire sur n points glissants

juste pour ma culture qu'es ce que ceci veux dire ?
=> si on a un vecteur de "n" composantes tu interpoles les composantes 1:10 10:20 20:30..... jusqu'à "n" c'est ça une regression glissante ?

Au fait j'ai une autre question (encore pour ma culture) :
=> savez vous ce qu'es le principe des moindres carré mobiles ?
j'ai souvant entendu se terme mais je ne sais pas vraiment qu'es ce que c'est

[le fab]

juste pour ma culture qu'es ce que ceci veux dire ?
=> si on a un vecteur de "n" composantes tu interpoles les composantes 1:10 10:20 20:30..... jusqu'à "n" c'est ça une regression glissante ?

plutôt 1:n puis 2:n+1 puis 3:n+2 jusqu'à la fin
et ce n'est pas de l'interpolation mais de la régression : pour chaque plage on cherche la droite qui passe au mieux au milieu des points

moindres carrés je connais, mobiles, jamais entendu parlé

[ol9245]

pourquoi
bonne question

en fait la seule régression linéaire que j'ai codé jusqu'à aujourd'hui est sur un polynôme d'ordre 1 (dérivée par régression linéaire sur n points glissants)
on est loin des 14, passer par une boite 0-1 n'est donc pas utile

cependant comme j'applique ma fonction sur un vecteur de plus de 300 000 points (issu d'une simulation sur un an dans simulink), j'ai été contraint de faire la moitié du travail, cad travailler sur un Xk = X - min(X) pour ramener les X à l'"origine" du repère temporelle et éviter des erreurs d'arrondis sur des nombres très grands

quand j'ai vu ton (beau) code, je me suis dis que j'allais faire de même et travailler sur un Xk = (X - min(X))/max(x) (même si au vue de tes dernières explications c'est loin d'être utile sur mon polynôme d'ordre 1)

or ce qui m’intéresse dans ma dérivée, c'est la pente, donc le coefficient a du polynôme a.x+b
d’où ma question "comment récupérer les coefficients ?"

mais bon, vu la complexité de mon probleme, ma solution initiale de ramener mes n points glissants à une origine est suffisante pour éviter les pb d'arrondis et plus simple pour récupérer a sans rajouter des risques d'erreur

m'enfin, c'était surtout pour ma culture générale

au degré 1, si tu divises tes X par k, tu multiplie la pente par k. C'est assez visuel/intuitif et facile à vérifier expérimentalement. Pour revenir à la vraie pente : diviser la pente calculée par ce même k. Là, pas besoin d'un aide mémoire...

C'est dès le degré 2 que tu as des coefficients qui sont distribués aux degrés inférieurs donc ça demande au moins papier-crayon et plus éventuellement.