Code pour curve fitting avec un modèle connu

Citation:

---

Envoyé par Ayvik

EDIT:

Code:

---

1 2 3 4 5 6 7 8 9 10 11 12

>> result = myFunction ([-2.521e+006 -0.03662 -0.0004335 -0.0009894 -0.0001161 0.0003684],x,y)   result =   2.6402e+009   >> result = myFunction ([-2.521e+006 -0.03662 -0.0004335 -0.0009894 -0.0001161 0.0003684],(x-mean(x))/std(x),y)   result =   1.0837e+007

---

Il semble y'avoir un problème dans mon Matlab . . .

---

En fait j'ai refait les calculs avec cftool et j'obtiens ceci :
b = 0.0195 (0.01924, 0.01976)
c = 0.0006975 (0.000319, 0.001076)
d = 0.002046 (0.001681, 0.002411)
e = -0.0008303 (-0.001197, -0.0004635)
f = -0.0001406 (-0.0005166, 0.0002354)
a = 3.543e+06 (3.543e+06, 3.543e+06)

Je crois que tu as du prendre de mauvais résultat car les paramètres que tu donnes ne correspondent pas à la courbe... (j'ai fait un plot et ça ne correspond pas, mais alors pas du tout!)

Ah oui je n'avais pas vu le FMINSEARCHBND
ah ba non ok c'est une fonction du site, pas une fonction intégré à matlab!

Ohlala, ça veut dire que je ne me sers même pas correctement de la toolbox ??
Voilà, ce que je fais :
- Curve fitting toolbox
- Data/ je spécifie des données
- Fitting/ je spécifie mon équation dans l'onglet "Linear equations"
- J'applique et j'obtiens ceci :

http://img24.imageshack.us/img24/2816/capturevhq.jpg

Attends, tu te rends compte que ce graph n'a plus grand chose à voir avec le graph que tu as mis en première page!!

Dans le fichier .mat tu as mis une matrice nommée file avec 4 colonnes.
moi j'ai regardé le graph de la première page et pour rester cohérent avec ce graph j'ai mis la première colonne dans x, et la quatrième dans y, apparemment toi tu as pris la deuxième colonne pour y...

quels sont les données à prendre en compte?!?

Ok, :D

En fait j'ai une vingtaine de fichiers .xls organisés comme suit :

Colonne 1: une date au format numérique
Colonne 2: nombre au format . . . peu importe je crois
Colonne 3: nombre au format . . . peu importe je crois
Colonne 4: nombre au format . . . peu importe je crois

Pour chacuns de ces .xls, je dois donc produire 3 graphs sous matlab.
Le set de données que tu as doit contenir 2 colonnes si je ne m'abuses ?
Col 1: dates en x
Col 2: valeurs en y
Au début du post j'utilisais AZU.xls et après j'ai changé d'exemple, puisque ce n'est pas tant la forme des données qui m'importe, mais bien la différence dans la façon dont les courbes sont calculées ( Toolbox & fminsearch ).

Edit:
Désolé pour mon manque de cohérence, je refais en suivant la méthode que tu as suivie. Col1 = x et Col4 = y. J'update dans une minute. Cependant, je ne trouve pas le post précédent où j'ai spécifié y = a (:,4)

P.S.: Merci pour ta patience . . . ;)

ben, non ton fichier contient 4 colonnes :

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

>> file   file =   1.0e+06 *   0.0369 -2.5211 -4.6573 3.5426 0.0369 -2.5211 -4.6573 3.5426 0.0369 -2.5211 -4.6573 3.5426 0.0369 -2.5211 -4.6573 3.5426 0.0369 -2.5211 -4.6573 3.5426 0.0369 -2.5211 -4.6573 3.5426 0.0369 -2.5211 -4.6573 3.5426 0.0369 -2.5211 -4.6573 3.5426 0.0369 -2.5211 -4.6573 3.5426 0.0369 -2.5211 -4.6573 3.5426 0.0369 -2.5211 -4.6573 3.5426 0.0369 -2.5211 -4.6573 3.5426 0.0369 -2.5211 -4.6573 3.5426 0.0369 -2.5211 -4.6573 3.5426 0.0369 -2.5211 -4.6573 3.5426 0.0369 -2.5211 -4.6573 3.5426 0.0369 -2.5211 -4.6573 3.5426 0.0369 -2.5211 -4.6573 3.5426 0.0369 -2.5211 -4.6573 3.5426 0.0369 -2.5211 -4.6573 3.5426 0.0369 -2.5211 -4.6573 3.5426 0.0369 -2.5211 -4.6573 3.5426 0.0369 -2.5211 -4.6573 3.5426 0.0369 -2.5211 -4.6573 3.5426 0.0369 -2.5211 -4.6573 3.5426 0.0369 -2.5211 -4.6573 3.5426 0.0369 -2.5211 -4.6573 3.5426 0.0369 -2.5211 -4.6573 3.5426 0.0369 -2.5211 -4.6573 3.5426 0.0369 -2.5211 -4.6573 3.5426 0.0369 -2.5211 -4.6573 3.5426 0.0369 -2.5211 -4.6573 3.5426 0.0369 -2.5211 -4.6573 3.5426 0.0369 -2.5211 -4.6573 3.5426 0.0369 -2.5211 -4.6573 3.5426 0.0369 -2.5211 -4.6573 3.5426

---

et j'ai mis la première dans x et la quatrième dans y....

Oui, c'est juste mon fichier contient 4 colonnes, mais je croyais qu'il n'allait exporter que les 2 que j'avais définies en x et y . . . no comment. :oops:

So let's be clear, on a les mêmes données :
Col1 = x
Col4 = y

Avec fminsearch j'obtiens :
http://img141.imageshack.us/img141/3...sfminearch.png

Avec la toolbox :
http://img24.imageshack.us/img24/4773/cftoolbox.png

La toolbox me donne les mêmes paramètres que toi, à savoir :

Code:

---

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

Linear model: f(x) = b*x + c*sin(2*pi*x) + d*cos(2*pi*x) + e*sin(4*pi*x) + f*cos(4*pi*x) + a where x is normalized by mean 3.769e+004 and std 467.3 Coefficients (with 95% confidence bounds): b = 0.0195 (0.01924, 0.01976) c = 0.0006975 (0.000319, 0.001076) d = 0.002046 (0.001681, 0.002411) e = -0.0008303 (-0.001197, -0.0004635) f = -0.0001406 (-0.0005166, 0.0002354) a = 3.543e+006 (3.543e+006, 3.543e+006)   Goodness of fit: SSE: 0.04322 R-square: 0.9325 Adjusted R-square: 0.9323 RMSE: 0.005259

---

Et j'ai toujours un truc fou 8O :

Code:

---

1 2 3 4 5

>> result = myFunction ([3.543e+006 0.0195 0.0006975 0.002046 -0.0008303 -0.0001406],(x-mean(x))/std(x),y)   result =   2.4890e+008

---

Dans le cas précédent j'avais noté que bien que l'on tape la même chose, nos résultats diffèrent :

Citation:

---

Envoyé par magelan

Code:

---

1 2 3

result = myFunction ([-2.521e+006 -0.03662 -0.0004335 ... -0.0009894 -0.0001161 0.0003684], ... (x-mean(x))/std(x),y)

---

On se retrouve avec un résultat qui est effectivement minimisé (0.0432)

---

Citation:

---

Envoyé par Ayvik

Code:

---

1 2 3 4 5 6

>> result = myFunction ([-2.521e+006 -0.03662 -0.0004335 -0.0009894 -0.0001161 0.0003684],(x-mean(x))/std(x),y)   result =   1.0837e+007

---

Il semble y'avoir un problème dans mon Matlab . . .

---

Comment est-ce possible ? :(

Oui, mais en fait vu que tu manipules des gros chiffres avec des petits, les décimales de tes nombres ont de l'importance, et dans l'affichage matlab ne t'en montre des arrondis, donc si toi tu récupères ces arrondis, forcément à la fin les calculs sont faussés...

pour faire le calcul, j'ai exporté le résultat de cftool (save to workspace)
puis :

Citation:

---

result = myModel2 ([fittedmodel1.a fittedmodel1.b fittedmodel1.c fittedmodel1.d fittedmodel1.e fittedmodel1.f],x2,y)

---

Ca va mieux ! Donc tout ce qu'il affiche est arrondi, je prends note ...

Code:

---

1 2 3 4 5 6

>> result = myFunction ([WRHSzfm.a WRHSzfm.b WRHSzfm.c WRHSzfm.d WRHSzfm.e WRHSzfm.f],((x-mean(x)))/std(x),y)   result =   0.0432

---

Bon déjà, c'est un résultat cohérent. :)
Mais bon, là, on l'a un peu aidé hein !

EDIT: Du nouveau,

Code:

---

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

>> filename = 'WRHS.xls'; file= xlsread (filename); x= file (:,1); y= file (:,4); scatter (x,y); hold;   bestcoeffs=fminsearch(@myFunction,[0 0 0 0 0 0],[],x,y); yfit=bestcoeffs(1)+ bestcoeffs(2)*x + bestcoeffs(3)*sin(2*pi*x) + bestcoeffs(4)*cos(2*pi*x)+ bestcoeffs(5)*sin(4*pi*x) + bestcoeffs(6)*cos(4*pi*x); %Now compare y with yfit plot(x,yfit,'r'); Current plot held   Exiting: Maximum number of function evaluations has been exceeded - increase MaxFunEvals option. Current function value: 0.047641

---

MaxFunEvals
MaxFunEvals, a field closely related to MaxIter, specifies the maximum number of loglikelihood objective function evaluations. The default value is 100 times the number of parameters estimated in the model. For example, the default model has four parameters, so the default value of MaxFunEvals for the default model is 400. When the estimation process terminates prematurely, it is usually because MaxFunEvals, rather than MaxIter, is exceeded. You can increase MaxFunEvals if you suspect that the estimation terminated prematurely.

Je fais :

Code:

---

MaxFunEvals = 600

---

J'obtiens :(:

Code:

---

1 2 3

Exiting: Maximum number of function evaluations has been exceeded - increase MaxFunEvals option. Current function value: 0.047641

---

C'est vraiment bizarre, parce qu'il suffit que je redémarre Matlab, et relance mon script, et ça fonctionne, sans retourner d'erreur.. . . .. :mur:

Citation:

---

Envoyé par Ayvik

EDIT: Du nouveau,

Code:

---

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

>> filename = 'WRHS.xls'; file= xlsread (filename); x= file (:,1); y= file (:,4); scatter (x,y); hold;   bestcoeffs=fminsearch(@myFunction,[0 0 0 0 0 0],[],x,y); yfit=bestcoeffs(1)+ bestcoeffs(2)*x + bestcoeffs(3)*sin(2*pi*x) + bestcoeffs(4)*cos(2*pi*x)+ bestcoeffs(5)*sin(4*pi*x) + bestcoeffs(6)*cos(4*pi*x); %Now compare y with yfit plot(x,yfit,'r'); Current plot held   Exiting: Maximum number of function evaluations has been exceeded - increase MaxFunEvals option. Current function value: 0.047641

---

MaxFunEvals
MaxFunEvals, a field closely related to MaxIter, specifies the maximum number of loglikelihood objective function evaluations. The default value is 100 times the number of parameters estimated in the model. For example, the default model has four parameters, so the default value of MaxFunEvals for the default model is 400. When the estimation process terminates prematurely, it is usually because MaxFunEvals, rather than MaxIter, is exceeded. You can increase MaxFunEvals if you suspect that the estimation terminated prematurely.

Je fais :

Code:

---

MaxFunEvals = 600

---

J'obtiens :(:

Code:

---

1 2 3

Exiting: Maximum number of function evaluations has been exceeded - increase MaxFunEvals option. Current function value: 0.047641

---

---

Pour modifier MaxFunEvals, tu dois le faire dans les options, par : exemple :

Code:

---

options = ('MaxFunEvals',600);

---

et ensuite quand tu appel fminsearch :

Code:

---

bestcoeffs=fminsearch(@myFunction,[0 0 0 0 0 0],[],x,y);

---

----------------------------------------------------

Quoi qu'il en soit, j'ai un peu regardé de mon côté, et tu devrais plutôt te pencher sur la fonction lsqcurvefit.
Les deux méthodes tentent de minimiser le carré de la différence entre ton modèle et l'expérimental, sauf que lsqcurvefit est particulièrement dédié à ce genre de problème.
De plus j'ai réussi à obtenir les mêmes résultat que cftool avec lsqcurvefit, en prenant un vecteur x centré et normalisé et en changeant l'algorithme utilisé par lsqcurvefit (il vaut mieux prendre l'algorithme Levenberg-Marquardt qui est bien plus rapide et performant dans ton cas, pour le changer, ça se passe dans optimset comme pour MaxFunEvals).

Pour écrire l'algo, si tu suis bien le lien que tu avais donné, ça ne devrait pas te poser trop de difficultés.

Voila, yapluka :mrgreen:

Bon bha ça précise plutôt bien ma solution :D
Je vais gratter de ce coté-là alors . . .
J'avais essayé hier de l'adapter (lsqcurvefit), d'après la recommandation de Caro-Line, sans succès. Mais je vais m'y pencher de plus près ;)

C'est cool !
Merci pour tout Magelan !

:ccool:

De rien!
bon courage pour la suite :D

Pour la version avec lsqcurvefit, j'ai :

myFunction.m:

Code:

---

1 2 3 4 5 6 7 8 9

function output= myFunction (param,X) a = param(1); b = param(2); c= param(3); d= param(4); e= param(5); f= param(6); % this is the 3rd order polynomial equation here output = a + b.*X+ c.*sin(2*pi*X)+ d.*cos(2*pi*X)+ e.*sin(4*pi*X)+ f.*cos(4*pi*X);;

---

Code:

Code:

---

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23

filename = 'WRHS.xls'; file= xlsread (filename); x= file (:,1); y= file (:,2);   lb = [ -10 -10 -10 -10]; %define the lower bound ub =[10 10 10 10]; %define the upper bound options = optimset('lsqcurvefit'); options = optimset(options,'LevenbergMarquardt' ,'on'); %options = optimset(options,'Algorithm','levenberg-marquardt'); options.MaxFunEvals = 600;   [newParameters,error] = lsqcurvefit(@myFunction, [1 1 1 1 1 1],(x-mean(x))/std(x),y,lb,ub,options);   figure scatter(x,y) %plot the scatter plot hold %hold the figure   %use new parameters to get new output values y2 = myFunction(newParameters,(x-mean(x))/std(x));   %plot the new data using the color red plot(x,y2,'r')

---

En executant ce code, je n'ai pas d'erreur, cependant comme tu peux le voir en comment, j'ai trouvé 2 syntaxes pour paramétrer l'algorithme utilisé. Celle en commentaire produisait une erreur.

J'ai l'impression que mon code est correct, cependant le graphique est . . . comment dire . . . :D . . . pas juste !

http://img18.imageshack.us/img18/7051/figlsq.png

En ayant paramétré TolFun à 1.10-10, ça me transforme la ligne rouge en une courbe rouge, et j'ai :

Code:

---

1 2	Maximum number of function evaluations exceeded; increase options.MaxFunEvals

---

Et plus j'augmente MaxFunEvals:

Code:

---

options.MaxFunEvals = 10000;

---

Plus la courbe se caractérise comme une sinusoïde . . .
http://img18.imageshack.us/img18/3662/sinusiode.png

:arf:

Graphiquement, vu de loin, c'est presque ça :mouarf:

Tu n'as pas besoin de définir lb et ub (de plus tu les as mal défini, lb et ub doivent être de la taille de param (là tu contrains juste les 5 premiers paramètres à être entre -10 et 10), et trouver des résultats entre -10 et 10, ce n'est pas très cohérent avec les paramètre que tu connais déjà).

Mais bon, ça marche aussi sans borne.

Pour optimset, je crois que j'avais mis qu'une seule option:
options = optimset('Algorithm','levenberg-marquardt');

Je ne voulais pas définir les bornes spécialement au début, mais je ne sais pas si cette syntaxe est correcte :

Code:

---

[newParameters,error] = lsqcurvefit(@myFunction, [1 1 1 1 1 1],(x-mean(x))/std(x),y,options);

---

Et je ne comprends pas comment le y que je spécifie ici, est utilisé dans le processus, vu que mon myFunction commence comme ça :

Code:

---

function output= myFunction (param,X)

---

Il y'a une erreur dans mon code là ?
Pour que le tout fonctionne il faut que j'envoie, mon X et mon Y dans la fonction pour qu'il puisse déterminer mes coeffs ?
L'exemple d'avant m'a embrouillé un peu, j'arrive pas à saisir où se fait le calcul des moindres carrés là-dedans !

EDIT

Citation:

---

Envoyé par magelan

Graphiquement, vu de loin, c'est presque ça :mouarf:

---

:ccool: Je touche au but ! :D

Quand j'utilise ça:

Citation:

---

Envoyé par magelan

Pour optimset, je crois que j'avais mis qu'une seule option:
options = optimset('Algorithm','levenberg-marquardt');

---

Matlab me retourne ceci:

Code:

---

1 2	??? Error using ==> optimset at 215 Unrecognized parameter name 'Algorithm'.

---

Citation:

---

Envoyé par Ayvik

Je ne voulais pas définir les bornes spécialement au début, mais je ne sais pas si cette syntaxe est correcte :

Code:

---

[newParameters,error] = lsqcurvefit(@myFunction, [1 1 1 1 1 1],(x-mean(x))/std(x),y,options);

---

Et je ne comprends pas comment le y que je spécifie ici, est utilisé dans le processus, vu que mon myFunction commence comme ça :

Code:

---

function output= myFunction (param,X)

---

Il y'a une erreur dans mon code là ?
Pour que le tout fonctionne il faut que j'envoie, mon X et mon Y dans la fonction pour qu'il puisse déterminer mes coeffs ?
L'exemple d'avant m'a embrouillé un peu, j'arrive pas à saisir où se fait le calcul des moindres carrés là-dedans !

EDIT

:ccool: Je touche au but ! :D

---

Le code est juste :
Tu as maintenant une fonction myFunction (param,X) où param représente les coeficients que tu cherches et X est un vecteur comprenant tes données expérimentales en x.

LSQCURVEFIT cherche à minimiser :
(myFunction (param,X)-Y).^2
(regarde le help pour plus de renseignement)
mais cette fois-ci le calcul de myFunction (param,X)-Y).^2 ne se fait plus dans myFunction comme précédemment mais c'est bien lsqcurvefit qui fait le calcul (pas très clair cette phrase, mais j'espère que tu as compris).

Donc, ça devrait marcher, maintenant, non?

Edit :

Citation:

---

Matlab me retourne ceci:

Code:

---

1 2	??? Error using ==> optimset at 215 Unrecognized parameter name 'Algorithm'.

---

---

Ca dépend peut-être de ta version de matlab...

Malheureusement non ! :cry:

Mon myFunction est le même que précédemment.
Mon code est :

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

filename = 'WRHS.xls'; file= xlsread (filename); x= file (:,1); y= file (:,2);   lb = [ -10 -10 -10 -10]; %define the lower bound ub =[10 10 10 10]; %define the upper bound %options = optimset('lsqcurvefit'); options = optimset('Algorithm','levenberg-marquardt'); %options = optimset(options,'LevenbergMarquardt' ,'on'); %options = optimset(options,'Algorithm','levenberg-marquardt'); options.MaxFunEvals = 10000;   [newParameters,error] = lsqcurvefit(@myFunction, [1 1 1 1 1 1],(x-mean(x))/std(x),y,options);   figure scatter(x,y,'.') %plot the scatter plot hold %hold the figure   %use new parameters to get new output values y2 = myFunction(newParameters,(x-mean(x))/std(x));   %plot the new data using the color red plot(x,y2,'r')

---

Avec le code que tu sembles avoir utilisé:

Code:

---

options = optimset('Algorithm','levenberg-marquardt');

---

J'obtiens:

Code:

---

1 2	??? Error using ==> optimset at 215 Unrecognized parameter name 'Algorithm'.

---

Avec le code, trouvé sur un site web:

Code:

---

options = optimset(options,'Algorithm','levenberg-marquardt');

---

J'obtiens:

Code:

---

1 2	??? Error using ==> optimset at 215 Unrecognized parameter name 'Algorithm'.

---

Avec cet autre code:

Code:

---

options = optimset(options,'LevenbergMarquardt' ,'on');

---

Cela fonctionne et produit le graph altéré. Mais seulement si j'ai les limites de spécifiées.

Code:

---

[newParameters,error] = lsqcurvefit(@myFunction, [1 1 1 1 1 1],(x-mean(x))/std(x),ub,lb,y,options);

---

Une fois les bornes supprimées :

Code:

---

[newParameters,error] = lsqcurvefit(@myFunction, [1 1 1 1 1 1],(x-mean(x))/std(x),y,options);

---

J'ai :

Code:

---

1 2	??? Error using ==> lsqcurvefit at 171 LSQCURVEFIT only accepts inputs of data type double.

---

Pour le coup, il semble que la syntaxe sans les bornes, mais avec les options ne soit pas correct, enfin de la façon dont je l'ai écrite du moins :cry:

Yé né sé plu ! :D

Citation:

---

Envoyé par Ayvik

Je ne voulais pas définir les bornes spécialement au début, mais je ne sais pas si cette syntaxe est correcte :

Code:

---

[newParameters,error] = lsqcurvefit(@myFunction, [1 1 1 1 1 1],(x-mean(x))/std(x),y,options);

---

[/COLOR][/CODE]

---

Il faut quand même préciser que tu ne mets pas de borne :

Code:

---

[newParameters,error] = lsqcurvefit(@myFunction, [1 1 1 1 1 1],(x-mean(x))/std(x),y,[],[],options);

---

-----------------------------------

EDIT
Il faut des majuscules dans le nom:
options = optimset('Algorithm','Levenberg-Marquardt');
EDIT2 : ah ba non, chez moi, ça marche aussi sans majuscule...

:D Tellement évident que j'y n'y avait pas pensé . . . :oops:

Après correction du nombres de bornes, et modification sus-dite:

Code:

---

1 2 3 4 5 6 7 8 9

lb = [ -1 -1 -1 -1 -1 -1]; %define the lower bound ub =[1 1 1 1 1 1]; %define the upper bound options = optimset('lsqcurvefit'); %options = optimset('Algorithm','levenberg-marquardt'); options = optimset(options,'LevenbergMarquardt' ,'on'); %options = optimset(options,'Algorithm','levenberg-marquardt'); options.MaxFunEvals = 10000;   [newParameters,error] = lsqcurvefit(@myFunction, [1 1 1 1 1 1],(x-mean(x))/std(x),y,[],[], options);

---

J'ai toujours un truc comme ça :
http://img199.imageshack.us/img199/5702/yadumieux.png

Cette courbe qui maintenant semble être la bonne, semble juste être centrée autour de zéro . . .
J'ai raison de faire ça ?

Code:

---

y2 = myFunction(newParameters,(x-mean(x))/std(x));

---

PS: Très important de spécifier correctement les bornes pour ne pas les utiliser ! 8-)

Il m'a demandé d'augmenter MaxIter, après avoir augmenté MaxFunEvals.
Le tout fonctionne, mais me donne toujours ce même graphique.

Code:

---

1 2 3 4 5

options = optimset('lsqcurvefit'); options = optimset(options,'LevenbergMarquardt' ,'on'); options.MaxFunEvals = 5000; options.MaxIter = 600; [newParameters,error] = lsqcurvefit(@myFunction, [1 1 1 1 1 1],(x-mean(x))/std(x),y,[],[], options);

---

Je n'arrive pas à savoir pourquoi mes deux données (origine et courbe produite) sont décalés, comme une simple translation apparement, bien que le zoom par défaut donne cette impression de données "aplaties" . . .