2015-11-17 5 views
6

Ich versuche einige Daten an eine Potenzgesetzfunktion mit exponentiellem Abschneiden anzupassen. Ich erzeuge einige Daten mit numpy und ich versuche, diese Daten mit scipy.optimization anzupassen. Hier ist mein Code:Falsche Anpassung mit scipy curve_fit

import numpy as np 
from scipy.optimize import curve_fit 

def func(x, A, B, alpha): 
    return A * x**alpha * np.exp(B * x) 

xdata = np.linspace(1, 10**8, 1000) 
ydata = func(xdata, 0.004, -2*10**-8, -0.75) 
popt, pcov = curve_fit(func, xdata, ydata) 
print popt 

Das Ergebnis erhalte ich ist: [1, 1, 1], die nicht mit den Daten entspricht. ¿Mache ich etwas falsch?

+1

Meine Methode bricht, weil Ihr letzten Datenpunkt mit (bei 655.642.210) hat einen Wert von 0. Wenn Sie das Protokoll davon nehmen, erhalten Sie ein NaN.Ich habe die Anpassung mit meiner Methode berechnet, die diesen Punkt ausschließt und etwas erhält, das vernünftig aussieht. A = 0,00326, Alpha = -0,767, B = -1,88e-8 –

+0

Ja, Sie haben recht! Ich habe meinen Fehler bemerkt und meine zweite Sekunde gelöscht. Ich danke dir sehr – ivangtorre

Antwort

4

Während xnx Sie auf die Antwort gab, warum curve_fit hier scheiterte ich dachte, dass ich d schlagen eine andere Art vor, sich dem Problem der Anpassung Ihrer funktionalen Form zu nähern, die nicht auf einem Gradientenabstieg beruht (und daher eine vernünftige anfängliche Schätzung).

Beachten Sie, dass, wenn Sie nehmen das Protokoll der Funktion, die Sie passend sind Sie das Formular

\log f = \log A + \alpha \log x + B x

erhalten, die linear in jedem der unbekannten Parameter (A einzuloggen, alpha, B)

Wir deshalb kann die Maschinerie der linearen Algebra verwenden, um dies zu lösen, indem die Gleichung in der Form einer Matrix zu schreiben als

log y = M p

wobei log y ein Spaltenvektor des Protokoll Ihrer YDATA Punkten Ist p ein Spaltenvektor der unbekannten Parameter und M die Matrix [[1], [log x], [x]]

Oder

enter image description here

explizit Der am besten passende Parametervektor kann dann effizient unter Verwendung np.linalg.lstsq

Ihr Beispiel finden Problem im Code könnte dann geschrieben werden als

import numpy as np 

def func(x, A, B, alpha): 
    return A * x**alpha * np.exp(B * x) 

A_true = 0.004 
alpha_true = -0.75 
B_true = -2*10**-8 

xdata = np.linspace(1, 10**8, 1000) 
ydata = func(xdata, A_true, B_true, alpha_true) 

M = np.vstack([np.ones(len(xdata)), np.log(xdata), xdata]).T 

logA, alpha, B = np.linalg.lstsq(M, np.log(ydata))[0] 

print "A =", np.exp(logA) 
print "alpha =", alpha 
print "B =", B 

, die die Anfangsparameter schön erholt:

A = 0.00400000003736 
alpha = -0.750000000928 
B = -1.9999999934e-08 

Beachten Sie auch, dass diese Methode um 20x schneller als curve_fit für das Problem bei der Hand

In [8]: %timeit np.linalg.lstsq(np.vstack([np.ones(len(xdata)), np.log(xdata), xdata]).T, np.log(ydata)) 
10000 loops, best of 3: 169 µs per loop 


In [2]: %timeit curve_fit(func, xdata, ydata, [0.01, -5e-7, -0.4]) 
100 loops, best of 3: 4.44 ms per loop 
2

Anscheinend ist Ihre anfängliche Schätzung (die standardmäßig [1,1,1] ist, da Sie nicht eine geben - siehe the docs) ist zu weit von den tatsächlichen Parametern, damit der Algorithmus konvergieren kann. Das Hauptproblem ist wahrscheinlich mit B, die, wenn positiv, wird Ihre exponentielle Funktion zu sehr großen Werten für Ihre bereitgestellten xdata senden.

Probieren Sie etwas ein wenig näher an den tatsächlichen Parametern bereitstellt und es funktioniert:

p0 = 0.01, -5e-7, -0.4 # Initial guess for the parameters 
popt, pcov = curve_fit(func, xdata, ydata, p0) 
print popt 

Ausgang:

[ 4.00000000e-03 -2.00000000e-08 -7.50000000e-01]