Discussion :

[PatDo10]

Bonjour,
Je ne sais pas si mon message est dans le bon forum.

Je dois donc programmer, en c++ avec Qt et openCv au besoin, un algorithme qui permettra de faire une approximation de courbe non linéaire.

Je m'explique:
J'ai un vecteur de points, les X sont représentés par des TI (un certains temps en ms) et les Y représente le signal S associé à chaque TI. J'ai donc des points et je voudrais faire une courbe ajustée.
J'ai une équation, S(TI) = S0 (1 - exp(-TI/T1)), je ne connais pas T1 et S0. Il faudrait donc que je détermine ces variables, à la manière du solveur d'excel par exemple. Pour ensuite recalculer chaque signal afin d'avoir une courbe ajustée.
Je ne sais pas comment m'y prendre et je ne comprends pas bien le fonctionnement de l'algorithme de Levenberg Marquardt malgré mais recherche sur internet.
Si quelqu'un pourrait m'aider ça serait très gentil.

Merci d'avance.
Cordialement.

[wiwaxia]

Bonjour,

... J'ai donc des points et je voudrais faire une courbe ajustée.
J'ai une équation, S(TI) = S0 (1 - exp(-TI/T1)), je ne connais pas T1 et S0. Il faudrait donc que je détermine ces variables ...

Pour commencer, tu devrais écrire la fonction paramétrée sous une forme plus appropriée:
F(x, a, b) = a*(1 - Exp(b*x))
dans laquelle a = S0 et b = 1/T1 .
Tu auras sans doute besoin des dérivées partielles correspondantes:
D_aF = (1 - Exp(b*x)) ; D_bF = -a*x*Exp(b*x) .

Il faut envisager la somme des carrés des différences
S = Som_k=1^N(y_k - F(x_k))²
dont tout minimum est caractérisé par l'annulation simultanée des dérivées partielles:
D_aS = 0 ; D_bS = 0 .
On obtient dans le cas présent, et après suppression du facteur (2), les équations:

Som_k=1^N(y_k - F(x_k))*(1 - Exp(b*x_k)) = 0 ;

Som_k=1^N(y_k - F(x_k))*a*x_k*Exp(b*x_k) = 0 .

Quel est le nombre total (N) de points ?
Quel est le domaine des abscisses ? Sont-elles toutes positives ?
Leur intervalle s'approche-t-il suffisamment de l'origine, et a-t-on Min(x_k) << Max(x_k) ? ou Min(x_k) << (1/b) = T1 ?
Observe-t-on un palier - (y_k) quasi-constant aux grades valeurs des abscisses ?

Voilà qui permettrait de mieux cerner le problème numérique.

Au sujet du procédé évoqué (algorithme de Levenberg-Marquard), peut-être recevrais-tu plus de réponses sur le forum Statistiques et Data Mining.

[PatDo10]

Tout d'abord merci beaucoup pour ta réponse

Ensuite est ce que la fonction paramétrée ne serait pas plutôt :
F(x, a, b) = a*(1 - Exp(-b*x))

Du coup les dérivées partielles seraient :
DaF = (1 - Exp(-b*x))
DbF = -a*x*Exp(-b*x)


Mon nombre total de points est de 8 avec des valeurs de x toujours positives.
Voici les points que j'obtiens (capture d'écran)



Je n'ai pas compris les autres questions désolé

Merci d'avance

[wiwaxia]

... ce que la fonction paramétrée ne serait pas plutôt :
F(x, a, b) = a*(1 - Exp(-b*x))

Du coup les dérivées partielles seraient :
DaF = (1 - Exp(-b*x))
DbF = -a*x*Exp(-b*x)

C'est exact, le signe m'avait échappé.
Un changement de signe intervient pour la seconde dérivée:
D_bF = +a*x*Exp(-b*x)

... Mon nombre total de points est de 8 ...

Alors pas besoin de se creuser la tête: on entreprend la recherche du minimum de S(a, b).

C'est par contre très décevant, en ce qui concerne la disposition des points, censés se rassembler autour d'un graphe présentant une asymptote horizontale:



Le traitement numérique reste possible, mais les résultats seront entachés d'une énorme incertitude.

Pire encore: je vois à l'instant que F(x, a, b) devrait être une fonction croissante de (x) - à moins d'admettre un coefficient (a) négatif:
D_xF = a*b*Exp(-b*x) - et qu'il devrait intervenir un troisième terme constant:
F(x, a, b, c) = a*(1 - Exp(-b*x)) + c .

La situation me paraît bloquée: il te faut reprendre soit la relation théorique, soit les données numériques.
Courage ...

[Jerome Briot]

je ne comprends pas bien le fonctionnement de l'algorithme de Levenberg Marquardt malgré mais recherche sur internet.

Voir peut être ici : Ajustements d'équations non linéaires avec des exemples MATLAB

[Invité]

bj,

pour ma part j'ai tout trouvé sur wiki https://en.wikipedia.org/wiki/Levenb...ardt_algorithm

vu que j'ai fait ca par curiosité...si ca peut servir
(le nom des variables suit wiki de fait...)

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
 
import numpy as np
def linear_solve(a,b):
    #not what is important
    return np.matmul(np.linalg.inv(a), b)
 
def S(f, y, x, beta):
    v = y-f(x, beta)
    return np.dot(v,v)
 
def run(y, x, f, jac):
    MAX_ITER = 50
    MAX_K = 20
    eps = 1e-7
 
    n_steps = 0
    err = 1e5
 
    lambd = 1 #should be 0.01 according to https://damien-forthomme.developpez.com/tutoriels/matlab/ajustement-equations-non-lineaires/ 
    mu = 2
 
    beta = [1,1]
    J = computeJacobian(x, beta)
    I = np.eye(len(beta))
    while n_steps<MAX_ITER and err > eps:
        left = np.matmul(J.T,J)+lambd*I
        right = np.dot(J.T, y-f(x, beta))
        delta = linear_solve(left, right)
        err1 = S(f,y,x, beta+delta)
 
        lambd2 = lambd / mu
        left2 = np.matmul(J.T,J)+lambd2*I
        delta2 = linear_solve(left2, right)
        err2 = S(f,y,x, beta+delta2)
 
        if err2 < err:
            lambd = lambd2
            err = err2
            beta += delta2
        elif err1 > err:
            lambdk = lambd2
            i = 0
            while i < MAX_K:
                lambdk = lambdk * mu
                leftk = np.matmul(J.T,J)+lambdk*I
                deltak = linear_solve(leftk, right)
                errk = S(f,y,x, beta+deltak)
                if errk < err:
                    lambd = lambdk
                    err = errk
                    beta += deltak
                    break
            if i == MAX_K:
                print('FAILED to reduce error')
        else:
            err = err1
            beta += delta
 
        n_steps += 1
        print('err : ', err, 'beta', beta, 'lambda', lambd)
    return [err, beta, n_steps]
 
 
def f(x, beta):
    return beta[0]*(1-np.exp(-x/beta[1]))
 
def computeJacobian(x, beta):
    J = np.zeros((len(x), 2))
    [a, t1] = beta;
    dfda = 1 - np.exp(-x/t1)
    dfdt1 = -a/(t1*t1)*np.exp(-x/t1)
    J[:,0] = dfda
    J[:,1] = dfdt1
    return J
x=np.array([1,2,3,4,5,6,7,8])
y=f(x, [5,3])
[err, beta, n_steps] = run(y=y, x=x, f=f, jac=computeJacobian)
#we expect beta to be [5,3]...

output

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
 
err :  2.6064559169149604 beta [3.72642092 1.70908193] lambda 0.5
err :  1.1536161441877706 beta [4.20474051 2.42671229] lambda 0.25
err :  0.567696931476128 beta [4.58666269 2.90653567] lambda 0.125
err :  0.20629510300362114 beta [4.87570217 3.15342765] lambda 0.0625
err :  0.028870090309115207 beta [5.04863703 3.18937032] lambda 0.03125
err :  0.010401188799042827 beta [5.10069331 3.09424863] lambda 0.015625
err :  0.00978856504109014 beta [5.06808875 3.0348687 ] lambda 0.125
err :  0.006206422242097047 beta [5.03073602 2.98871621] lambda 0.0625
err :  0.0014755456545637981 beta [5.00027612 2.97273909] lambda 0.03125
err :  0.00012160259413442261 beta [4.98679749 2.98284017] lambda 0.015625
err :  0.00010334527301431457 beta [4.9884182  2.98645611] lambda 0.5
err :  9.416785799788405e-05 beta [4.99074049 2.99151212] lambda 0.25
err :  8.33810386447726e-05 beta [4.99388766 2.99704601] lambda 0.125
err :  4.939921187477068e-05 beta [4.99734552 3.00115299] lambda 0.0625
err :  1.1018041768446205e-05 beta [5.0000594  3.00247847] lambda 0.03125
err :  1.0459898580773698e-06 beta [5.00121378 3.00152233] lambda 0.015625
err :  8.442290813224142e-07 beta [5.00103989 3.00119752] lambda 0.5
err :  7.51575485580231e-07 beta [5.00082327 3.00074819] lambda 0.25
err :  6.574937889473048e-07 beta [5.00054081 3.00025851] lambda 0.125
err :  3.8678653268172065e-07 beta [5.00023362 2.99989605] lambda 0.0625
err :  8.53050124303239e-08 beta [4.99999354 2.99978028] lambda 0.03125

ca a l'air ok

(je me suis pas embarassé avec les variables temporaires)


edit: apres relecture ci-dessus c'est juste levenberg, j'ai pris I au lieu de diag(J^tJ) pour le calcul de left

[Nebulix]

Comme l'a remarqué wiwaxia, la fonction proposée ne peut pas être la bonne. Avant de te lancer dans la programmation il est impératif de réfléchir aux mécanismes en jeu.

Il n'y a que deux points significativement différents de l'asymptote, et le deuxième est plus haut que le premier. Il y a un bruit très important dans ton signal ce qui rend une exploitation quantitative hasardeuse.

Le mieux qu'on puisse tirer de ces données est qu'il y a quelque chose pour t<=200 et rien pour t>500.

Aucun traitement informatique ne peut compenser des données de mauvaise qualité. Ni le fait qu'on n'a pas compris le phénomène.

Pour la programmation, quand tu en seras arrivé là, je te conseille le "downhill simplex" que tu trouveras dans Numerical Recipes.

[PatDo10]

Je ne peux donc pas utiliser cette méthode pour trouver T1 et S0 ?

J'ai vraiment du mal à comprendre

[Invité]

Code :

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

 Pour la programmation, quand tu en seras arrivé là, je te conseille le "downhill simplex" que tu trouveras dans Numerical Recipes.  *

en tt cas d'apres
http://docs.mantidproject.org/v3.7.2...inimizers.html

marquardt levenberg semble qd même plus accurate et plus rapide
mais bon je connaissais pas dc c'est cool

Code :

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

J'ai vraiment du mal à comprendre

ya pas 56 trucs à comprendre.
On t'a dit que ta fonction qd bien même tu lui trouves les meilleurs T1 et S0, si tu la traces, la courbe passe pas trop par tes points. t'as juste à regarder le graphique que wiwaxia t'as copié..

[Nebulix]

Je ne peux donc pas utiliser cette méthode pour trouver T1 et S0 ?

J'ai vraiment du mal à comprendre

Je suppose que tu faisais référence à mon post.
AUCUNE méthode ne peut tirer de paramètres d'un jeu de données si l'information pertinente n'est pas dans les données ou est brouillée par un bruit trop important (ou autres pb).
AUCUNE méthode informatique ne peut comprendre le phénomène que tu veux étudier et proposer une fonction qui puisse le décrire mathématiquement.
Dans les deux cas, un programme te donnera des valeurs ... qui ne voudront rien dire et pourront entrainer de gros dégâts si tu y crois.

[wiwaxia]

PatDo10 est manifestement déçue et choquée qu'un algorithme pointu, et difficile à mettre en oeuvre, ne puisse conduire à de bons résultats.
Cela vient pourtant d'être très clairement exprimé: la relation sur laquelle on travaille peut provenir d'un autre domaine des mathématiques, ou des sciences expérimentales - elle est en tout cas étrangère à la statistique, qui ne se prononce que sur la précision avec laquelle la relation en question est vérifiée.

Exemple emprunté à la physique: dans l'équilibre liquide-vapeur d'un corps pur, les variations de la pression de vapeur saturante en fonction de la température sont représentées par une relation du type:
Log(P) = A - B/T + C*Ln(T) .
Ce n'est pas le fruit d'une lubie des expérimentateurs, mais une conséquence des lois de la Thermodynamique: une étude statistique d'une série de données (T_k, P_k) portera sur cette relation, à l'exclusion de toute autre - par exemple
P = A + B*T + C*T² plus conforme à l'intuition, mais dépourvue d'intérêt.

De plus, le procédé statistique peut s'appliquer à tout nuage de points; mais les résultats auxquels il peut conduire n'ont aucun sens si la forme du nuage est incompatible avec celle du graphe de la fonction utilisée - et c'est malheureusement ce qui se passe ici.

Réfléchis un peu: le graphe d'équation y = a(1 - Exp(-x/b))
# passe par l'origine F(0, a, b) = 0
# présente une pente positive D_xF = (a/b)*Exp(-x/b) (si a et b sont positifs)
# et un palier horizontal Lim(F)_{(quand x tend vers l'inf.)} = a .
Comment peux-tu espérer le superposer à ton nuage de points ?
Prend s'il le faut une feuille de papier et un crayon (ce n'est pas déshonorant ...)

[PatDo10]

Je veux principalement trouver S0 et T1, si la courbe ne passe pas par la plupart des points c'est pas grave.

Je vais essayer de bien comprendre l'algorithme de Levenberg Marquardt et je reviendrais vers vous si j'ai des questions.

Merci à tous.

[tbc92]

8 points, c'est peu. Plus tu as de points, plus tu auras des résultats corrects.

L'autre truc, c'est que tu pars d'un postulat : Ces 8 points sont proches d'une courbe de la forme y = S0 (1 - exp(-x/T1))

D'où vient ce postulat ? Es-tu sûr de cette étape ?

A partir des 8 points en question, on peut par exemple imaginer une courbe de la forme y = a+ b/(x+c) et une courbe de cette forme collerait très bien à tes points.

[PatDo10]

C'est ce genre de courbe que je cherche à avoir :

[tbc92]

Donc on a bien mis le doigt sur un problème.
La courbe que tu as dessinée, son équation n'est pas de la forme y = S0 (1-exp(x/T1))

Fais l'expérience, prend des valeurs de S0 et de T1 comme tu veux, et dessines les courbes correspondantes y = S0 (1-exp(x/T1)) ; tu obtiendras toujours des courbes croissantes, ou bien décroissantes pour d'autres valeurs de S0 et de T1, mais tu n'obtiendras jamais des courbes avec une partie plus ou moins rectiligne, puis une partie 'asymptotique' comme sur ton dessin.

[wiwaxia]

On peut même incorporer la première branche descendante en considérant que le graphe est celui de la fonction:
y = F(x) = Abs(a*(1 - 2*Exp(-x/b))) , qui vérifie: F(0) = a et LimF_{(quand x tend vers l'inf.)} = a .
Les pentes des tangentes au point de raccordement sont identiques en valeur absolue.

L'inversion de signe des ordonnées des 2 premiers points permet de rattacher l'ensemble à une équation unique:
G(x) = a*(1 - 2*Exp(-x/b)) , dont découle la précédente par la relation: F(x) = Abs(G(x)) .

PS; le graphe présente en abscisse la mention "Temps d'inversion"; il y a donc bien en ordonnée une grandeur qui change de signe, et dont on s'est contenté de représenter la valeur absolue. Vouloir intégrer les deux premiers points à l'ensemble sans aucune transformation relève de l'aberration.

[wiwaxia]

Ça saute aux yeux.

[wiwaxia]

On peut se limiter à la partie du graphe d'ordonnée positive (ce qui est du gaspillage) à condition de prendre comme nouvelle origine le point d'intersection avec l'axe horizontal.

Ce point vérifie G(x₀) = 0 = a*(1 - 2*Exp(-x₀/b))
d'où: Exp(x₀/b) = 2
et par conséquent: G(x) = a*(1 - Exp(x₀/b)*Exp(-x/b)) = a*(1 - Exp((x₀ - x)/b)) ;

Si l'on pose: x' = x - x₀ , il vient: G(x) = a*(1 - Exp(-x'/b)) = H(x') . CQFD

# Et si des esprits pointilleux arguent du fait que F(0) n'est peut_être pas égal à (-a), il suffit de passer à une équation plus générale de la forme :
G(x) = a*(1 - K*Exp(-x/b)) , qui conduit par les mêmes considérations à un résultat analogue:
H(x') = a*(1 - Exp(-x'/b)) , avec x' = x - x₀ et Exp(x₀/b) = K .
Seul ennui: il y a cette fois 3 paramètres à ajuster, au lieu de deux.
C'est peut-être bien la raison secrète de l'élimination des deux premiers points.

Et l'on reste toujours tributaire de la localisation incertaine du point d'intersection ...

[wiwaxia]

On peut envisager le traitement statistique des (N - 1) points, sur la même relation générale
y = G(x) = a*(1 - K*Exp(-x/b)) , à condition de se référer au premier; ce dernier joue ainsi un rôle privilégié - mais c'est le prix à payer pour en revenir à une fonction à deux paramètres ajustables.

Le premier couple de coordonnées vérifie en effet: y₁ = a*(1 - K*Exp(-x₁/b)) ,
et il vient dans ces conditions:
(y/a) - 1 = -K * Exp(-x/b)
(y₁/a) - 1 = - K*Exp(-x₁/b) , ce qui donne: (y - a)/(y₁ - a) = Exp((x₁ - x)/b) ;

si l'on pose: x' = x - x₁ , y' = y - y₁ et a' = a - y₁ , il vient tout naturellement:

y' = (y - a) - (y₁ - a) = (y₁ - a)*Exp(-x'/b) - (y₁ - a) = (a - y₁)*(1 - Exp(-x'/b)) ,

soit finalement: y' = a'*(1 - Exp(-x'/b)) .

Autrement dit, toute portion du graphe se ramène à la même équation, moyennant un changement approprié d'origine.

Elle est tout de même sympa, la branche d'exponentielle !

[wiwaxia]

Revenons au nuage initial.



Le remplacement des deux premiers points par leur symétrique ne permet pas un démarrage immédiat des calculs, compte tenu de leur disposition: vis-à-vis du graphe moyen aisément visualisable , et qui passe entre les deux, ils apparaissent aussi mauvais l'un que l'autre.

Le moindre mal consiste donc à les remplacer par leur milieu (M₁₂), de coordonnées négatives:
x_m = (x₁ + x₂)/2 , y_m = (y₁ + y₂)/2
avant de procéder au changement de coordonnées:
x' = x - x_m , y' = y - y[SUB]m[/SUB
pour les six points restants.