Ajustement de la courbe de décroissance exponentielle en numpy et scipy

Question

J'ai quelques difficultés à ajuster une courbe à certaines données, mais je n'arrive pas à déterminer où je me trompe.

Dans le passé, je l'ai fait cela avec numpy.linalg.lstsq pour les fonctions exponentielles et scipy.optimize.curve_fit pour les fonctions sigmoïde. Cette fois-ci, j'ai souhaité créer un script qui me permettrait de spécifier différentes fonctions, de déterminer les paramètres et de tester leur adéquation aux données. En faisant cela, j'ai remarqué que Scipy leastsq et Numpy lstsq semblent fournir des réponses différentes pour le même ensemble de données et la même fonction. La fonction est simplement y = e^(l*x) et est contrainte de telle sorte que y=1 à x=0.

La ligne de tendance Excel est en accord avec le résultat Numpy lstsq, mais comme Scipy leastsq peut prendre n'importe quelle fonction, il serait bon de déterminer quel est le problème.

import scipy.optimize as optimize 
import numpy as np 
import matplotlib.pyplot as plt 

## Sampled data 
x = np.array([0, 14, 37, 975, 2013, 2095, 2147]) 
y = np.array([1.0, 0.764317544, 0.647136491, 0.070803763, 0.003630962,  0.001485394,  0.000495131]) 

# function 
fp = lambda p, x: np.exp(p*x) 

# error function 
e = lambda p, x, y: (fp(p, x) - y) 

# using scipy least squares 
l1, s = optimize.leastsq(e, -0.004, args=(x,y)) 
print l1 
# [-0.0132281] 


# using numpy least squares 
l2 = np.linalg.lstsq(np.vstack([x, np.zeros(len(x))]).T,np.log(y))[0][0] 
print l2 
# -0.00313461628963 (same answer as Excel trend line) 

# smooth x for plotting 
x_ = np.arange(0, x[-1], 0.2) 

plt.figure() 
plt.plot(x, y, 'rx', x_, fp(l1, x_), 'b-', x_, fp(l2, x_), 'g-') 
plt.show() 

Edition - informations supplémentaires

Le MWE ci-dessus comprend un petit échantillon de l'ensemble de données. Lors de l'ajustement des données réelles, la courbe scipy.optimize.curve_fit présente un R^2 de 0,82, tandis que la courbe numpy.linalg.lstsq, qui est la même que celle calculée par Excel, a un R^2 de 0,41 .

Answer 1

Vous minimisez différentes fonctions d'erreur.

Lorsque vous utilisez numpy.linalg.lstsq, la fonction d'erreur étant réduite au minimum est

np.sum((np.log(y) - p * x)**2) 

tout scipy.optimize.leastsq minimise la fonction

np.sum((y - np.exp(p * x))**2) 

Le premier cas nécessite une dépendance linéaire entre les variables dépendantes et indépendantes, mais la La solution est connue par analogie, tandis que la seconde peut gérer n'importe quelle dépendance, mais repose sur une méthode itérative.

Sur une note séparée, je ne peux pas le tester en ce moment, mais lorsque vous utilisez numpy.linalg.lstsq, je vous ne pas besoin de vstack une rangée de zéros, les travaux suivants ainsi:

l2 = np.linalg.lstsq(x[:, None], np.log(y))[0][0] 

Answer 2

Pour exposer un peu sur le point de Jaime, toute transformation non linéaire des données conduira à une fonction d'erreur différente et donc à des solutions différentes. Cela conduira à des intervalles de confiance différents pour les paramètres d'ajustement. Vous avez donc trois critères possibles à utiliser pour prendre une décision: quelle erreur voulez-vous minimiser, quels paramètres vous voulez plus de confiance, et enfin, si vous utilisez le raccord pour prédire une valeur, quelle méthode donne moins d'erreur dans l'intéressant valeur prédite. En jouant autour d'un bit analytiquement et dans Excel, on suggère que différents types de bruit dans les données (par exemple si la fonction de bruit ajuste l'amplitude, affecte la constante de temps ou est additif) conduit à différents choix de solution. J'ajouterai également que bien que cette astuce "fonctionne" pour la décroissance exponentielle à 0, elle ne peut pas être utilisée dans le cas plus général (et commun) des exponentielles amorties (ascendantes ou descendantes) à des valeurs qui ne peuvent pas être supposé être 0.