MATLAB: méthode de Newton, ou du gradient, (ou autre) pour une fonction à plusieurs variables

[fiatlux]

Bonjour,
J'ai une courbe de points expérimentaux (x,y) sur lesquels j'aimerais coller un modèle de courbe F dépendant de 4 paramètres: y=F(x,p1,p2,p3,p4). Par exemple F(x,p1,p2,p3,p4)=p1*x^3+p2*sin (x)+p3*x+p4. Je connais les vecteurs de points x et y, je cherche les p.
Je voudrais faire ça avec Matlab. Je sais qu'il y a des fonctions (genre lsqcurvefit) qui permettent de faire ça, mais j'aimerais:
1) pouvoir donner "plus d'importance" à certaines zones, donc par exemple fenêtrer l'erreur entre le modèle et la réalité à chaque itération pour que ce soit surtout l'erreur sur les zones importantes qui compte. Et ça je ne peux pas le faire avec lsqcurvefit
2) implémenter ça moi-même.
J'ai d'abord pensé à la méthode de Newton. Mais je vois 2 problèmes:
1) la méthode sert à minimiser une fonction. Ici je devrais donc tenter de minimiser l'erreur entre mon modèle et la réalité, mais du coup ce n'est plus une fonction, or avec la méthode de Newton il nous faut ensuite calculer analytiquement la dérivée de la fonction....
2) je sais utiliser cette méthode pour des fonctions à une variable (genre uniquement p1) mais pas quand il y en a 4... Genre comment est-ce que j'update les variables à chaque pas d'itération ?
Autre méthode possible: algorithme du gradient. Mais là aussi, comment résoudre cela avec une fonction à 4 variables au lieu d'une ?

En gros je suis un peu perdu :-S
Merci à quiconque pouvant m'aider un peu ou simplement me donner une piste.
-----

[ABN84]

Quelle est ta fonction modèle exactement?
selon qu'elle est linéaire par rapport aux paramètres p ou pas, la mmethode à choisir sera differente

[Dlzlogic]

Bonsoir,
Je ne suis pas sûr d'avoir tout compris.
1- vous parlez de points expérimentaux. Donc, on peut en déduire que les points forment un nuage. Alors, d'où vient votre fonction f(x) ?
2- on utilise dans un domaine précis une fonction de la forme
Q = K . A^a . C^c . I^i
où K, a, c, i sont des paramètres fixés d'après des éléments précis. C'est assez comparable à vos p, puisque ce sont des paramètres calculés une fois pour toutes, suivant les conditions.
et A, C, I sont les variables
Cette forme offre l'intérêt de pouvoir être calculée par un système linéaire, en prenant son logarithme.

Par ailleurs, j'ai mis au point dernièrement une méthode d'ajustement avec 3 variables, c'est à dire qu'on détermine une fonction f(x, y, z).
Sauf le complexité, on pourrait l'étendre à une fonction à 4 variables.
Je ne pense pas qu'il soit intéressant de chercher une solution itérative, mais plutôt vers l'établissement d'un système linéaire à 4 inconnues.

[fiatlux]

Bonjour, merci à vous.
La fonction modèle est en réalité F(x,p)=p1*tanh((x-p2)/p3) + p4. La raison pour laquelle j'ai ça c'est parce que mon signal a plus ou moins cette forme-là (je parlais de points mais en faits c'est simplement un signal y qui a à peu près la forme d'une tanh sur x=3:0.01:7 (par exemple)). Donc en fait j'aimerais trouver p1,p2,p3 et p4 qui font que mon modèle F(x,p)=p1*tanh((x-p2)/p3) + p4 avec x=3:0.01:7 colle au mieux à mon signal.

[A voir en vidéo sur Futura]

[ABN84]

methode de newton:

tu parts d'un p0 arbitaire:

p_i+1=x_i-f(p_i)* pinv(g(p_i))

avec f(p)=p(1)*tanh((x-p(2))/p(3)) + p(4)
et g(p)=
[-tanh((x-p(2))/p(3));
p(1)*(1-tanh((x-p(2))/p(3))^2)/p(3);
p(1)*(1-tanh((x-p(2))/p(3))^2)*(x-p(2))/p(3)^2;
-1]

[Dlzlogic]

Bonjour,
Soit cette fonctions f.
Elle comporte 4 paramètres à définir.
4 valeurs de (x,y) permettent d'établir 4 équations.
Le système à 4 inconnues se résout sans difficulté, au prix d'un ou deux changements de variable peut-être

[fiatlux]

@ABN84
Merci! C'est exactement ce que je cherchais. Donc en gros t'as remplacé la dérivée de f par le gradient ? ça paraît effectivement plus ou moins "logique" mais tu as trouvé ça où ? Parce que j'ai cherché et à part tomber sur la méthode de Newton avec 1 fonction à 1 variable ou la méthode de Newton à N fonctions à N variables, j'ai trouvé nulle part pour 1 fonction à N variables.

2ème chose: g(p) c'est bien g(p)=[df(p)/dp1 ; df(p)/dp2 ; df(p)/dp3 ; df(p)/dp4] ? Si oui d'où viennent les signes moins ? Genre dérivée par rapport à p4, pourquoi -1 et pas 1 ? Idem pour la dérivée par rapport à p1 et p2. Y'a juste celle par rapport à p3 où y'a pas ça (d'après mes calculs).

dernière chose: je vois que tu utilises pinv au lieu de inv. J'ai regardé l'help de pinv... ça ne m'inspire pas énormément. Pourquoi avoir utilisé pinv et non inv? ça marcherait avec inv ?
Merci!

@Dlzlogic
Merci, mais mon signal n'a pas exactement la forme d'une tanh, et il peut être un peu bruité. Il suffit de choisir 4 mauvais points et.... Ou alors (ça me donne une idée) j'essaie avec genre 1000 quadruplets de 4 points choisis au hasard, je calcule à chaque fois l'erreur et je regarde pour quel quadruplet elle est la plus petite ? ça peut être pas mal non ?

[ABN84]

4 valeurs de (x,y) permettent d'établir 4 équations.
Le système à 4 inconnues se résout sans difficulté, au prix d'un ou deux changements de variable peut-être

non, tu n'as pas pris en compte l'erreur de mesure. il faut avoir un système d'equation surdeterminé pour esperer rejeter toute incertitude

[fiatlux]

@ABN84
J'oubliais: puisque la méthode de newton sert à approcher un zéro d'une fonction, ça va pas marcher non ? Je ne veux pas approcher le zéro de f, je veux approcher le zéro de l'erreur entre f et mon signal plutôt... ou bien ?

[Dlzlogic]

4 valeurs de (x,y) permettent d'établir 4 équations.
Le système à 4 inconnues se résout sans difficulté, au prix d'un ou deux changements de variable peut-être
non, tu n'as pas pris en compte l'erreur de mesure. il faut avoir un système d'equation surdeterminé pour esperer rejeter toute incertitude

Parfaitement d'accord. Je voulais simplement noter que le texte précisait 3 valeurs pour x, alors qu'il en faut au moins 4, et en pratique une bonne quinzaine.
En fait, il me semble qu'on est dans le cas typique de calcul de valeurs par la méthode des moindres carrés.
Les 4 valeurs à calculer sont p1, p2, p3, p4. La valeurs les meilleures sont obtenues si la somme des carrés des écarts constatés est minimale. Ceci est obtenu en annulant la dérivée de cette somme.
En vertu de l'indépendance des écarts on peut écrire que la somme des dérivées partielles est nulle.
On peut établir ainsi un système de 4 équations où l'on trouvera des termes de la forme Somme(X), Somme(X²), Somme(XY) etc.
Un millier de couples (x,y) me parait beaucoup pour la capacité d'une machine, une cinquantaine, c'est déjà pas mal, mais rien n'empêche de traiter successivement plusieurs paquets et comparer les résultats.

[ABN84]

Bonsoir,

dernière chose: je vois que tu utilises pinv au lieu de inv. J'ai regardé l'help de pinv... ça ne m'inspire pas énormément. Pourquoi avoir utilisé pinv et non inv? ça marcherait avec inv ?
Merci!

l'inverse d'un matrice n'est definie que pour une matrice carré.
dans ton cas tu auras besoin d'une matrice 2 x n, c'est pourquoi on utilise la pseudo inverse.

J'oubliais: puisque la méthode de newton sert à approcher un zéro d'une fonction, ça va pas marcher non ? Je ne veux pas approcher le zéro de f, je veux approcher le zéro de l'erreur entre f et mon signal plutôt... ou bien ?

regarde bien la fonction dont j'ai cherché d'erreur: il s'agit de y-f(p) pas uniquement f(p). Je n'ai pas le temps de le démontrer là, mais defacto, A=Y*pinv(X) donne A qui minimise la norme L2 de l'erreur Y-A*X, d'ou l'interret de la pseudo inverse pour la methode de newton à plusieurs paramètres.
Je vérifierais demain la jacobienne de ta fonction.
++

[fiatlux]

@Dlzlogic: il me semble que justement ce système de 4 équations à 4 inconnues n'est pas si trivial que ça à résoudre analytiquement. Du coup je suis obligé de faire appel à une méthode numérique....

@ABN84
Dans ton post est-ce qu'il n'y aurait donc pas des erreurs ?
p_i+1=x_i-f(p_i)* pinv(g(p_i)) ---> x_i ? plutôt p_i non ? et f(p_i) ? pas justement (y-f(pi)) ?

avec f(p)=p(1)*tanh((x-p(2))/p(3)) + p(4) ----> ok
et g(p)=
[-tanh((x-p(2))/p(3)); ----> ok
p(1)*(1-tanh((x-p(2))/p(3))^2)/p(3); ----->ok
p(1)*(1-tanh((x-p(2))/p(3))^2)*(x-p(2))/p(3)^2; ----> manque un signe moins devant non ?
-1] ---> ok


Merci

[ABN84]

Je corrige:
p_i+1=p_i-f(p_i)* pinv(g(p_i))

avec

f=y-(p1*tanh((x-p2)/p3) + p4)
g=[
tanh((p2 - x)/p3);
-(p1*(tanh((p2 - x)/p3)^2 - 1))/p3;
(p1*(tanh((p2 - x)/p3)^2 - 1)*(p2 - x))/p3^2;
-1]

[fiatlux]

Merci.

Juste une question encore si je peux:
on est bien d'accord que p est un vecteur 1x4, x est un vecteur tel que x=3:0.01:7, donc vecteur 1x401, y vaut s(x) où s est mon signal donc vecteur 1x401, f = vecteur 1x401 et g est une matrice 4x401 ?
La pseudo inverse de g me donne une matrice 401x4. Donc le produit matriciel f(p_i)*pinv(g(p_i)) me donne un vecteur 1x4.

[ABN84]

c'est exact

[fiatlux]

Ça marche ABN84 ! Merci beaucoup pour votre aide à tous les deux !
Bonne soirée