03-Csúcs keresés¶
A ../data/pupos file tartalmaz egy zajos adatsort.
Illesszünk erre az adatsorra két Gauss-görbe összegét! $$ f_{2G}(x)=A_1 \mathrm{e}^{-(x-e_1)^2/s_1^2}+A_2 \mathrm{e}^{-(x-e_2)^2/s_2^2} $$
Illesszünk erre az adatsorra két Lorentz-görbe összegét is! $$ f_{2L}(x)=\frac{B_1}{w_1^2+(x-p_1)^2}+\frac{B_2}{w_2^2+(x-p_2)^2} $$
Ábrázoljuk az eredeti adatsort, illetve mind a két esetben ábrázoljuk az illesztés eredményét is! Ha az illesztés során esetleg problémákba ütközünk, akkor próbálkozzunk az illesztési paraméterek kezdeti értékének becslésével!
Megoldás¶
Először betöltjük a szükséges csomagokat
%pylab inline
from scipy.optimize import curve_fit
Majd az adatokat is a fájlból:
adatok = loadtxt("../data/pupos.txt")
plot(adatok[:,0], adatok[:,1])
xlabel("x tengely", size=12)
ylabel("y tengely", size=12)
title("A fájl adatai", size=18, y=1.05)
Ezek után definiáljuk a szükséges függvényeket, amiket később használni fogunk:
def f1g(x, A1, e1, s1):
'''Visszaadja a megadott paraméterek szerinti Gauss-görbét'''
return A1 * exp(-(x - e1)**2/s1**2)
def f1l(x, B1, p1, w1):
'''Visszaadja a megadott paraméterek szerinti Lorentz-görbét'''
return B1/(w1**2 + (x-p1)**2)
def f2g(x, A1, e1, s1, A2, e2, s2):
'''Visszaadja a megadott paraméterek szerinti Gauss-görbék összegét'''
return A1 * exp(-(x - e1)**2/s1**2) + A2 * exp(-(x - e2)**2/s2**2)
def f2l(x, B1, p1, w1, B2, p2, w2):
'''Visszaadja a megadott paraméterek szerinti Lorentz-görbék összegét'''
return B1/(w1**2 + (x-p1)**2) + B2/(w2**2 + (x-p2)**2)
Gauss¶
Először illesszük meg a Gauss-görbe-összeget! Ehhez meg kell becsüljük a paraméterek értékét. Első ránézésre van egy 4.5 körüli magasság vékonyabb csúcs 3 körül, és egy 2 magasságú szélesebb csúcs 5.5 körül:
plot(adatok[:,0], adatok[:,1], label="Az adatfájl adatai")
plot(adatok[:,0], f1g(adatok[:,0], 4.5, 3, 1), label="Az egyik becsült Gauss-görbe")
plot(adatok[:,0], f1g(adatok[:,0], 2, 5.5, 3), label="A másik becsült Gauss-görbe")
xlabel("x tengely", size=12)
ylabel("y tengely", size=12)
yticks(linspace(-1, 7, 9))
title("Az adatokra becsült Gauss-görbék", size=18, y=1.05)
legend(loc=0)
Ezekkel a paraméterekkel elvégezzük az illesztést:
gfit = curve_fit(f2g, adatok[:,0], adatok[:,1], p0=[4.5, 3, 1, 2, 5.5, 3])
print("illesztési paraméterek: {}".format(gfit[0]))
plot(adatok[:,0], adatok[:,1], label="Az adatfájl adatai")
plot(adatok[:,0], f2g(adatok[:,0], gfit[0][0], gfit[0][1], gfit[0][2], gfit[0][3], gfit[0][4], gfit[0][5]), "c-", lw=2.5, label="Illesztett görbe")
plot(adatok[:,0], f1g(adatok[:,0], gfit[0][0], gfit[0][1], gfit[0][2]), "--", label="Egyik tag")
plot(adatok[:,0], f1g(adatok[:,0], gfit[0][3], gfit[0][4], gfit[0][5]), "--", label="Másik tag")
xlabel("x tengely", size=12)
ylabel("y tengely", size=12)
title("Az adatokra illesztett Gauss-görbe", size=18, y=1.05)
legend(loc=0)
Lorentz¶
Ezután illesszük a Lorentz-görbe összeget! Itt is szükség van értékbecslésre, a csúcsok pedig ugyanott vannak, mint az előbb, de a paraméterek értékei kicsit módosulnak. A becsléshez most is támaszkodhatunk az ábrára:
plot(adatok[:,0], adatok[:,1], label="Az adatfájl adatai")
plot(adatok[:,0], f1l(adatok[:,0], 4.5, 3, 1), label="Az egyik becsült Lorentz-görbe")
plot(adatok[:,0], f1l(adatok[:,0], 8, 5.5, 2), label="A másik becsült Lorentz-görbe")
xlabel("x tengely", size=12)
ylabel("y tengely", size=12)
yticks(linspace(-1, 7, 9))
title("Az adatokra becsült Lorentz-görbék", size=18, y=1.05)
legend(loc=0)
A paraméterek megbecslése után már illeszthetjük a görbét:
lfit = curve_fit(f2l, adatok[:,0], adatok[:,1], p0=[4.5, 3, 1, 8, 5.5, 2])
print("illesztési paraméterek: {}".format(gfit[0]))
plot(adatok[:,0], adatok[:,1], label="Az adatfájl adatai")
plot(adatok[:,0], f2l(adatok[:,0], lfit[0][0], lfit[0][1], lfit[0][2], lfit[0][3], lfit[0][4], lfit[0][5]), "c-", lw=2.5, label="Illesztett görbe")
plot(adatok[:,0], f1l(adatok[:,0], lfit[0][0], lfit[0][1], lfit[0][2]), "--", label="Egyik tag")
plot(adatok[:,0], f1l(adatok[:,0], lfit[0][3], lfit[0][4], lfit[0][5]), "--", label="Másik tag")
xlabel("x tengely", size=12)
ylabel("y tengely", size=12)
title("Az adatokra illesztett Lorentz-görbe", size=18, y=1.05)
legend(loc=0)
