y = A*x^B und scipy.curve_fit

mit matplotlib, NumPy, pandas, SciPy, SymPy und weiteren mathematischen Programmbibliotheken.
Antworten
zweihorn
User
Beiträge: 25
Registriert: Donnerstag 11. Mai 2017, 17:09

Montag 27. November 2017, 16:43

Hallo!

Ich versuche gerade Daten zu fitten, welche generell der Form y = A*x^B folgen.

Ich habe bisher drei Varianten getestet:

Variante 1: Linearisieren und mit nump.polyfit fitten.
Variante 2: stumpf scipy.curve_fit anwenden (Startwerte geschätzt, ich gebe ja die analytische form oben im code vor) <- erscheint mir am elegantesten.
Variante 3: Adaption des Beispiels von hier http://scipy-cookbook.readthedocs.io/it ... gData.html (ganz unten).

Code: Alles auswählen

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

powerlaw = lambda x, amp, index: amp * (x**index)

num_points = 20

xdata = np.linspace(1.1, 10.1, num_points)
ydata = powerlaw(xdata, 10.0, -2.0)     # simulated perfect data
yerr = 0.1 * ydata                      # simulated errors (10%)

ydata += np.random.randn(num_points) * yerr# 

########### Variante 1: ###########

A = np.polyfit(np.log10(xdata), np.log10(ydata), 1)

amp1 = 10**A[1]

yn1=amp1*xdata**A[0]

plt.plot(xdata,ydata,'ro')
plt.plot(xdata,yn1,'g-')

########### Variante 2: ###########

B = optimize.curve_fit(lambda t,a,b: a*t**b,  xdata,  ydata,  p0=(9, -1))
yn2=B[0][0]*xdata**B[0][1]

plt.plot(xdata,yn2,'b--')

########### Variante 3: ###########

logx = np.log10(xdata)
logy = np.log10(ydata)
logyerr = yerr / ydata

# define our (line) fitting function
fitfunc = lambda p, x: p[0] + p[1] * x
errfunc = lambda p, x, y, err: (y - fitfunc(p, x)) / err

pinit = [1.0, -1.0]
out = optimize.leastsq(errfunc, pinit,
                       args=(logx, logy, logyerr), full_output=1)

pfinal = out[0]

index = pfinal[1]
amp = 10.0**pfinal[0]

plt.plot(xdata, powerlaw(xdata, amp, index),'r-')

plt.show()
Betrachtet man die Graphen, sind sich Variante 1 und 3 recht einig (solange yerr != 0). Variante 2 weicht etwas davon ab. Setze ich yerr = 0, dann spielt Variante 3 komplett verrückt.

Daher habe ich ein paar Fragen:

1. Weshalb liefert mir curve_fit einen subjektiv schlechteren fit? (Bei echten Daten sieht man es deutlicher als in diesem Testfall)
2. Weshalb steigt Variante 3 aus, sobald der yerr = 0 ist?
3. Was wäre denn der ideale Weg des fittings?

Danke!
Benutzeravatar
Sr4l
User
Beiträge: 1091
Registriert: Donnerstag 28. Dezember 2006, 20:02
Wohnort: Kassel
Kontaktdaten:

Montag 27. November 2017, 19:36

Ich kann dir nicht ganz Folgen. Mit deinen original Script kann ich qualitativ keinen unterschied der Verfahren sehen. Natürlich kommt es je "Zufall" zu leichten Variationen, aber es sind ja auch unterschiedliche Methoden.

Das es yerr = 0 ist bei Variante bestimmt ein Problem, weil Teilen durch Null nie eine gute Idee ist ('errfunc = lambda p, x, y, err: (y - fitfunc(p, x)) / err').
narpfel
User
Beiträge: 183
Registriert: Freitag 20. Oktober 2017, 16:10

Montag 27. November 2017, 20:08

@zweihorn: Was das beste Vorgehen ist, hängt IMHO von den Daten ab. Beim Logarithmieren hast du das Problem, dass deine Fehlerbalken mit logarithmiert werden. Dadurch werden beim Fitten einige Werte stärker gewichtet als andere. Eventuell macht es also Sinn, die Standardabweichungen der Werte an die Fitfunktion zu übergeben.

Zu Frage 3: Du kannst ja mal ein paar Testdaten generieren und jeweils den mittleren relativen Fehler zu den Originaldaten (also den nicht verrauschten) berechnen. Würde mich auch interessieren, was dabei rauskommt.
zweihorn
User
Beiträge: 25
Registriert: Donnerstag 11. Mai 2017, 17:09

Dienstag 28. November 2017, 11:22

Sr4l hat geschrieben:Ich kann dir nicht ganz Folgen. Mit deinen original Script kann ich qualitativ keinen unterschied der Verfahren sehen. Natürlich kommt es je "Zufall" zu leichten Variationen, aber es sind ja auch unterschiedliche Methoden.

Das es yerr = 0 ist bei Variante bestimmt ein Problem, weil Teilen durch Null nie eine gute Idee ist ('errfunc = lambda p, x, y, err: (y - fitfunc(p, x)) / err').
Ja, das Beispiel ist schlecht, ich habe den Code verändert. Die Daten werden nun von zwei sich überlappenden Potenzfunktionen dargestellt. Das ist nur um das Problem zu verdeutlichen, in der Realität kommt der unterschied aus Messungenauigkeiten!

Code: Alles auswählen

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

powerlaw = lambda x, amp, index: amp * (x**index)

num_points = 20

xdataa = np.linspace(1.1, 10.1, num_points)
xdatab = np.linspace(5.1, 15.1, num_points)
ydataa = powerlaw(xdataa, 124, -0.7999)     # simulated perfect data
ydatab = powerlaw(xdatab, 160, -1.2)
xdata = np.concatenate((xdataa,xdatab))
ydata = np.concatenate((ydataa,ydatab))

data = np.array([xdata,ydata]).T

data = data[data[:,0].argsort()]

xdata = data[:,0]
ydata = data[:,1]

yerr = 0.05 * ydata                      # simulated errors (10%)

ydata += np.random.randn(len(ydata)) * yerr

# Variante 1:

A = np.polyfit(np.log10(xdata), np.log10(ydata), 1)

amp1 = 10**A[1]

yn1=amp1*xdata**A[0]

plt.plot(xdata,ydata,'ro')
plt.plot(xdata,yn1,'g-')

# Variante 2:

B = optimize.curve_fit(lambda t,a,b: a*t**b,  xdata,  ydata,  p0=(9, -1))
yn2=B[0][0]*xdata**B[0][1]

plt.plot(xdata,yn2,'b--')


# Variante 3:

logx = np.log10(xdata)
logy = np.log10(ydata)
logyerr = yerr / ydata

# define our (line) fitting function
fitfunc = lambda p, x: p[0] + p[1] * x
errfunc = lambda p, x, y, err: (y - fitfunc(p, x)) / err

pinit = [10, -2.0]
out = optimize.leastsq(errfunc, pinit,
                       args=(logx, logy, logyerr), full_output=1)

pfinal = out[0]

index = pfinal[1]
amp = 10.0**pfinal[0]

plt.plot(xdata, powerlaw(xdata, amp, index),'r-')

plt.show()
Die curve_fit Variante läuft sozusagen zwischen den beiden Potenzfunktionen entlang. Die beiden anderen Varianten schaffen es, bei größeren x-Werte den untern Verlauf zu treffen. Daher meine Vermutung, dass curve_fit "schlechter" passt. Bei niedrigen x-Werten ist es hingegen anders herum. Müssten diese beiden Methoden nicht ziemlich ähnlich verlaufen (von den Rundungsfehlern und etwas Zufall mal abgesehen)?


Den "/ err" Teil in Variante 3 habe ich gestern total übersehen... Es war wohl an der Zeit eine Pause zu machen :mrgreen: Blöde Frage: Weshalb ist das überhaupt auf den Fehler normiert? Also was ist der genaue Grund für das Teilen durch "err"?

narpfel hat geschrieben:@zweihorn: Was das beste Vorgehen ist, hängt IMHO von den Daten ab. Beim Logarithmieren hast du das Problem, dass deine Fehlerbalken mit logarithmiert werden. Dadurch werden beim Fitten einige Werte stärker gewichtet als andere. Eventuell macht es also Sinn, die Standardabweichungen der Werte an die Fitfunktion zu übergeben.
Betrifft das nicht nur die Wahl des Modells? Also welche Funktion ich zum fitten nutze. Bei gleicher Funktion müssten die doch auf einen ziemlich ähnlichen Verlauf kommen oder nicht?

Was meinst du denn mit "Fehlerbalken mit logarithmiert werden. Dadurch werden beim Fitten einige Werte stärker gewichtet als andere", kannst du das etwas genauer erläutern?


Kurze Zwischenfrage: curve_fit gibt mir ja pcov aus. Daraus müsste ich doch R^2 berechnen können, nach der Formel R^2 = (cov(x,y)*cov(x,y))/(var(x)*var(y)), oder? (Ja ich weiß, dass R^2 nur bei einem linearen Modell anwendbar ist, hat jetzt nur indirekt mit dem Problem von oben zu tun).

Danke für eure Hilfe!
Benutzeravatar
Sr4l
User
Beiträge: 1091
Registriert: Donnerstag 28. Dezember 2006, 20:02
Wohnort: Kassel
Kontaktdaten:

Mittwoch 29. November 2017, 23:59

zweihorn hat geschrieben:Die curve_fit Variante läuft sozusagen zwischen den beiden Potenzfunktionen entlang. Die beiden anderen Varianten schaffen es, bei größeren x-Werte den untern Verlauf zu treffen. Daher meine Vermutung, dass curve_fit "schlechter" passt. Bei niedrigen x-Werten ist es hingegen anders herum. Müssten diese beiden Methoden nicht ziemlich ähnlich verlaufen (von den Rundungsfehlern und etwas Zufall mal abgesehen)?
Wenn du was haben willst was Rauschen reduziert dann willst du z.B Tiefpassfiltern.

Wenn du was fittest dann willst du die Daten möglichst einfach beschreiben. Die einfachste Beschreibung ist eine Konstante, dann vll eine Gerade usw. Du kannst also nur so gut an deinem Geraden liegen wie die Funktion, die du gewählt hast, es ermöglicht. Beim Rest geht es den Verfahren immer nur darum die Kosten zu minimieren. Das heißt das die Funktion meist da am besten passen, wo die meisten Datenpunkte liegen oder Ausreißer Datenpunkte zu besonders hohen Kosten führen.

Das ist also alles sehr stark abhängig von dem was du erreichen willst.

Hier mal ein paar Beispiele zu deinen xdata, ydata aus dem letzen Beispiel:

Code: Alles auswählen

from numpy.polynomial import Polynomial
p = Polynomial.fit(xdata, ydata, 1)
plt.plot(*p.linspace())

p = Polynomial.fit(xdata, ydata, 2)
plt.plot(*p.linspace())

p = Polynomial.fit(xdata, ydata, 10)
plt.plot(*p.linspace())
zweihorn
User
Beiträge: 25
Registriert: Donnerstag 11. Mai 2017, 17:09

Sonntag 3. Dezember 2017, 12:41

Sorry für die späte Rückmeldung.

Ich kenne den Unterschied zwischen Filtern und Fitten.

Mir geht es um den Unterschied, zwischen der Polyfit und der curve_fit Methode. Wie ich mittlerweile herausgefunden habe, scheint es am von narpfel geschilderten Problem mit dem Logarithmieren der Fehler zu liegen, bzw. der Übergewichtung kleiner Werte durch das Logarithmieren.

Bei Polyfit wird der Datensatz erst linearisiert und dann gefittet. Hierbei werden anscheinend die kleinen Werte zu stark, bzw. stärker gewichtet. Mein Verständnisproblem ist nun: Warum passt der Fit mit diesem vermeintlichen Fehler besser als der nicht-lineare Fit ohne diesen Fehler?


Meine Zusatzfrage mit pcov hat sich erledigt, das sind nur die Werte der Parameter und nicht des ganzen Fits. Ich habe einen anderen Weg gefunden.

LG
zweihorn
narpfel
User
Beiträge: 183
Registriert: Freitag 20. Oktober 2017, 16:10

Montag 4. Dezember 2017, 10:17

@zweihorn: `polyfit` ist außerdem noch exakt (liefert immer das beste Polynom), `curve_fit` dagegen für *alle* Funktionen benutzbar.

Dass `curve_fit` tendenziell schlechtere Fits liefert als die `polyfit`-Variante könnte daran liegen, dass deine Fehlerbalken proportional zum Funktionswert sind und damit beim Logarithmieren alle gleich groß werden. Für vergleichbarere Ergebnisse könntest du die y-Fehler an `curve_fit` übergeben und gucken, ob der Fit besser wird.
Antworten