A természeti jelenségek és fizikai törvények általában folytonos függvényekkel írhatók le. A méréseket vagy megfigyeléseket azonban csak diszkrét pontokban tudjuk végrehajtani. Ezért az adatok mindig csak véges számban állnak rendelkezésre, miközben az ezek közötti értékekre is gyakran szükség van, például további számításokhoz, vizualizációhoz vagy modellezéshez.
Az interpoláció célja, hogy ismert adatpontok alapján becslést adjunk közbenső, ismeretlen értékekre. Számos interpolációs módszer létezik: lineáris, polinom-, spline- és kernelalapú interpoláció.
A lineáris interpoláció a legegyszerűbb módszer, amelyben a szomszédos pontokat egyenes szakaszokkal kötjük össze.
A polinominterpoláció esetén egyetlen, magas fokszámú polinomot illesztünk az összes n mintapontra. A legismertebb változatok a Lagrange- és a Newton-féle interpolációk.
A spline interpoláció a polinomalapú módszer egy fejlettebb változata: az adatsort kisebb szakaszokra osztja, és minden szakaszra külön, alacsony fokszámú (gyakran harmadfokú) polinomot illeszt. A szomszédos szakaszok találkozási pontjainál biztosítja a folytonosságot és a deriváltak egyezőségét is.
A kernelalapú interpoláció olyan eljárás, amely súlyozott átlagot használ az interpolált érték kiszámítására. Minden mintapont egy kernel (magfüggvény) segítségével hat a környezetére. Az interpolált érték a mintapontok súlyozott összege, ahol a súlyokat a kernelfüggvény adja meg a távolság függvényében: a közelebbi pontok nagyobb, a távolabbiak kisebb súllyal szerepelnek.
Az interpolációs módszerek eltérnek a becslés pontosságában, a görbe simaságában és a számítási igényben. A lineáris interpoláció gyors, de töréspontokat eredményez a mintapontoknál. A többi módszer simább görbét ad, viszont több mintapont esetén a számításigény nő.
A továbbiakban a kernelalapú interpolációval foglalkozunk.
Ez a módszer számos előnnyel rendelkezik:
- A kernel alakja és paraméterei testreszabhatók az adott probléma jellegéhez.
- Az íveltebb görbéjű (például Gauss-jellegű) kernelek általában simább — matematikai nyelven: folytonosan differenciálható — interpolációs függvényt eredményeznek.
- Nem szükséges előzetesen megadni a modellstruktúrát, mint például a polinominterpolációnál.
- Kis mintaszám esetén is jól működhet, ha a kernel megfelelően van megválasztva.
Hátránya, hogy minden új becsült pont kiszámításakor a módszer újra meghatározza a kernel súlyokat az összes releváns mintaponthoz viszonyítva. Ez jelentős számítási igénnyel jár, ezért nagy adathalmazok esetén a módszer lassú lehet.
Bizonyos kernelfüggvények csak egy szűk tartományban vesznek fel nem nulla értéket. Ezeket a matematikai szakzsargonban kompakt tartójú kernelfüggvényeknek nevezik. Az interpoláció eredménye érzékenyen függ a kernelfüggvény típusától és a kompakt tartó hosszától is.
A kernelalapú interpoláció számításigényének csökkentésére most a következő megközelítést alkalmazzuk:
- egyenközű mintavétel.
- két minta közötti interpolációhoz csak a két minta értékét vesszük figyelembe.
Ehhez olyan K(x) kerneleket alkalmazunk, amelyek csak a [−1,1] intervallumon belül vesznek fel nem nulla értéket, azaz K(-1)=K(1)=0, továbbá a nullában 1 értékűek, azaz K(0)=1. Ekkor, ha a kerneleket az egyes adatpontokra (x= 0, 1, 2, 3…) helyezzük, akkor azok a szomszédos pozíciók mintáinak értékeire nem hatnak. Így, ha minden minta értékével megszorozzuk az adott pontra illesztett kernelt, és ezeket összegezzük, olyan interpolációs függvényt kapunk, amely az eredeti mintáknál azok pontos értékét adja vissza. - Annak érdekében, hogy egy konstans függvényt is helyesen interpoláljunk a kernelnek páros szimmetriájúnak kell lennie, valamint teljesülnie kell, hogy K(0.5) = K(-0.5) = 0.5. Ez biztosítja, hogy a két szomszédos pont súlyozott átlaga valóban az eredeti (konstans) értéket adja vissza középen is.
E feltételeknek megfelelő három eltérő kernelfüggvényt mutat az alábbi ábra. Ezek a háromszög alakú (triangular), a köbös (cubic) és az emelt koszinusz (raised cosine).
Az egyes grafikonokon egy adott kernelt két szomszédos mintapontra (x=0 és x=1) illesztve láthatunk, ahol az egyik minta 4, a másik 2 értékű. A két minta közötti értékeket a kernelek összegével becsüljük. A kiadódó interpolációs függvénygörbe kék színnel lett kirajzolva. Megfigyelhetjük, hogy a háromszög alakú kernel a lineáris interpolációt valósítja meg, ami egyébként képletekkel is egyszerűen bizonyítható.
A következőkben azt vizsgáljuk, hogy tetszőlegesen kiválasztott interpolálandó tesztfüggvényekre milyen pontos az interpoláció a három kernel esetén. A pontosságot az egzakt függvényértékek és az interpolált értékek négyzetes eltérésének átlagával (mean squared error, MSE) fogjuk mérni, amihez egy mse() nevű függvényt írunk.
Az egyenközű mintavétel után előálló x és y értékeket egy sorozattípusú konténerben tároljuk, amelyeket jelölje xs és ys. Ahhoz, hogy megtudjuk, hogy egy adott vizsgált x érték mely két szomszédos minta által meghatározott intervallumba, vagyis az xs sorozat mely két indexe közé esik, szükségünk van egy gyors indexkereső függvényre. Ennek megvalósításában az O(log n) időkomplexitású bináris keresést megvalósító bisect() függvényt használjuk fel, amely a szabványos könyvtár bisect moduljában található.
Az interpoláció eredményét grafikonon is meg akarjuk jeleníteni. Ehhez is készítünk egy függvényt.
Az eddig említett függvények még nem az interpolációt végzik, hanem csak közreműködők a célunk elérésében. Ezért ezeket a segédfüggvényeket külön, egy util nevű modulban gyűjtjük össze.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 |
# modul: util # Python 3.12+ from typing import Sequence, Callable import matplotlib.pyplot as plt from bisect import bisect from math import isclose, gamma, exp type RealNum = int | float type UnaryRealFunction = Callable[[RealNum], RealNum] def find_interval_indexes(seq: Sequence[RealNum], x: RealNum) -> tuple[int, int] | None: """Meghatározza, hogy a seq rendezett, ismétlődő elemet nem tartalmazó sorozat mely két szomszédos indexű eleme közé esik a megadott x érték. Ha ilyen elemet nem talál, a visszatérési érték None. Ha az x az egyik elem értékével egyenlő, akkor a visszaadott két index azonos. """ if not seq or x < seq[0] or x > seq[-1]: return None index_high = bisect(seq, x) # Gyors indexmeghatározást tesz lehetővé (időkomplexitás O(log n)). index_low = index_high - 1 if isclose(seq[index_low], x): index_high = index_low return index_low, index_high def mse(seq1: Sequence[RealNum], seq2: Sequence[RealNum]) -> float: """Az argumentumként megadott sorozatok értékeinek átlagos négyzetes hibájával (mean squared error) tér vissza.""" return sum((y1 - y2) ** 2 for y1, y2 in zip(seq1, seq2)) / len(seq1) def gamma_dist(x, a, b): if x <= 0: raise ValueError('Az x pozitív kell, hogy legyen.') return b ** a / gamma(a) * x ** (a - 1) * exp(-b * x) def plot_interpolation(_ys_exact, _ys_interpolated, _x_values, kernel_func): plt.plot(_x_values, _ys_exact, color='red', marker='o') plt.plot(_x_values, _ys_interpolated, linestyle='', marker='.', markerfacecolor='black') plt.title(f'Kernel: {kernel_func.__name__}') plt.grid(True) plt.xlabel('X') plt.ylabel('Y') plt.show() def plot_adjacent_kernels_with_sum(kernel_func: UnaryRealFunction, amp1: RealNum, amp2: RealNum, a: RealNum = -2, b: RealNum = 2, n: int = 200): delta = (b - a) / (n - 1) _xs = [a + i * delta for i in range(n)] _ys1 = [amp1 * kernel_func(x) for x in _xs] _ys2 = [amp2 * kernel_func(x - 1) for x in _xs] _xs_chunk = [x for x in _xs if 0 <= x <= 1] _ys1_chunk = [amp1 * kernel_func(x) for x in _xs_chunk] _ys2_chunk = [amp2 * kernel_func(x - 1) for x in _xs_chunk] _ys12sum = [y1 + y2 for y1, y2 in zip(_ys1_chunk, _ys2_chunk)] plt.plot(_xs, _ys1, linestyle='solid', color='red', label='Kernel az első adatponton') plt.plot(_xs, _ys2, linestyle='solid', color='green', label='Kernel a második adatponton') plt.plot(_xs_chunk, _ys12sum, linestyle='solid', color='blue', linewidth=2, label='Kernelek összege a [0,1] intervallumban') plt.title(f'Kernel: {kernel_func.__name__}') plt.xlabel('X') plt.ylabel('Y') plt.minorticks_on() plt.grid(which='major', linestyle='-', linewidth='0.5', color='black') # Fő rácsvonalak plt.grid(which='minor', linestyle=':', linewidth='0.5', color='gray') # Segéd rácsvonalak plt.legend(loc='center left', fontsize='x-small') plt.show() if __name__ == '__main__': from kernels import triangular, cubic, raised_cosine for kernel in (triangular, raised_cosine, cubic): plot_adjacent_kernels_with_sum(kernel, 4, 2) |
Szintén külön modulban definiáljuk a kernelfüggvényeket. E modul neve kernels lesz.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 |
# modul: kernels # Az iterpoláció kernelfüggvényei. from math import pi, cos def triangular(x): """Háromszög kernel, amelynek a tartója a [–1,+1] zárt intervallum.""" return 1 - abs(x) if abs(x) <= 1 else 0 def raised_cosine(x): """Emelt koszinusz kernel, amelynek a tartója a [–1,+1] zárt intervallum.""" return (cos(x * pi) + 1) / 2 if abs(x) <= 1 else 0 def cubic(x): """Harmadfokú polinom kernel, amelynek a tartója a [–1,+1] zárt intervallum.""" abs_x = abs(x) return 1 - abs_x + abs_x * (1 - abs_x) * (0.5 - abs_x) if abs(x) <= 1 else 0 |
A main modulban definiáljuk az egyenközű mintavételt és az interpolációt végző függvényeket, valamint a tesztelő függvényt. A main természetesen használja a kernels és a util modulokat, ezért ezeket beimportáljuk. Az eddig leírtak és a részletes kommentek segítik a megértést.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 |
# modul: main # Python 3.12+ from math import pi, sin, exp from util import RealNum, UnaryRealFunction, Sequence, find_interval_indexes, mse, gamma_dist, plot_interpolation from kernels import triangular, cubic, raised_cosine def evenly_sample(fn: UnaryRealFunction, a: RealNum, b: RealNum, n: int) -> tuple[list[RealNum], list[RealNum]]: """Az fn egyváltozós valós függvényből n darab (n>=2) egyenközű mintát vesz a megadott [a, b] zárt intervallumban. Visszatérési értéke két lista: a független változó (x) értékeinek listája és a hozzájuk tartozó függvényértékek (fn(x)) listája. """ if n <= 1: raise ValueError("A minták száma (n) legalább kettő kell, hogy legyen.") delta = (b - a) / (n - 1) _xs = [a + i * delta for i in range(n)] _ys = [fn(x) for x in _xs] return _xs, _ys def compact_kernel_interpolation(x, _xs: Sequence[RealNum], _ys: Sequence[RealNum], kernel_func: UnaryRealFunction): """Az egyenközű mintavételezéssel kapott x és y értékeket tartalmazó, egyenlő hosszú, rendezett _xs és _ys sorozatok alapján interpolálja az x-hez tartozó értéket a megadott kernel függvényt használva. A kernel a [-1, 1] intervallumon kívül zérus értékű. """ if len(_xs) != len(_ys): raise ValueError('A két sorozatnak egyenlő hosszúnak kell lenni.') # Meghatározzuk, hogy az x melyik intervallumba esik. index_low, index_high = find_interval_indexes(_xs, x) # Ha x pontosan egybeesik egy mintaponttal nincs interpoláció, hanem az egzakt mintaértéket adjuk vissza. if index_low == index_high: return _ys[index_low] # Ha x két minta közé esik, akkor interpolálunk a két minta alapján az adott kernellel. z = (x - _xs[0]) / (_xs[1] - _xs[0]) return _ys[index_low] * kernel_func(z - index_low) + _ys[index_high] * kernel_func(z - index_high) def test_compact_kernel_interpolation(a, b, fn_to_interpolate): xs, ys = evenly_sample(fn_to_interpolate, a, b, 15) # Az interpolálandó függvény mintái. xs_test, _ = evenly_sample(fn_to_interpolate, a, b, 24) # A vizsgált x értékek. ys_exact = [fn_to_interpolate(x) for x in xs_test] # Az interpolálandó függvény értékei a vizsgált x értékeknél. # Eltérő kernelekkel interpolálunk. Megjelenítjük a függvényt és az iterpolációs pontokat, és # kiszámítjuk a négyzetes hibák átlagát. for kernel in (triangular, cubic, raised_cosine): # Interpolált értékek meghatározása. ys_interp = [compact_kernel_interpolation(x, xs, ys, kernel) for x in xs_test] # Megjelenítés grafikonon. plot_interpolation(ys_exact, ys_interp, xs_test, kernel) # Négyzetes hibaátlag számítás. print(f'{kernel.__name__:>15}: {mse(ys_exact, ys_interp) * 10000:.2f}') # TESZT for a, b, target_func in [(-pi, +pi, sin), (-5, 5, lambda x: exp(-x ** 2 / 4)), (0.1, 15, lambda x: 10 * gamma_dist(x, 7, 1))]: test_compact_kernel_interpolation(a, b, target_func) |
Három tesztfüggvényre végezett interpolációk eredményét láthatjuk a következő ábrán grafikonon ábrázolva, ahol feltüntettük az egyes kernelekkel elért négyzetes hibaátlagokat is. (A könnyebb olvashatóság és összehasonlítás érdekében a tényleges értéket egy adott konstanssal szorozva)
Az eredmények azt mutatják, hogy minden vizsgált esetben a háromszög alakú kernel – azaz a lineáris interpoláció – adja a legkisebb átlagos négyzetes hibát. Más tesztfüggvényekkel végzett próbák is ugyanezt az eredményt hozzák.
Ez elsőre meglepő lehet, mivel a kernelalapú interpolációban általában a lekerekített görbéjű kernelek simább és pontosabb interpolációs függvényt eredményeznek. Ez a simaság azonban nagyobb számításigénnyel jár.
Mivel jelen esetben a kernelek csak két mintapontot használnak az interpolációhoz, a simább kernelek előnye nem tud érvényesülni; így a legegyszerűbb, lineáris megközelítés bizonyul a leghatékonyabbnak.
E bejegyzésben a matematika egy erős elméleti alapokon nyugvó, és a gyakorlatban számos területen alkalmazott módszerének kis szeletét érintetettük, amelynek numerikus vizsgálatához a Python megfelelő nyelvi elemeit használtuk. Ezekről a Python tudásépítés lépésről lépésre című e-könyv „Készétel fogyasztás – a szabványos könyvtár moduljainak használata” fejezet „Matematikai számítások támogatása” alfejezetben lehet bővebben olvasni.