ajustement exponentiel

el_douwen

salutations.

En pleine créations d'énoncé pour des étudiants qui apprennent R, je tombe sur ce souci que je n'ai jamais vraiment résolu.

J'ai des données (x,y) que j'ai générées par y=exp(a*x+b) mais en rajoutant du bruit autour de y, tant et si bien que y contient des valeurs négatives.

Comment, à partir des données .csv, retrouver a et b ?

Si je translate y par y=y-min(y) puis à la fin:

curve(exp(a*x+b)+min(y))

la corrélation n'est visuellement pas bonne. Logique puisque en moyenne $y$ ne tend pas vers 0 si x tend vers + l'infini

Si je fait log(y) tel quel avec les NaN engendrés au final la corrélation est mauvaise ; logique puisqu'il manque des valeurs dans le calcul de a,b

Comment faire ?

gerard0

Bonjour.

Si je comprends bien, tu veux trouver le modèle exponentiel qui s'adapte le mieux à tes données. Il n'y en a pas. Il n'y a pas de modèle exponentiel qui modélise à la fois des valeurs positives et des valeurs négatives.

Après, tu peux t'appuyer sur ton procédé de construction et rechercher a et b tels que la moyenne quadratique des erreur soit minimale. Mais alors ce n'est plus un fit automatique.

Il me semble qu'il serait plus sain de générer un échantillon de la forme $y=\exp(ax+b+\varepsilon(x))$ où le bruit est dans l'exponentielle.

Cordialement.

el_douwen

Ok je comprend bien et je ne génèrerai plus que des bruits de la forme "fois (1+ bruit)" au lieu de "+bruit" quand je ferai des énoncés.
Ma problématique, derrière, est que, aussi, je n'ai aucune connaissance du terrain, et de où trouver des données réelles modélisables par loi expo.
C'est pourquoi je cherche peut être des choses un peu artificielles.

lourrran

1) tu as créé des données fictives avec un process 'anormal'... fallait pas. Fallait éventuellement mettre un bruit du type y=exp(ax+b+random() ).
2) de façon plus constructive, cherchons une approche.
Si on 'supprime' les données avec y<0, on biaise le résultat, on garde majoritairement les lignes où le bruit était 'positif'.
Si on traite toutes les lignes, en remplaçant les valeurs y=0 par une valeur constante arbitraire, positive, idem, on biaise le résultat.
En dernier recours, j'appliquerais une espèce de moyenne mobile. L'utilisateur sait que les données $y$ représentent une mesure, qui doit toujours être positive, il sait qu'il doit récupérer d'une façon où d'une autre les valeurs négatives. Il sait que le bruit provient d'une variable aléatoire qui a été ajoutée à la vraie mesure. Il fait en sorte de lisser les valeurs. Chaque valeur $a_i$ est remplacée par la moyenne de $a_{i-1}$, $a_i$ et $a_{i+1}$ jusqu'à ce que toutes les valeurs soient positives. On s'assure que la somme de toutes les valeurs reste constante, on ne fait pas ce traitement uniquement 'à la hausse'. 
C'est du bricolage, ça va améliorer les données sur certains intervalles, mais ça va les polluer sur d'autres intervalles. 
Bof. 100 fois bof.

Clairement, dans le cadre de formation pour des étudiants, ok pour leur donner des valeurs produites artificiellement, mais pas avec des y<0 si ils sont sensés trouver une loi exponentielle.

gerard0

Tu trouveras des ajustements exponentiels en physique, dans le cadre de la mise en évidence de "lois de puissances". Peut-être en astronomie, dans la classification des étoiles. Mais je n'ai jamais vu les données brutes.

Une autre idée, en démographie, par exemple l'évolution de la population de la Chine, ou de l'Inde, de 1960 à aujourd’hui.

Cordialement.

el_douwen

merci bcp pour vos réponses !

leon1789

Bonjour
Une idée pour réaliser des jeux de données à partir d'un modèle quelconque ( ici Y = exp(aX+b) ) :
on peut choisir des x_i , des valeurs des paramètres (ici a, b), 
puis on perturbe légèrement les paramètres : des a_i proches de a,  et des b_i proches de a
et on calcule les y_i = exp(a_i x_i+ b_i) correspondants.

(on ne modifie pas les x_i, ni des y_i calculés avec des paramètres fixés)



En l’occurrence ici, cela revient au conseil déjà donné y=exp(ax+b+ε(x)) 

JJ

 Bonjour,

Une possibilité (à vérifier) serait de considérer la fonction $y=a+b\:e^{cx}$ au lieu de $y=e^{\alpha x+\beta}$.

Le résultat sera $c\simeq\alpha$ et $b\simeq e^\beta$. Si la dispersion des $y$ (qui en rend certains négatifs) est faible, $a$ sera petit et pourrait éventuellement être négligeable.

En contre-partie, cette fonction $y=a+b\:e^{cx}$ nécessite une méthode de régression non-linéaire, ce qui est un peu plus compliqué. Les logiciels usuels effectuent un calcul itératif initialisé par des valeurs "supposées" ou "estimées" des paramètres. C'est une difficulté car si les valeurs initiales sont trop mal choisies le calcul itératif peut échouer.

Il existe une méthode (non usuelle) non itérative et ne nécessitant pas de donner de valeurs initiales aux paramètres.

Dans le cas de la fonction $y=a+b\:e^{cx}$ le schéma de calcul est très simple :

Un exemple de calcul est donné ci-dessous. C'est un cas où certaines valeurs de $y$ sont négatives en raison de la dispersion qui a été appliquée.

Cette méthode est , en fait, une régression linéaire mais portant sur les paramètres d'une équation intégrale linéaire dont la fonction à ajuster est solution. Dans le cas présent :

Pour plus d'explication voir : https://fr.scribd.com/doc/14674814/Regressions-et-equations-integrales

leon1789

Bonjour
Bien sûr JJ, ce que je tentais d'expliquer hier à votre ami Dlz, il existe des ajustements sur y = c + exp(a*x+b).
(loin de la loi exponentielle en probabilité. Malheureusement, il n'a rien compris, comme souvent)

Mais il semble bien que le souci de el_douwen est de traiter un ajustement sur le modèle y = exp(a*x+b) :
<< comment, à partir des données .csv, retrouver a et b ? >> via une méthode classique (programme de son enseignement)

leon1789

Sur votre exemple, avec le modèle y = a + b*exp(c*x),
je trouve cette solution en voulant minimiser l'erreur quadratique.


Dans votre méthode, on voit bien que vous avez approché une intégrale par la méthode des trapèzes pour déterminer c (que l'on peut exprimé comme un déterminant d'une matrice 3x3).
Puis méthode des moindres carrés classique pour déterminer a et b.
Au final, votre méthode optimise-t-elle une "fonction de contrôle" ?

JJ

Bonjour leon1789,

Je suis d'accord avec votre commentaire. Vous avez bien décrit ce qui est fait: approximation de l'intégrale puis calcul des coefficients. Bien entendu ce n'est pas une méthode usuelle qui minimiserait l'erreur quadratique moyenne. C'est une méthode non itérative ne nécessitant pas de donner des valeurs initiales estimées à priori. C'est son avantage. Bien entendu, le résultat est plus ou moins approximatif selon la précision du calcul numérique d'intégration. Si le résultat n'est pas suffisamment précis, il peut servir comme excellente donnée initiale pour le calcul itératif classique. On notera que le critère d'ajustement (LMSE, ou LMSRE, ou LMAE, ou etc.) n'était pas spécifié dans la question posée. On peut donc choisir le critère que l'on veut et (pourquoi pas ?) celui utilisé dans ma réponse. Néanmoins il est vrai que LMSE est le plus courant et souvent sous-entendu. Dans ce cas ma réponse aurait du être complétée par la régression non-linéaire usuelle (Comme vous l'avez fait très à propos).

el_douwen

En reprenant l'idée de JJ il me semble que dans le $a+exp(b x+c)$ on pourrait estimer empiriquement le $a$ :

par un truc du style :

mean(df[df$x>M,][,2])

et on cherche la "limite" de cette moyenne glissante lorsque M grandit en s'approchant de max(df[1,])
puisque mathématiquement $a$ est la limite de la fonction si $x\to+\infty$ (j'imagine ici $b<0$).

Je n'ai pas le temps d'essayer, mais j'aurai le temps cet été

JJ

En complément à ma réponse précédente.

Comme attendu, le calcul par régression non-linéaire ( itératif, en partant des valeurs trouvées antérieurement) donne un RLMSE un peu amélioré (Root Least Mean Square Error). Il y a très peu de différence en représentation graphique.

leon1789

Merci pour votre réponse très claire JJ.
Je m'interroge sur le modèle y(x) = a + b*exp(c*x)
En reprenant votre idée d'approximation des intégrales sur les intervalles $[x_{k-1}, x_{k}]$ : 
$\int_{x_{k-1}}^{x_k} y(x) \ dx = (\ x_k - x_{k-1} \ )a + (\ y(x_k) - y(x_{k-1})\ )/c$. 

J'en suis arrivé à cette conclusion : pour $k=2,..,n$, on pose
$\alpha_k = x_k - x_{k-1}$ et
$\beta_k = y_k - y_{k-1}$ et
$\gamma_k = (y_k + y_{k-1}).(x_k-x_{k-1}) / 2$   (aires des trapèzes)

on calcule alors $c = \frac{ det(M) }{ det(N) }$  où $M = \left( \begin{matrix} 
\sum_{k=2}^n \alpha_k^2 & \sum_{k=2}^n \alpha_k \beta_k \\
\sum_{k=2}^n \alpha_k \beta_k & \sum_{k=2}^n \beta_k^2
\end{matrix} \right)$ et $N = \left( \begin{matrix} 
\sum_{k=2}^n \alpha_k^2 & \sum_{k=2}^n \alpha_k \gamma_k \\
\sum_{k=2}^n \alpha_k \beta_k & \sum_{k=2}^n \beta_k \gamma_k
\end{matrix} \right)$

Puis on termine en déterminant a et b par moindres carrés habituels avec la dernière formule que vous avez rappelée : 


Numériquement, on obtient c = 1.804 , puis a = -0.470 et b = 0.006423

leon1789

el_douwen ,
il y a en fait deux stratégies :
déterminer a, puis résoudre usuellement y - a = exp(b*x +c)
ou bien déterminer c, puis résoudre usuellement y = a + b*exp(c*x) 

Personnellement, je préfère la seconde car il n'y a pas d'usage de logarithme (créant des contraintes).

JJ

@ leon1789,

C'est du beau travail.

En fait, il y a beaucoup de variantes qui donnent des résultats voisins. 

Je pense que la méthode que vous proposez devrait donner un RLMSE un peu meilleur que la mienne  car votre approche pour l'intégration est plus élaborée. Néanmoins c'est un peu plus compliqué que le simple calcul des $S_k$. De mon coté, je recherchais vraiment la plus grande simplicité "quoi qu'il en coûte" sur l'ajustement. De toutes façons, si l'on veut définitivement le meilleur ajustement (relativement à un critère bien spécifié), que ce soit ma méthode ou la vôtre ou d'autres variantes, on ne peut pas faire mieux que la régression non-linéaire itérative (Quand elle veut bien converger ! ). 

JJ

Bonjour,

Je crains que l'on se soit éloigné de la question initialement posée par el_douwen. 

Supposant que sa préoccupation est surtout d'avoir une méthode de calcul la plus simple possible, voici une version élémentaire qui ne met même pas en oeuvre de calcul matriciel (du moins en apparence).

EXPLICATIONS :

Il s'agit de l'ajustement de la fonction $y(x)=b\;e^{\,c\; x}$ par la méthode  dont le principe général est décrit dans  https://fr.scribd.com/doc/14674814/Regressions-et-equations-integrales

Dans le cas présent l'équation intégrale est  $y(x)=y_1+\frac{1}{c}\int_{x_1}^{x}y(\chi)d\chi$

Intégration numérique :  $\int_{x_1}^{x_k}y(\chi)d\chi \simeq S_k=S_{k-1}+\frac12 (y_{k-1}+y_k)(x_k-x_{k-1})$

Pour l'ajustement du paramètre $c$ , régression linéaire  :  $y_k-y_1\simeq\frac{1}{c}S_k$

Pour l'ajustement du paramètre $b$ , régression linéaire  :  $y_k\simeq b\;e^{\,c\;x_k}$.

leon1789

Bonjour

@JJ
Je trouve mon approche (qui est totalement inspirée de la vôtre) plus naturelle, plus simple, mais c'est subjectif.
Par ailleurs, comme vous dites, il y a beaucoup de variantes pour déterminer une valeur convenable du facteur c. 
Je suis d'accord que les méthodes itératives vont donner de meilleurs résultats, mais elles sont plus lourdes à implémenter et plus longues en calcul.

leon1789

J'en reviens à la question initiale de el_douwen demandant un ajustement suivant le modèle  y=exp(a*x+b) lorsque les données en y_k pouvent être un peu négatives. 
On reprend toujours la même idée, $\int_{x_{k-1}}^{x_k} y(x) \ dx = (\ y(x_k) - y(x_{k-1})\ )/a$, mais les choses se simplifient :

pour $k=2,..,n$, on pose
$\beta_k = y_k - y_{k-1}$ et
$\gamma_k = (y_k + y_{k-1}).(x_k-x_{k-1}) / 2$   (aires des trapèzes)

Alors on peut prendre $a = \frac{\sum_{k=2}^n \beta_k^2}{ \sum_{k=2} \beta_k .\gamma_k} $, puis $b = \ln \left( \frac{\sum_{k=1}^n e^{a.x_k}.y_k }{\sum_{k=1}^n e^{2a.x_k}} \right)$

Le logarithme existe bien car au moins un $y_k$ est suffisamment grand positif...

En recopiant JJ;
Pour l'ajustement du paramètre 1/a , régression linéaire  :  $\gamma_k ≃ (y_k−y_{k-1}) . 1/a$
Pour l'ajustement du paramètre b, régression linéaire  :  $y_k≃b. e^{c.x_k}$

Numériquement, on obtient a=1.813 et b=-5.103  ( pour faire le lien avec les résultats précédents : exp(b) = 0.006076 )

RLMSE =1.347 , donc légèrement moins efficace que la solution de JJ ci-dessus (qui est très très proche de la solution optimale !!) .

JJ

Je pense que el_douwen devrait être satisfait, que ce soit avec ma réponse ou la vôtre. L'écart d'ajustement est insignifiant. Les courbes sont quasiment superposées.