Függvény illesztés

</h1>

Ebben a notebookban bemutatjuk, hogyan kell az illesztési feladatokat elvégezni.

In [1]:
%pylab inline
figsize(6,6) #Képméret megváltoztatása

Populating the interactive namespace from numpy and matplotlib

Függvény illesztés adatokra

Az adatainkra is ugyanolyan függvényt kell illeszteni, mint amit fentebb megismertünk. A függvények illesztéséhez a curve_fit eljárást használjuk, a scipy csomagból, mely a Gnuplot-hoz hasonló elven keresi meg a paraméterek ideális értékét.

In [2]:
adat=loadtxt("sinusadatok.dat")
x=adat[:,0]
y=adat[:,1]
z=adat[:,2]
In [3]:
from scipy.optimize import curve_fit
figsize(10,10)

def func(x, a,b,c,d):
    return a*cos((d*(x+c)))+b

parameter, covariance_matrix = curve_fit(func, x, z)

Az illesztett paraméterek:

In [4]:
print(parameter)

[-1.21224409  2.60067796  1.22001586  1.00095743]

A kovariancia mátrix diagonális elemei megegyeznek a hiba négyzetével, tehát az illesztési hiba így számolandó:

In [5]:
print(sqrt(diag(covariance_matrix)))

[ 0.00795894  0.00558654  0.01376697  0.0013635 ]

Kicsit "szofisztikáltabban" így lehet leríni:

In [6]:
print("Amplitúdó:",parameter[0], "+/-", sqrt(diag(covariance_matrix)[0]))
print("y eltolás:",parameter[1], "+/-", sqrt(diag(covariance_matrix)[1]))
print("x eltolás:",parameter[2], "+/-", sqrt(diag(covariance_matrix)[2]))
print("Fázisnyújtás:",parameter[3], "+/-", sqrt(diag(covariance_matrix)[3]))

Amplitúdó: -1.21224408966 +/- 0.00795893637343
y eltolás: 2.6006779641 +/- 0.00558654343834
x eltolás: 1.22001586347 +/- 0.0137669722977
Fázisnyújtás: 1.000957427 +/- 0.00136350325803
In [7]:
plot(x,z, "go")
xlabel("Idő [sec]")
ylabel("Elmozdulás [m]")
title("Lengés")
grid(True)
plot(x, func(x, *parameter), 'r-', label='fit',linewidth=5.0)
Out[7]:
[<matplotlib.lines.Line2D at 0x7f3ec90d0898>]

Kezdeti érték és hiba megadása

Lehetőség van kezdőfeltételt is megadni az illesztéshez. Erre a "p0" parameter szolgál. A mérési hibát a "sigma" paraméter után kell megadni.

In [8]:
error=randn(len(z)) # normál eloszlású random tömböt választunk hibának
print(error[0])

1.27196942432
In [9]:
# Megadunk kezdőfeltétlet, és azt mondjuk, hogy az előbb számolt hiba a mérési hibánk.
param, cov = curve_fit(func, x, z, p0=(1.,3.,0,1), sigma=error)
In [10]:
figure(figsize=(6,6))
plot(x,z, "go")
xlabel("Idő [sec]")
ylabel("Elmozdulás [m]")
title("Lengés")
grid(True)
plot(x, func(x, *param), 'r-', label='fit',linewidth=5.0)
Out[10]:
[<matplotlib.lines.Line2D at 0x7f3ef6381240>]

Az eredmények: (biztos egyeznek?)

In [11]:
print("Amplitúdó:",param[0], "+/-", sqrt(diag(cov)[0]))
print("y eltolás:",param[1], "+/-", sqrt(diag(cov)[1]))
print("x eltolás:",param[2], "+/-", sqrt(diag(cov)[2]))
print("Fázisnyújtás:",param[3], "+/-", sqrt(diag(cov)[3]))

Amplitúdó: 1.24397081059 +/- 0.0118207573735
y eltolás: 2.66489448411 +/- 0.00773304298654
x eltolás: -1.80530982223 +/- 0.0123644586756
Fázisnyújtás: 0.974388916202 +/- 0.00147777187424

Paraméter becslés

Gyakran előfordul, hogy olyan függvényt kell illeszteni, ami nehezen találja meg a megfelelő illeszkedést (erősen függ a kezdeti értékek helyes megadásától). Az előre definiált kezdő értékekből gyakran rossz eredményt fogunk kapni. Ekkor szükséges olyan, közel helyes paramétereket beállítani a függvénynél, ahonnan könnyen megtalálja az algoritmus a helyes megoldást.

Található a mappában egy "fitting.csv" állomány, ami egy csillapított rezgőmozgás adatsorát tartalmazza. Az illesztendő függvény alakja: $$f(x)=A*sin(2\pi\cdot\frac{x}{T} -\phi)*exp(\frac{-x}{\tau})+F0 , $$ ahol A az amplitúdót, $\phi$ a kezdő fázist, T a periódusidőt, $\tau$ a csillapítást és F0 a nyugalmi kitérést mutatja.

A függvény bonyolultsága miatt az illesztést csakis a kezdőparaméterek helyes becslése után lehet jól elvégezni. A tapasztalatok azt mutatják, hogy a periódusidő becslése a legfontosabb az ilyen periodikus jelek illesztése során, így a T jó megválasztása nélkülözhetetlen a sikeres illesztéshez.

Első lépésben nézzünk bele a fájlba:

In [12]:
%%bash
head fitting.csv

0.0000e+00 1.8744e+02
2.5025e-02 2.0521e+02
5.0050e-02 2.1968e+02
7.5075e-02 1.9634e+02
1.0010e-01 2.0968e+02
1.2513e-01 1.9958e+02
1.5015e-01 2.2859e+02
1.7518e-01 1.9386e+02
2.0020e-01 1.8352e+02
2.2523e-01 1.6855e+02

Rajzoljuk ki az adatokat, hogy meg tudjuk becsülni a paraméterek értékét:

In [13]:
fitting=loadtxt('fitting.csv')
x=fitting[:,0]
data=fitting[:,1]
plot(fitting[:,0],fitting[:,1], ".-")
grid()

Próbáljuk meg a szokott módon az illesztést:

In [14]:
def func(x, amp,t,phi,tau,f0):
    return amp*sin(x/t*2*pi-phi)*exp(-x/tau)+f0

parameter, covariance_matrix = curve_fit(func, x, data)
print(parameter)

figure(figsize=(10,10))
plot(x,data, "g.")
xlabel("Idő [sec]")
ylabel("Elmozdulás [m]")
title("Lengés")
grid(True)
plot(x, func(x, *parameter), 'r-', label='fit',linewidth=1.0)

[  1.65149700e+00   2.48210323e+01  -1.03505300e+02  -1.38811967e+05
   1.19819880e+02]
Out[14]:
[<matplotlib.lines.Line2D at 0x7f3ecbd19518>]

Mivel a fenti illesztés láthatóan nem helyes, adjunk meg közel helyes kezdőértéket az illesztéshez. Ehhez a curve_fit parancs p0 kapcsolóját kell ismét használni.

In [15]:
def func(x, amp,t,phi,tau,f0):
    return amp*sin(x/t*2*pi-phi)*exp(-x/tau)+f0

parameter, covariance_matrix = curve_fit(func, x, data, p0=[95., 2.5, -pi/2, 50.,120.])
print(parameter)

figure(figsize=(10,10))
plot(x,data, "g*")
xlabel("Idő [sec]")
ylabel("Elmozdulás [m]")
title("Lengés")
grid(True)
plot(x, func(x, *parameter), 'r-', label='fit',linewidth=2.0)

[  89.42072117    2.49985936   -1.58176518   61.27361296  119.74623238]
Out[15]:
[<matplotlib.lines.Line2D at 0x7f3eca6ca7b8>]