Schwierigkeiten mit curve_fit

[meitnerium]

Schönen Abend euch allen!

Ich arbeite daran eine Funktion auf Daten aus einer Titrationskurve zu fitten, zum Glück weiß ich, dass die Funktion in etwa eine logistische Funktion sein sollte. Leider passen die ausgegebenen Parameter nicht wirklich zu den Datenpunkten (s. Abbildung).

Code: Alles auswählen

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

xdata = numpy.array([0, 5,   7,   9.9, 10, 10.1, 13,   15,   20])
ydata = numpy.array([1, 1.5, 1.7, 4,   7,  10,   12.3, 12.5, 13])
xaxis = numpy.array([0, 20])

def func(x, L, k, x0):
    return L/(1 + numpy.exp(-k*(x-x0)))
popt, pcov = curve_fit(func, xdata, ydata, bounds=([12., 1., 7.5], [14., numpy.inf, 8.5]))

plt.plot(xaxis, func(xaxis, *popt))
plt.plot(xdata, ydata, "bo")
print(popt)
plt.show()

Die optimierten Parameter für die Funktion, die sich dabei ergeben, sind dann:

Code: Alles auswählen

[12.   1.   8.5]

Der Graph der Funktion im Vergleich mit den Datenpunkten sieht dann so aus:


Das ist natürlich nicht das, was ich wollte... Hat jemand eine Idee, was ich falsch gemacht habe? Vielen Dank schon einmal!

[narpfel]

Moin,

dein `xaxis` hat nur zwei Punkte (bei x = 0 und x = 20), da kann also nur eine Gerade bei rauskommen. `np.linspace` erstellt dir ein Array aus gleichmäßig verteilten Werten in einem Intervall.

[meitnerium]

Vielen Dank, was für ein blöder Fehler...

Ich habe das behoben, leider ist die gefittete Funktion dennoch viel zu flach:

Dazu hast du nicht zufällig auch eine Idee?

[narpfel]

Das sieht so aus, als wenn deine Messwerte gegen ein echtes logistisches Wachstum nach oben verschoben sind. Eventuell hilft ein `+ y_0` in der Modellfunktion (je nach dem, ob das physikalisch Sinn macht), ansonsten müsstest du dir Gedanken darüber machen, was die richtige Modellfunktion ist.

[webbygirl21]

Hat das dann noch geklappt mit dem y_0? Falls es nicht die richtige Modellfunktion wäre, wie würde man denn da jetzt weiter vorgehen?

[meitnerium]

Vielen Dank nochmal für die Hinweise!
Ich habe leider erst jetzt nochmal Zeit gefunden mich damit auseinanderzusetzen und hmmm... Die Modellfunktion ist durchaus geeignet; ich hatte das nochmal nachgeschaut und man sieht es auch an der grünen Funktion (mit manuell ausprobierten Parametern). Leider liefert curve_fit dennoch nur einen sehr eckigen Verlauf, egal ob die Bestimmung ganz frei oder mit relativ engen Grenzen stattfindet (s. Grafik). Der aktuelle Stand meiner Versuche ist also:

Code: Alles auswählen

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

xdata = numpy.array([0, 5,   7,   9.9, 10, 10.1, 13,   15,   20])
ydata = numpy.array([1, 1.5, 1.7, 4,   7,  10,   12.3, 12.5, 13])
xaxis = numpy.linspace(0, 20)

def func(x, L, a, k, x0, y0):
    return (L/(a + numpy.exp(-k*(x-x0))))+y0
popt, pcov = curve_fit(func, xdata, ydata)
popt2, pcov2 = curve_fit(func, xdata, ydata, bounds=([0, 0, 0, 9, 0.], [30, 10, 10, 11, 2.]))

def func2(x):
        return (10/(0.8+numpy.exp(-3*(x-10))))+1
        
plt.plot(xaxis, func(xaxis, *popt), 'r')
plt.plot(xaxis, func(xaxis, *popt2), 'orange')
plt.plot(xaxis, func2(xaxis), 'g')

plt.plot(xdata, ydata, "bo")
print(popt)
plt.show()

und gibt mir diesen Plot:


Auf die Idee, die Parameter einzugrenzen, kam ich vor allem wegen dieser Warnung

Code: Alles auswählen

runfile('C:/Users/TS/Desktop/untitled1.py', wdir='C:/Users/TS/Desktop')
[ 2.96832762  0.26502925 11.96250756  9.88899366  1.39999999]
C:\Users\TS\Desktop\untitled1.py:10: RuntimeWarning: overflow encountered in exp
  return (L/(a + numpy.exp(-k*(x-x0))))+y0

das hat aber wie man sieht nichts gebracht, die Funktion ist dennoch schrecklich eckig (wie auch immer das sein kann...)

Ich gebe dann jetzt wohl mal auf...



Und @webbygirl21: Ich bin (wie man leider merkt) kein Numeriker, aber wenn man gar keine Ahnung hat, welche Funktion das modelliert, dann kann man immer zu einer Polynom- oder Splineinterpolation greifen, denn jede Funktion kann ja durch ein Polynom hinreichend hohen Grades angenähert werden - das braucht aber dann seeeehr viel Rechenleistung. Das habe ich übrigens auch ausprobiert:

Code: Alles auswählen

import numpy
import matplotlib.pyplot as plt
from scipy.interpolate import Rbf, InterpolatedUnivariateSpline
x = numpy.array([0, 5, 7, 9.9, 10, 10.1, 13, 15, 20])
y = numpy.array([1, 1.5, 1.7, 4, 7, 10, 12.3, 12.5, 13])
xaxis = numpy.linspace(0, 20, num=101, endpoint=True)

#k gibt quasi an, welcher Grad an Glattheit gefordert ist: für k=1
#folgt einfach ein Linienzug
f1=InterpolatedUnivariateSpline(x, y, k=1)
f2=InterpolatedUnivariateSpline(x, y, k=2)
f3=InterpolatedUnivariateSpline(x, y, k=3)

plt.plot(x,y, 'o')
plt.plot(xaxis, f1(xaxis), 'b')
plt.plot(xaxis, f2(xaxis), 'r')
plt.plot(xaxis, f3(xaxis),'y')

plt.show()

Man bekommt dann diese wunderschönen Ergebnisse:


Außerdem kann man natürlich die "üblichen Verdächtigen" durchprobieren oder jemanden mit genug Erfahrung (bzw das Internet) um Hilfe bitten; manche Menschen schauen sich die Datenpunkte an und wissen, welche Funktion sie gut beschreibt...

[narpfel]

@meitnerium: Das ist eckig, weil dein `linspace` zu wenige Punkte hat. Ändert allerdings nichts daran, dass die Werte nicht wirklich gut auf eine logistische Kurve passen. `curve_fit` findet da schon mehr oder weniger die besten Parameter. Vielleicht hilft es, noch die Fehlerbalken der Messwerte anzugeben.

manche Menschen schauen sich die Datenpunkte an und wissen, welche Funktion sie gut beschreibt...

Eher andersrum. Die Funktion kommt normalerweise nicht aus den Daten, sondern aus der Physik/Chemie/was auch immer hinter den Daten. In anderen Worten: Man sollte die Funktion kennen, bevor man sich daran macht, die Daten auszuwerten. Das sieht man sehr gut an den Splines, die du gefunden hast: Sie lassen sich perfekt an die Daten anfitten, aber man lernt daraus nichts über das vermessene System, weil die Funktionen, die da rauskommen, nichts mit der Realität zu tun haben.

Ich bin kein Chemiker (und ich kenne den Versuch nicht) und kann deswegen nicht beurteilen, ob eine logistische Kurve den Zusammenhang beschreibt, ob das eine Näherung ist und es nicht richtig passt, weil sich das System nicht wie idealisiert verhält, oder ob die logistische Kurve (wie die Splines) gar keinen Zusammenhang mit dem System hat. Vielleicht ist das Ergebnis also genau so, wie man das erwarten würde?

[meitnerium]

Vielen vielen Dank narpfel, damit ist mein Problem jetzt endlich gelöst! Genau so hatte ich mir das vorgestellt, wenn ich die echten Daten auswerte, werde ich auch das mit den Fehlerbalken machen.



Falls es noch wen interessiert, hier mein endgültiger Code (ich habe nur noch einen linearen Term hinzugefügt, damit der Anstieg neben dem Sprung abgebildet werden kann):

Code: Alles auswählen

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

xdata = numpy.array([0, 5,   7,   9.9, 10, 10.1, 13,   15,   20])
ydata = numpy.array([1, 1.5, 1.7, 4,   7,  10,   12.3, 12.5, 13])
xaxis = numpy.linspace(0, 20, num=500)

def func(x, L, a, k, x0, y0, t):
    return (L/(a + numpy.exp(-k*(x-x0))))+y0+t*x
popt2, pcov2 = curve_fit(func, xdata, ydata, bounds=([0, 0, 0, 9, 0., 0.], [30, 10, 10, 11, 2., numpy.inf]))
      
plt.plot(xaxis, func(xaxis, *popt2), 'b')

plt.plot(xdata, ydata, "bx")
print(popt2)
plt.show()