Code pour curve fitting avec un modèle connu

**Ayvik** · 22/07/2009, 19h28

Bonjour tout le monde,

Cela fait maintenant un sacré bout de temps que j'essaie d'écrire un code basique pour produire un graphique avec mes données d'origines, et la courbe représentant au mieux ces données. Je connais déjà la forme de l'équation, et donc il me faut trouver la valeur des six inconnues, pour la meilleur courbe.

Mes données sont dans un fichier excel.
Jusqu'a présent j'ai trouvé plusieurs exemples, mais après avoir essayé de les transposer à mon application . . . ça marche pas

Pour l'instant j'ai ça :

Code :

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

1
2
3
4
5
6
7
8
9
10
filename = 'input.xls';
 a= xlsread (filename);
 x= a (:,1);
 y= a (:,2);
 
 bestcoeffs=fminsearch(@fun,[1 1 1 1 1 1],[],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);
 
 scatter (x,y);
 plot(x,yfit);

Ce qui est adapté d'un exemple censé marcher.
Matlab me retourne :

??? Undefined function or method 'fun' for input arguments of type 'double'.

Error in ==> fminsearch at 205
fv(:,1) = funfcn(x,varargin{:});

En tant que bon newbie, je ne comprends pas ce qui ne va pas . . .
J'ai déjà utilisé la toolbox "Curve fitting", et ça fonctionne plutôt bien, mais j'aimerai arriver au même résultat avec un bout de code !

Si vous avez une idée . . .
Merci d'avance !

Ayvik.

**Caro-Line** · 22/07/2009, 19h33

As-tu regardé l'aide de FMINSEARCH pour comprendre ce que cette fonction fait ?
Là il te dit simplement que la fonction fun que tu lui passes en paramètre n'existe pas.

**Ayvik** · 25/07/2009, 18h27

J'ai effectivement potassé la section sur fminsearch dans la doc en ligne de Matlab, ai testé plusieurs fonctions, mais je suis sur de devoir utiliser fminsearch maintenant.

Le nom de la fonction dont tu parles devrait être reconnu, puisque j'ai un fichier .m avec ma fonction déclarée dedans. Le problème c'est que je ne trouve pas d'exemple sur la syntaxe a utiliser quand on veut passer beaucoup de paramètres à la fonction.

J'ai un set de donnée X et Y dans un xls. et je connais le modèle qui régit la relation, cependant je doit trouver la valeur des 6 paramètres, et ce, pour la meilleur courbe. Le but serait d'arriver à reproduire le resultat fourni par la toolbox, mais en le codant.

Ma fonction est :
a + b.x + c.sin(2.pi.x) + d.cos(2.pi.x) + e.sin(4.pi.x) + f.cos(4.pi.x)

J'ai créé un fichier .m avec ce code, trouvé dans un exemple:

Code :

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

1
2
3
4
5
6
7
8
9
function f = myModel (x,input)
a = param(1);
b = param(2);
c=  param(3);
d=  param(4);
e=  param(5);
f=  param(6);
% this is the model equation here
output = a + b.*input.^1+ c.*sin(2*pi*input.^1)+ d.*cos(2*pi*input.^1)+ e.*sin(4*pi*input.^1)+ f.*cos(4*pi*input.^1);

Je ne suis évidemment pas du tout sur de la validité du code, vu que je ne sais pas ce que fais !

(J'aurai pas du faire confiance à mon tuteur sur les objectifs :-S ) De plus, j'ai très peu de temps, trop peu pour passer un mois à apprendre les bases de la programmation en tout cas. J'essaie donc de me baser sur des exemples.

En ce qui concerne, le code en lui même, j'ai ça (toujours adapté d'un exemple) :

Code :

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

1
2
3
4
5
6
7
8
9
10
filename = 'CCCS.xls';
 a= xlsread (filename);
 x= a (:,1);
 y= a (:,2);
 
 bestcoeffs=fminsearch(@myModel,[1 1 1 1 1 1],[],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);
 
 scatter (x,y);
 plot(x,yfit);

En executant ce code, Matlab me retourne :

??? Error using ==> myModel
Too many input arguments.

Error in ==> fminsearch at 205
fv(:,1) = funfcn(x,varargin{:});

Je ne sais pas comment lui faire accepter tous les arguments que je lui passe. Et en ce qui concerne la deuxième erreur, je n'ai aucune idée de la fonction dont il parle, vu que je ne lui ai pas demandé de l'utiliser . . .

Toute aide me serait vraiment vraiment précieuse . . .
Ayvik,

P.S: J'ai du faire de la prog 'régulière' (CShell) il y'a quelque jours, et je trouve que les ressources sur Matlab sont vraiment . . . nulles, surtout pour les vrais débutants ! Aucun exemple vraiment simple, ou d'explication claire des exemples . . . Je suppose que le top est d'avoir des notions de prog avant de s'attaquer au Matlab !

**Caro-Line** · 25/07/2009, 18h59

Sans notion de programmation cela va être vraiment compliqué.

1. Ta fonction myModel
A pour paramètres d'entrée x et input (soit 2 entrées)
Dès la 1ère tu essayes d'utiliser une variable param qui n'est mentionnée nulle part => cela ne peut pas fonctionner.
Comme tu n'utilises x nulle part je suppose qu'en fait param=x et/ou inversement.
Dans un 1er temps écris donc correctement ta fonction myModel.
Pour cela teste-la d'abord.
Par exemple là si tu fais :

Code :

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

myModel(1,2)

ça va planter tout de suite.

2. Utilisation de FMINSEARCH
Le 1er argument est le nom de la fonction : maintenant ça c'est bon il trouve bien myModel, par contre le message d'erreur te dit que tu as un souci au niveau du nombre d'arguments d'entrée passés : il y en a trop.
Le 2ème argument de FMINSEARCH correspond à des paramètres d'entrée de la fonction myModel (la/les valeurs à partir desquelles FMINSEARCH fait la recherche du minimum local) =>
Tu lui en donnes 6 alors que myModel n'a que 2 arguments d'entrée => la problème vient donc de là.

Nota : FUNFCN est une fonction privée de MATLAB qui vérifie la cohérence du nombre d'arguments.

**Ayvik** · 25/07/2009, 19h32

Merci Caro-line pour ta réponse.

Effectivement, il n'y avait pas assez de paramètres déclarés.
Je me sers de ces deux exemples pour faire l'adaptation :
Code en général:
http://blinkdagger.com/matlab/matlab...-lsqcurvefit/2
Fminsearch:
http://www-h.eng.cam.ac.uk/help/tpl/...e_fitting.html

Est-ce une bonne idée ?

J'ai modfié en ceci,qui semble être correct d'après l'exemple 1.

Code :

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

1
2
3
4
5
6
7
8
9
function output = myModel (param,input)
a = param(1);
b = param(2);
c=  param(3);
d=  param(4);
e=  param(5);
f=  param(6);
% this is the model equation here
output = a + b.*input.^1+ c.*sin(2*pi*input.^1)+ d.*cos(2*pi*input.^1)+ e.*sin(4*pi*input.^1)+ f.*cos(4*pi*input.^1);

Et Matlab me retourne :

??? Error using ==> times
Matrix dimensions must agree.

Error in ==> myModel at 9
output = a + b.*input.^1+ c.*sin(2*pi*input.^1)+ d.*cos(2*pi*input.^1)+ e.*sin(4*pi*input.^1)+
f.*cos(4*pi*input.^1);
Error in ==> fminsearch at 205
fv(:,1) = funfcn(x,varargin{:});

J'espère qu'il est capable de multiplier X par 2pi ! :-)

Est-ce que mes deux sources ont l'air correctes ?

Merci,

EDIT:

J'ai finalement ré-adapté le tout, donc : (j'espère que tout ça servira à quelqu'un . . . même si très basique)

Ma fonction myModel.m :

Code :

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

1
2
3
4
5
6
7
8
9
function out = myModel (coeff,X,Y)
a = coeff(1);
b = coeff(2);
c=  coeff(3);
d=  coeff(4);
e=  coeff(5);
f=  coeff(6);
% this is the model equation here
Y_fun = a + b.*X+ c.*sin(2*pi*X)+ d.*cos(2*pi*X)+ e.*sin(4*pi*X)+ f.*cos(4*pi*X);

Mon bout de code :

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
 filename = 'CCCS.xls';
 a= xlsread (filename);
 x= a (:,1);
 y= a (:,2);

 bestcoeffs=fminsearch(@myModel,[1 1 1 1 1 1],[],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);

 scatter (x,y);
 hold;
 plot(x,yfit,'r');

Le tout ne produit pas d'erreur à l' éxecution, mais me produit un graphique totalement différent de celui produit par la tool box.

Produit par la "Curve fitting" Toolbox :

Produit par fminsearch :

Il semble que la fonction ne produise une courbe qu'avec les faibles oscillations, comme si c.sin(2.pi.x) + d.cos(2.pi.x) n'etaient pas pris en compte . . .

La question s'éloigne de son point de départ . . .

Any ideas ?
Merci beaucoup . . .

**Caro-Line** · 25/07/2009, 19h46

Le message d'erreur indique un problème dans les dimensiosn des matrices multipliées.

As-tu essayé de lancer ta function myModel avec en arguements d'entrée ce que tu vas lui donner avec FMINSEARCH ?

Peux-tu montrer les dimensions de tes variables x et y ?

**Ayvik** · 25/07/2009, 20h08

X et Y étant deux colonnes de mon fichier excel, chacune fait 1540*1.

X est une date, au format numérique, e.g. 36891
Y est un chiffre comme :- 2472979.37630603
Le format des nombres n'a pas d'incidence à ce que je crois comprendre, les dimensions des matrices oui. Mais je ne vois quelle matrice il multiplient entre elles.

Voir l'update de mon post précédent, le code semble fonctionner, mais il doit y'avoir quelque chose à paramétrer en plus je pense . . . Je ne comprends pas pourquoi ça ne produit qu'un résultat 'altéré'.

**Caro-Line** · 25/07/2009, 23h20

Tu as essayé avec LSQCURVEFIT plut^t que FMINSEARCH comme proposé dans ton 1er lien ?

**Ayvik** · 26/07/2009, 00h31

J'avais essayé précedemment, et ça produisait le même resultat je crois . . .
Enfin ca ne donnait pas la courbe obtenue avec la toolbox.

J'essaie de ré-adapter mon cas à l'example lsqcurvefit, mais sans succès.
Et bizarrement, ce qui marchait tout à l'heure avec fminsearch, ne marche plus !

Anyway, il ne me semble pas avoir vu de paramètres supplémentaires pour la fonction fminsearch . . . mis à part :
Ces quelques paramètres

Je vais re-vérifier ce que donne lsqcurvefit. Et au-dela de ça . . .
J'editerai le post pour spécifier ce que je j'obtiens.

Ciao !

EDIT:

[U]

En ce qui concerne fminsearch, j'ai réussi à obtenir ce que j'avais eu au post précédent, en revanche pour lsqcurvefit . . . Je n'ai pas réussi à l'adapter.
J'essaie de trouver quelqu'un qui saurait m'expliquer le phénomène décrit plus haut.

Ciao !

**Jean Dumoncel** · 26/07/2009, 14h40

Salut,

dans ton code qui utilise fminsearch, une fois que bestcoeffs a été calculé, que retourne la ligne suivante :

Code :

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

result = myModel (bestcoeffs,x,y)

**Ayvik** · 26/07/2009, 16h31

Mon myFunction.m est comme suit, à l'instant présent :
calqué sur l'exemple précité, sous la section "Arbitrary functions".

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
13
function out= myFunction(coeff,X,Y)
 a = coeff(1);
 b = coeff(2);
 c = coeff(3);
 d=  coeff(4);
 e=  coeff(5);
 f=  coeff(6);
 % this is the model equation here
 Y_fun = a + b.*X+ c.*sin(2*pi*X)+ d.*cos(2*pi*X)+ e.*sin(4*pi*X)+ f.*cos(4*pi*X);
 DIFF = Y_fun - Y; 
 SQ_DIFF = DIFF.^2;
 
 out = sum(SQ_DIFF);

J'ai testé l'exemple du lien, et il fonctionne parfaitement.
Après avoir fait l'adaptation, j'obtiens ceci :

Code :

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

1
2
3
4
5
>> result = myFunction (bestcoeffs,x,y)
 
result =
 
    0.0293

0.0293 = sum(SQ_DIFF) je suppose.

Matlab passe probablement cette valeur unique a tous mes paramètres ?

Je pensais qu'une simple substitution suffirait, mais j'ai l'air de m'être trompé dans l'adaptation . . .

Code :

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

1
2
3
4
5
6
7
>> bestcoeffs
 
bestcoeffs =
 
  1.0e+006 *
 
   -1.8003   -0.0000    1.5790   -1.3857    0.5903    0.6649

Bien que le modèle ne soit pas linéaire à proprement parler, (il est "Statistically Linear" ou aussi dit "Linear in the parameters"), je rentre mon modèle dans "Custom equations/Linear equations", et la courbe ressort montrant les tendances du nuage de points.

Pour les meme données avec la toolbox :

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
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.03662  (-0.03683, -0.03641)
       c =  -0.0004335  (-0.0007393, -0.0001278)
       d =  -0.0009894  (-0.001284, -0.0006946)
       e =  -0.0001161  (-0.0004124, 0.0001802)
       f =   0.0003684  (6.459e-005, 0.0006721)
       a = -2.521e+006  (-2.521e+006, -2.521e+006)

Goodness of fit:
  SSE: 0.02821
  R-square: 0.9868
  Adjusted R-square: 0.9867
  RMSE: 0.004248

**Jean Dumoncel** · 26/07/2009, 17h22

Envoyé par Ayvik

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
13
function out= myFunction(coeff,X,Y)
 a = coeff(1);
 b = coeff(2);
 c = coeff(3);
 d=  coeff(4);
 e=  coeff(5);
 f=  coeff(6);
 % this is the model equation here
 Y_fun = a + b.*X+ c.*sin(2*pi*X)+ d.*cos(2*pi*X)+ e.*sin(4*pi*X)+ f.*cos(4*pi*X);
 DIFF = Y_fun - Y; 
 SQ_DIFF = DIFF.^2;
 
 out = sum(SQ_DIFF);

J'ai testé l'exemple du lien, et il fonctionne parfaitement.
Après avoir fait l'adaptation, j'obtiens ceci :

Code :

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

1
2
3
4
5
>> result = myFunction (bestcoeffs,x,y)
 
result =
 
    0.0293

0.0293 = sum(SQ_DIFF) je suppose.

Oui, donc 0.0293 représente le carré de la différence entre ton Yfit et Yexperimental.
fminsearch cherche donc à minimiser cette différence, et apparemment il n'a pas été assez loin dans la minimisation. Y-a-t-il un message renvoyé par la fonction fminsearch du genre "maximum d'itération atteinte" ou autre?

J'aurais peut-être quelque idées pour améliorer les résultats mais comme je ne suis pas sûr de moi, je ne voudrais pas te donner de mauvaises indications et te faire perdre du temps....
Serait-il possible que tu envoies en pièce jointes un set (ou une partie) de ton tableau a pour que je puisse exécuter le code avec tes données?
(pour enregistrer en .mat : save matriceA a, et tu joins le fichier matriceA.mat dans un de tes messages).

**Jean Dumoncel** · 26/07/2009, 17h38

Au fait que te renvoie : result = myModel (coeff,x,y) mais cette fois en utilisant les coeff retourné par la toolbox curve fitting?

**Ayvik** · 26/07/2009, 17h54

Je viens de relancer Matlab, parce que je pense que j'avais modifié MaxIter et Maxfun.evals. J'ai executé le code, et ça fonctionne, comme décrit précémment: pas d'erreur, mais courbe pas juste.

J'ai fais un .mat de mon tableau (fichier non valide pour le mettre en pièce jointe).

Si j'ai bien compris c'est comme ça qu'on entre les coeffs, donc :

Code :

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

1
2
3
4
5
>> result = myFunction ([-2.521e+006 -0.03662 -0.0004335 -0.0009894 -0.0001161 0.0003684],x,y)

result =

  2.6402e+009

2.6402e+009 - - > Avec les "bons" coeffs de la tool box, donc censé être "minimisé" ?

Je ne comprends pas vraiment la mécanique de tout ça, j'aimerai bien faire les choses proprement, mais malheureusement je n'en ai pas le temps . . .

**Jean Dumoncel** · 26/07/2009, 18h26

Oui, ce dernier résultat m'étonne aussi...

Je vais regardé de plus près, je ne te garantis rien mais je posterais si j'arrive à comprendre le pourquoi du comment des 2 types de méthodes!

**Ayvik** · 26/07/2009, 18h29

Ca éclairerait bien ma lanterne ! Merci !

Pendant ce temps j'essaie de trouver d'autres ressources, voir si mon code est mal écrit, ou s'il y'aurait des paramètres supplémentaires à définir . . .

**Jean Dumoncel** · 26/07/2009, 18h57

Bon, déjà une petite précision :

les données en x dans la toolbox sont centrées et normalisées (par la moyenne et l'écart-type), ce qui veut dire que le vecteur x utilisé par la toolbox est en fait :

Code :

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

(x-mean(x))/std(x)

Du coup i on calcule :

Code :

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

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)

**Ayvik** · 26/07/2009, 19h08

Très juste, j'avais oublié de préciser ça . . .

Quand je cochais cette option dans la toolbox, je pensais que ça ne faisait que centrer mes données sur le graphique, et pas qu'il fallait le comprendre dans le sens mathématique.

Donc le résultat peu être considéré comme correct maintenant, et le but du jeu c'est arriver à obtenir le même avec fminsearch ? La spécifications de bornes pourraient elle aider en quoi que ce soit, mis à part la rapidité d'éxecution ? (ce dont je ne bénificierai pas trop, puisque le calcul est relativement rapide)

EDIT:

Code :

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

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 . . .

**Jean Dumoncel** · 26/07/2009, 19h17

Pour les bornes, sûrement mais je ne crois pas qu'on ait la possibilité de les spécifier dans fminsearch?

**Ayvik** · 26/07/2009, 19h22

Pour les bornes, il semble que si, mais c'est une fonction différente :

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
function [x,fval,exitflag,output]=fminsearchbnd3(fun,x0,LB,UB,options,varargin)
% FMINSEARCHBND: FMINSEARCH, but with bound constraints by transformation
% usage: x=FMINSEARCHBND(fun,x0)
% usage: x=FMINSEARCHBND(fun,x0,LB)
% usage: x=FMINSEARCHBND(fun,x0,LB,UB)
% usage: x=FMINSEARCHBND(fun,x0,LB,UB,options)
% usage: x=FMINSEARCHBND(fun,x0,LB,UB,options,p1,p2,...)
% usage: [x,fval,exitflag,output]=FMINSEARCHBND(fun,x0,...)
% 
% arguments:
%  fun, x0, options - see the help for FMINSEARCH
%
%  LB - lower bound vector or array, must be the same size as x0
%
%       If no lower bounds exist for one of the variables, then
%       supply -inf for that variable.
%
%       If no lower bounds at all, then LB may be left empty.
%
%       Variables may be fixed in value by setting the corresponding
%       lower and upper bounds to exactly the same value.
%
%  UB - upper bound vector or array, must be the same size as x0
%
%       If no upper bounds exist for one of the variables, then
%       supply +inf for that variable.
%
%       If no upper bounds at all, then UB may be left empty.
%
%       Variables may be fixed in value by setting the corresponding
%       lower and upper bounds to exactly the same value.

P.S: As-tu vu l'update de mon post précédent ? Il semble que je n'obtienne pas le même résultat que toi . . .

Code pour curve fitting avec un modèle connu [Débutant]

Signal

Discussions similaires

Partager

Partager