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.