Ich versuche, eine Verteilung mit scipy's curve_fit anzupassen. Ich habe versucht, eine Einkomponenten-Exponentialfunktion zu verwenden, die zu einer fast geraden Linie führte (siehe Abbildung). Ich habe auch eine exponentielle Zwei-Komponenten-Anpassung versucht, die gut zu funktionieren schien. Zwei Komponenten bedeutet nur, dass ein Teil der Gleichung mit verschiedenen Eingabeparametern wiederholt wird. Wie auch immer, hier ist die eine Komponente Fit-Funktion:scipy curve_fit seltsames Ergebnis
def Exponential(Z,w0,z0,Z0):
z = Z - Z0
termB = (newsigma**2 + z*z0)/(numpy.sqrt(2.0)*newsigma*z0)
termA = (newsigma**2 - z*z0)/(numpy.sqrt(2.0)*newsigma*z0)
return w0/2.0 * numpy.exp(-(z**2/(2.0*newsigma**2))) * (numpy.exp(termA**2)*erfc(termA) + numpy.exp(termB**2)*erfc(termB))
und die Armatur erfolgt mit
fitexp = curve_fit(Exponential,newx,y2)
Dann habe ich versucht, etwas, nur um es auszuprobieren. Ich habe zwei Parameter der Zweikomponentenanpassung genommen, aber sie nicht in der Berechnung verwendet.
def ExponentialNew(Z,w0,z0,w1,z1,Z0):
z = Z - Z0
termB = (newsigma**2 + z*z0)/(numpy.sqrt(2.0)*newsigma*z0)
termA = (newsigma**2 - z*z0)/(numpy.sqrt(2.0)*newsigma*z0)
return w0/2.0 * numpy.exp(-(z**2/(2.0*newsigma**2))) * (numpy.exp(termA**2)*erfc(termA) + numpy.exp(termB**2)*erfc(termB))
Und plötzlich funktioniert das.
Nun, meine quation ist. WARUM? Wie Sie sehen können, gibt es absolut keinen Unterschied in der Berechnung der Anpassung. Es erhält nur zwei zusätzliche Variablen, die nicht verwendet werden. Sollte dies nicht das gleiche Ergebnis erzielen?
@Andras_Deak Ein aktuelles Beispiel:
from scipy.special import erfc
import numpy
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
#setup data
x = [-58.,-54.,-50.,-46.,-42.,-38.,-34.,-30.,-26.,-22.,-18.,-14.,-10.,-6.,-2.,2.,6.,10.,14.,18.,22.,26.,30.,34.,38.,42.,46.,50.,54.,58.]
y = [23.06763817, 16.89802085, 17.83258379, 16.63446237, 13.81878965, 12.97965839, 14.30451789, 16.98288216, 22.26811491, 28.56756908, 33.06990344, 38.59842098, 54.19860393, 86.37381604, 137.47253315, 199.49724512, 238.66047662, 219.89405445, 160.68820199, 103.88901303, 65.92405727, 43.84596266, 31.5395342, 25.9610156, 22.71683709, 18.06740651, 13.85362374, 11.12867065, 10.36502799, 11.31855619]
y_err = [17.9823065, 4.13684885, 1.66490726, 2.4109372, 2.93359141, 1.9701747, 3.19214881, 3.65593012, 2.89089074, 3.58922121, 4.25505348, 4.72728874, 6.77736567, 11.3888196, 21.87771722, 39.0087495, 56.6910311, 51.7592369, 26.39750958, 10.62678862, 7.85893395, 8.11741621, 7.91731416, 7.07739132, 5.41818744, 6.11286843, 8.27070757, 7.85323065, 4.26885499, 0.9047867]
#function to fit
def Exponential2(Z, w0, z0, w1, z1, Z0):
z = Z - Z0
s = 3.98098937586
a = z**2/(2.0*s**2)
b = (s**2 + z*z0)/(numpy.sqrt(2.0)*s*z0)
c = (s**2 - z*z0)/(numpy.sqrt(2.0)*s*z0)
d = (s**2 + z*z1)/(numpy.sqrt(2.0)*s*z1)
e = (s**2 - z*z1)/(numpy.sqrt(2.0)*s*z1)
return w0/2.0 * numpy.exp(-a) * (numpy.exp(c**2)*erfc(c) + numpy.exp(b**2)*erfc(b)) + w1/2.0 * numpy.exp(-a) * (numpy.exp(e**2)*erfc(e) + numpy.exp(d**2)*erfc(d))
#derive and set initial guess
ymaxpos = x[numpy.where(y==numpy.max(y))[0]]
p0_2 = [numpy.max(y),5,numpy.max(y)/2.0,20,ymaxpos]
#fit
fitexp2 = curve_fit(Exponential2,x,y,p0=p0_2,sigma=y_err)
#get results
w0err = numpy.sqrt(numpy.diag(fitexp2[1]))[0]
z0err = numpy.sqrt(numpy.diag(fitexp2[1]))[1]
w1err = numpy.sqrt(numpy.diag(fitexp2[1]))[2]
z1err = numpy.sqrt(numpy.diag(fitexp2[1]))[3]
w0 = fitexp2[0][0]
z0 = fitexp2[0][1]
w1 = fitexp2[0][2]
z1 = fitexp2[0][3]
Z0 = fitexp2[0][4]
#new x array for smoother curve
smoothx = numpy.arange(-58,59,0.1)
y2 = Exponential2(smoothx,w0,z0,w1,z1,Z0)
print 'Exponential 2: w0: '+str(w0.round(3))+' +/- '+str(w0err.round(3))+' \t z0: '+str(z0.round(3))+' +/- '+str(z0err.round(3))+' \t w1: '+str(w1.round(3))+' +/- '+str(w1err.round(3))+' \t\t z1: '+str(z1.round(3))+' +/- '+str(z1err.round(3))
#plot
fig = plt.figure()
ax = fig.add_subplot(111)
ax.errorbar(x,y,y_err,fmt='o',markersize=2,label='data')
ax.plot(smoothx,y2,label='fit',color='red')
ax.grid()
ax.legend()
plt.show()
Wie Sie sehen können, ist die Handlung gut aussehen, aber die zurückgegebene Wert z1 ist total unrealistisch.
Das ist ein ziemlich kompliziertes Modell, mit dem ich nicht vertraut bin. Entspricht es nicht zwei Höckern? Ihre Eingabedaten haben einen großen Buckel und einen kleinen Buckel auf der linken Seite, aber es sollte schwer sein, den letzteren zu finden. Mein Verdacht ist, dass einer der Höcker in Ihrem Modell weit entfernt von Ihren Daten von Interesse ist, also im Wesentlichen passen Sie einen einzigen Buckel an Ihre Daten (daher die sinnlosen Parameter). Du könntest versuchen, diesen einen großen Buckel anzupassen, ihn von deinen Daten zu subtrahieren und dann einen anderen Buckel an den Rest anzupassen. –
Tatsächlich hinterlässt das Subtrahieren dieses Hauptpeaks einen Rauschdatensatz, der keine deutlich bemerkbaren Merkmale zu haben scheint (innerhalb von Fehlerbalken). Sind Sie sicher, dass das Modell gut zu Ihrem Problem passt? Es kann der Fall sein, dass ein Teil der Information irrelevant/unnötig ist (z. B. Anpassen der Summe von zwei Gauß-Werten an einen einzelnen Peak). –
Ja. Dies ist die Lichtverteilung entlang der Nebenachse einer Galaxie, die durch dieses Modell beschrieben werden soll. Ich bemerkte jedoch, dass die Ergebnisse glaubwürdiger sind, wenn ich die ersten 5 Datenpunkte loswerde. Vielleicht kann curve_fit nicht mit dieser ersten Beule umgehen ... – Pythoneer