Základy programování v Matlabu - část 7
Ing. Petr Pokorný, Ph.D.
Základy aproximace dat
V řadě praktických aplikací je třeba diskrétní data (např. měřené body) vyjádřit vhodnou funkční závislostí, která bude daný jev (model) co nejlépe vyjadřovat. K tomu slouží tzv. aproximace dat.
Aproximovat data můžeme řadou způsobů. Asi nejčastěji používaným nástrojem je ale tzv. metoda nejmenších čtverců (MNČ), která minimalizuje sumu čtverců odchylek zadaných diskrétních dat od výsledného modelu. Pro podrobnější popis a odvození odkážeme na dokument metoda_nejmensich_ctvercu.html. Jiným často používaným typem je tzv. Čebyševova aproximace, která minimalizuje maximum absolutní hodnoty rozdílu mezi daty a hledanou funkcí.
Pro aproximaci je nutné znát (předpokládat), jaký model měřená data vyjadřují, resp. jakou funkční závislost chceme danými diskrétními daty prokládat. Obecně můžeme úlohu rozdělit na:
- Aproximace lineární kombinace bázových funkcí,
- Nelineární aproximace.
V prvním případě předpokládáme, že daná závislost je vyjádřena lineární kombinací vybraných funkcí a cílem aproximace je určit koeficienty této lineární závislosti.
Nejjednodušší případ je zřejmě polynomická aproximace N-tého stupně pro funkci jedné proměnné, kdy předpokládáme, že diskrétní data budou reprezentována funkcí
,
,kde
jsou hledané koeficienty. Pro výpočet zmíněné interpolace v Matlabu je implementována funkce polyfit.
y_original = 0.01*x_original.^3 - 0.1*x_original.^2 + 2*x_original + 0.001 + randn( size(x_original) );
% aproximace polynomem 3. stupne
p = polyfit( x_original, y_original, 3 );
% vycisleni aproximovanych hodnot
y_fit = polyval( p, x_original );
% rozdil aproximovanych a puvodnich hodnot
rms = sqrt( sum( e.^2 )/length( e ) );
plot( x_original, y_original, 'x', x_original, y_fit, '-' )
title( [ 'RMS = ', num2str( rms, '%.2f' ) ] )
Jako bázové funkce ale můžeme obecně volit libovolné druhy polynomů (Legendreovy, Čebyševovy, apod.) nebo funkcí. Pro nalezení koeficientů lineární kombinace poté ručně tvoříme algoritmus MNČ.
% Legendreovy polynomy Pn(x) do stupne 3
% (https://cs.wikipedia.org/wiki/Legendreovy_polynomy)
P2 = @(x) 0.5*( 3*x.^2 - 1 );
P3 = @(x) 0.5*( 5*x.^3 - 3*x );
% sestaveni matic pro MNC
A = [ P0(x_original'), P1(x_original'), P2(x_original'), P3(x_original') ];
L_fit_Legendre = A*p_Legendre;
e_legendre = L - L_fit_Legendre;
rms_legendre = sqrt( sum( e_legendre.^2 )/length( e_legendre ) );
plot( x_original, L, 'x', x_original, L_fit_Legendre, '-' )
title( [ 'RMS = ', num2str( rms_legendre, '%.2f' ) ] )
V případě aproximace nelineární závislosti řešíme úlohu buď iterativně s využitím Taylorova rozvoje, kdy je třeba základní rovnice MNČ modifikovat, nebo s použitím optimalizačních algoritmů (viz dále).
Interpolace v systému Matlab
Pro interpolaci dat lze využít výše uvedené aproximace, a poté vyčíslit hodnotu funkce v požadovaném bodě (interpolační aproximace). Existuje opět velká řada interpolačních algoritmů, pro jejich podrobný výčet odkážeme např. na http://kfe.fjfi.cvut.cz/~limpouch/numet/aprox.pdf . V prostředí Matlab je k dispozici řada předpřipravených nástrojů pro interpolaci (viz doc MATLAB/Mathematics/Interpolation): - interp1 - 1-D interpolace (vyhledání v tabulce)
- interp2 - interpolace 2-D dat v meshgrid formátu
- interp3 - interpolace 3-D dat v meshgrid formátu
- interpn - interpolace N-D dat v meshgrid formátu
- pchip - po částech kubická Hermiteovská interpolace
- makima - modifikovaná Akima po částech kubická Hermiteovská interpolace
- spline - kubická spline interpolace
- ...
y_i1 = interp1( x_original, y_original, x_i )
y_is = spline( x_original, y_original, x_i )
x_fine = linspace( min(x_original), max(x_original), 100 );
y_interp1 = interp1( x_original, y_original, x_fine );
y_spline = spline( x_original, y_original, x_fine );
plot( x_original, y_original, 'x', x_fine, y_spline, '-', x_i, y_is, 'o', x_fine, y_interp1, '--', x_i, y_i1, '^' )
Příklady
Příklad A1
Data, která by měla odpovídat kvadratické závislosti, jsou dána vektory x a y. Aproximujte je vhodným polynomem a určete, od jakého
je sklon závislosti menší než
. x_original = [ 0, 1.1111e-01, 2.2222e-01, 3.3333e-01, 4.4444e-01, 5.5556e-01, 6.6667e-01, 7.7778e-01, 8.8889e-01, 1.0000e+00 ];
y_original = [ -2.8934e-02, -1.1082e-03, 5.7065e-02, 4.3140e-02, 8.9658e-03, -2.3682e-02, -5.6995e-02, -5.3039e-02, -7.1372e-02, -1.1040e-01 ];
Příklad A2
Zadaná data x a y zobrazte na intervalu
s dělením po 1/100 pomocí kubické interpolace a porovnejte s interpolací makima a spline. Příklad A3
Níže uvedený obrázek peppers.png je načtený (imread), má snížené rozlišení (imresize) a je převeden do odstínů šedi (rgb2gray). Pomocí kubické interpolace a funkce interp2 vypočtěte obrázek nový, který bude mít rozlišení trojnásobné. Využijte také jiné metody interpolace (např. lineární, spline nebo makima).
Nápověda: Pro použití funkce interp2 je třeba mít data určená na vhodném gridu souřadnic. Ten můžete vytvořit např. pomocí funkce ndgrid nebo meshgrid.
I = double( rgb2gray( imresize( imread( 'peppers.png' ), [ 50 100 ] ) ) );
Úvod do optimalizačních (minimalizačních) algoritmů
Optimalizační úlohy jsou numerické metody, které umožňují nalézt lokální (někdy i globální) minimum (maximum) dané funkce. Tyto algoritmy nacházejí výrazné uplatnění v řadě moderních technických aplikací a je jich vyvinuto velké množství. Matlab má několik těchto algoritmů implementovaných (Optimization Toolbox). Uvažujeme-li funkci jedné proměnné danou předpisem
, potom cílem algoritmu je nalézt takové x, pro které bude funkce y minimální. Pro většinu algoritmů je nutné zadat počáteční hodnoty, ze kterých daný algoritmus vychází. Příklady
Příklad O1: Metoda největšího spádu (příklad pro funkci jedné proměnné)
V tomto příkladu budeme demonstrovat zjednodušený případ optimalizace metodou největšího spádu.
Metoda největšího spádu je jedna ze základních optimalizačních úloh, která vychází z následujícího jednoduchého iteračního alogritmu pro funkci
:
,kde
je předchozí hodnota nezávisle proměnné (na počátku volený odhad polohy minima), e je kladné reálné číslo,
je gradient funkce (diference) v bodě
. Napište funkci pro optimalizaci (hledání minima funkce) výše popsanou metodou. K testování využijte funkci
. Schéma možného algoritmu největšího spádu:
- Nastavení počáteční hodnoty
, tolerance
, po kterou má probíhat výpočet, a koeficientu e (optimální hodnotu koeficientu lze také vypočítat, pro potřeby tohoto demonstračního příkladu budeme ale hodnoty volit manuálně), - výpočet gradientu v bodě
(nejprve v počátečním bodě, v dalších iteracích v následujících vypočtených bodech), - výpočet
a odpovídající funkční hodnoty
, - testování tolerance:
.
Příklad je možné řešit s pomocí níže uvedených dvou uživatelských funkcí:
function [xmin,ymin,it] = stdesc(f,x0,tol,e)
% Funkce řešící optimalizaci: [xmin,ymin,it] = stdesc(f,x0,tol,e)
% Vstupní parametry: f reference na uživatelem definovanou funkci k minimalizaci, tato funkce musí vracet jak hodnotu v bodě x, tak hodnotu gradientu v tomto bodě,
% x0 počáteční odhad polohy minima,
% Výstupní parametry: xmin poloha minima,
% it počet iterací k výpočtu.
% Funkce k minimalizaci: [M, G] = fun(x)
% Vstupní parametry: x nezávisle proměnná
% Výstupní parametry: M hodnota funkce v bodě x,
% G hodnota gradientu v bodě x.
V těle programu budeme volat funkci optimalizace stdesc, ve které bude zakomponován iterativní výše zmíněný proces, pomocí příkazu [xmin,ymin,it] = stdesc(@fun,x0,tol,e), kde x0, tol a e jsou předem definované parametry.
Příklad O2: Optimalizace - kořeny funkcí
Pomocí funkce fzero nalezněte kořeny funkce
.Vyzkoušejte jak zadání anonymní funkce ( fzero(@(x) (cos(exp(x)) …), 1) ), tak definici uživatelsky definované funkce. Nalezené kořeny ověřte ve grafu.
Příklad O3: Optimalizace - minimalizace na omezeném intervalu
Na intervalu
nalezněte minimum funkce jedné proměnné (fminbnd)
.Pomocí grafu zkontrolujte, zda se opravdu jedná o minimum.
Příklad O4: Optimalizace - minimalizace na neomezeném intervalu
Pomocí funkce fminsearch hledejte minimum funkce
.pro počáteční body
,
, 0 a 1. Výsledky zobrazte do přehledného grafu. Co z těchto výsledků plyne? Komentář: Jak je patrné, nalezení globálního minima je zejména pro rychle oscilující funkce komplikované. Zkuste modifikovat řešení takovým způsobem, že na počátku vygenerujete N náhodných počátečních bodů, ze kterých bude optimalizace vycházet. Z výsledků N optimalizačních procedur jako minimum vyberte to, pro které je funkční hodnota nejmenší. S narůstajícím počtem generovaných počátečních bodů bude narůstat pravděpodobnost, že výsledek optimalizace bude hledaným minimem na zadaném intervalu. I přesto ale nemusí být řešení ideální. Pozn.: Vhodné je opatřit program takovou podmínkou, která kontroluje, že se v dané sekvenci vygenerovaných bodů neopakují stejné souřadnice.
Příklad O5: Optimalizace - aproximace funkcí
Pomocí vhodného optimalizačního algoritmu nalezněte koeficienty mocninné řady
takové, aby řada co nejlépe aproximovala funkci
na intervalu
pro
ve smyslu minimalizace absolutních hodnot rozdílů řady od původní funkce. Výsledek optimalizace porovnejte se skutečnými hodnotami v grafu. Návod k řešení: Koeficienty řady budou nyní tvořit proměnné pro cílovou funkci M, které se budou měnit. Cílovou funkci M můžeme napsat ve tvaru
,kde první sumace probíhá přes všechny hodnoty x. Cílem je nalézt takové koeficienty
, pro které bude cílová funkce minimální. Schematické znázornění algoritmu může být následující: - definice globálních proměnných (global) x a N,
- tvorba dostatečně jemně děleného vektoru x a zadání maximálního stupně polynomu N,
- volba počátečních hodnot koeficientů řady (na této volbě závisí úspěšnost algoritmu, pro jednoduchost můžete volit počáteční koeficienty nulové),
- tvorba uživatelsky definované funkce M = merit(a), na jejímž vstupu je vektor koeficientů řady, dále definovány globální proměnné x a N (předávané ze základní funkce); funkce merit() bude vracet hodnotu cílové funkce podle výše zmíněného vzorce,
- minimalizace pomocí vhodného optimalizačního algoritmu (např. fminsearch), nalezeny proměnné koeficienty řady takové, aby cílová funkce byla minimální.
Curve Fitting Toolbox
V systému Matlab je aproximace dat řešena také v Curve Fitting Toolbox (viz dokumentace). Grafický nástroj pro aproximaci lze otevřít pomocí ze záložky Apps/Curve Fitting (doc cftool).