Exponentielle Abklingkurve, die in numpy und scipy passt

Question

Ich habe ein paar Probleme mit der Anpassung einer Kurve an einige Daten, aber kann nicht herausfinden, wo ich falsch liege.

In der Vergangenheit habe ich diese für Exponentialfunktionen mit numpy.linalg.lstsq getan und scipy.optimize.curve_fit für Sigmoidfunktionen. Dieses Mal wollte ich ein Skript erstellen, mit dem ich verschiedene Funktionen spezifizieren, Parameter bestimmen und ihre Anpassung an die Daten testen kann. Dabei habe ich festgestellt, dass Scipy leastsq und Numpy lstsq unterschiedliche Antworten für den gleichen Datensatz und die gleiche Funktion zu liefern scheinen. Die Funktion ist einfach y = e^(l*x) und ist so eingeschränkt, dass y=1 bei x=0.

Excel-Trendlinie stimmt mit dem Numpy lstsq Ergebnis überein, aber wie Scipy leastsq ist in der Lage, jede Funktion zu übernehmen, wäre es gut zu erarbeiten, was das Problem ist.

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() 

Bearbeiten - zusätzliche Informationen

Die MWE über eine kleine Probe des Datensatzes enthält. Bei der Anpassung der tatsächlichen Daten zeigt die scipy.optimize.curve_fit Kurve ein R^2 von 0,82, während die numpy.linalg.lstsq Kurve, die die gleiche ist wie die von Excel berechnet, eine R^2 von 0,41 aufweist .

Answer 1

Sie minimieren verschiedene Fehlerfunktionen.

Wenn Sie numpy.linalg.lstsq verwenden, wird die Fehlerfunktion minimiert wird, ist

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

während scipy.optimize.leastsq die Funktion minimiert

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

Der erste Fall eine lineare Abhängigkeit zwischen den abhängigen und unabhängigen Variablen erfordert, aber die Lösung ist analytisch bekannt, während die zweite jede Abhängigkeit bewältigen kann, aber auf einer iterativen Methode beruht.

Auf einem separaten Notiz, Ich kann es nicht richtig testen jetzt, aber wenn numpy.linalg.lstsq verwenden, ich Sie nicht brauchen, um eine Reihe von Nullen, die folgenden Werke als auch auf vstack:

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

Answer 2

zu Wenn man ein wenig auf Jaimes Punkt eingeht, führt jede nichtlineare Transformation der Daten zu einer anderen Fehlerfunktion und damit zu anderen Lösungen. Dies führt zu unterschiedlichen Konfidenzintervallen für die Anpassungsparameter. Sie haben also drei mögliche Kriterien, um eine Entscheidung zu treffen: Welchen Fehler möchten Sie minimieren, welchen Parametern möchten Sie mehr Vertrauen schenken, und schließlich, wenn Sie die Anpassung verwenden, um einen Wert vorherzusagen, welche Methode ergibt weniger Fehler in der interessanten vorhergesagter Wert. Etwas analytisch und in Excel herumzuspielen legt nahe, dass verschiedene Arten von Rauschen in den Daten (z. B. wenn die Rauschfunktion die Amplitude skaliert, die Zeitkonstante beeinflusst oder additiv ist) zu unterschiedlichen Lösungsoptionen führt.

Ich füge auch hinzu, dass, während dieser Trick für exponentiellen Zerfall zu 0 "arbeitet", kann es nicht im allgemeineren (und gemeinsamen) Fall von gedämpften Exponentialen (steigend oder fallend) zu Werten verwendet werden, die nicht sein können angenommen, 0 zu sein.