PA081: Programování numerických výpočtů 8. Modelování experimentálních dat, metoda nejmenších čtverců Aleš Křenek jaro 2019 PA081: Programování numerických výpočtů A. Křenek Modelování experimentálních dat Metoda nejmenších čtverců Lineární regrese Obecný model Lineární modely SAXS Vícerozměrná data Metoda nejmenších absolutních odchylek Regularizace 1/34 Modelování experimentálních dat ► v experimentu naměříme v bodech xi hodnoty y i ► x může být libovolná veličina: čas, napětí, poloha, ... ► chování systému popisujeme modelem y = M(x) ► model závisí na sadě parametrů au tj. y = M(x, ai,..., um) ► hledáme takové hodnoty au pro něž model odpovídá nejlépe experimentu PA081: Programování numerických výpočtů A. Křenek Modelování experimentálních dat Metoda nejmenších čtverců Lineární regrese Obecný model Lineární modely SAXS Vícerozměrná data Metoda nejmenších absolutních odchylek Regularizace Modelování experimentálních dat Příklad řln2 ► radioaktivní rozpad N = Noe t ► N je počet atomů ve vzorku (No v čase t = 0), T poločas rozpadu 25 20 15 - 10 5 - \ \ "priklady/exp.dat" + 20*exp(-x/20*log(2)) ----- + + J_L 10 20 30 40 50 60 70 80 90 100 PA081: Programování numerických výpočtů A. Křenek Modelování experimentálních dat Metoda nejmenších čtverců Lineární regrese Obecný model Lineární modely SAXS Vícerozměrná data Metoda nejmenších absolutních odchylek Regularizace 3/34 Metoda nej menších čtverců Odvození ► „Jaká je pravděpodobnost, že konkrétni sada parametru at je správná?" ► špatně položená otázka ► neexistuje „náhodná veličina modelů" ► naopak, náhodnou veličinou jsou měřená data zatížená chybou ► tedy „Při daných parametrech au jaká je pravděpodobnost měření (xuyt)?" PA081: Programování numerických výpočtů A. Křenek Modelování experimentálních dat Metoda nejmenších čtverců Lineární regrese Obecný model Lineární modely SAXS Vícerozměrná data Metoda nejmenších absolutních odchylek Regularizace Metoda nej menších čtverců Odvození ► „Jaká je pravděpodobnost, že konkrétni sada parametru at je správná?" ► špatně položená otázka ► neexistuje „náhodná veličina modelů" ► naopak, náhodnou veličinou jsou měřená data zatížená chybou ► tedy „Při daných parametrech au jaká je pravděpodobnost měření (xuyt)?" ► nulová, je-li y spojitá veličina ► musíme přidat „plus/minus odchylka měření Ay" ► model považujeme za správný, maximalizuje-li tuto pravděpodobnost ► i tak je to velmi intuitivní konstrukce PA081: Programování numerických výpočtů A. Křenek Modelování experimentálních dat Metoda nejmenších čtverců Lineární regrese Obecný model Lineární modely SAXS Vícerozměrná data Metoda nejmenších absolutních odchylek Regularizace Metoda nejmenších čtverců Odvození PA081: Programování numerických výpočtů A. Křenek ► předpokládáme normální rozložení chyby měření ► pravděpodobnost výskytu dané sady měření n e 2 v a ) Ay ► maximalizace odpovídá minimalizaci logaritmu, tj Il)>-,y)-»My 2o-2 ► JV, O", jsou konstanty Modelování experimentálních dat Metoda nejmenších čtverců Lineární regrese Obecný model Lineární modely SAXS Vícerozměrná data Metoda nejmenších absolutních odchylek Regularizace Metoda nejmenších čtverců Poznámky ► rozložení chyby všech měření nemusí být stejné ► používá se modifikovaná funkce x2 = S (yí-M(Xj))2 2 a? PA081: Programování numerických výpočtů A. Křenek Modelování experimentálních dat Metoda nejmenších čtverců Lineární regrese Obecný model Lineární modely SAXS Vícerozměrná data Metoda nejmenších absolutních odchylek Regularizace Metoda nej menších čtverců Poznámky ► rozložení chyby všech měření nemusí být stejné ► používá se modifikovaná funkce 2 v (yi-M(Xi))2 ► rozložení chyby nemusí být normální ► počet měření v jednom bodě bývá příliš malý ► zatížení chybou typu „někdo kopl do váhy" ► metoda je na takové chyby nepřiměřeně citlivá ► tzv. robustní statistiky Programování numerických výpočtů A. Křenek Modelování experimentálních dat Metoda nejmenších čtverců Lineární regrese Obecný model Lineární modely SAXS Vícerozměrná data Metoda nejmenších absolutních odchylek Regularizace Metoda nej menších čtverců Poznámky ► rozložení chyby všech měření nemusí být stejné ► používá se modifikovaná funkce 2 v (yi-M(Xi))2 X - ^ 2a\ ► rozložení chyby nemusí být normální ► počet měření v jednom bodě bývá příliš malý ► zatížení chybou typu „někdo kopl do váhy" ► metoda je na takové chyby nepřiměřeně citlivá ► tzv. robustní statistiky ► systematická chyba ► např. špatně kalibrovaný přístroj, závislost na související veličme Programování numerických výpočtů A. Křenek Modelování experimentálních dat Metoda nejmenších čtverců Lineární regrese Obecný model Lineární modely SAXS Vícerozměrná data Metoda nejmenších absolutních odchylek Regularizace Metoda nejmenších čtverců Zhodnocení vypočtených parametrů ► minimalizací x2 se téměř vždy hodnot at dopočítáme ► nevypovídá to ještě nic o kvalitě modelu ► např. lineární model radioaktivního rozpadu 25 20 15 10 0 10 20 30 40 50 60 70 80 90 100 PA081: Programování numerických výpočtů A. Křenek Modelování experimentálních dat Metoda nejmenších čtverců 1 1 i i i i i i i "exp.dat" + 13.849155 -0.168764*x ------- regrese + Obecný model Lineární + + + + modely SAXS Vícerozměrná + ^ data Metoda i i + ++ # + + + + ++++>-++ ++ + i i i i i i i nejmenších absolutních odchylek Regularizace Metoda nejmenších čtverců Zhodnocení vypočtených parametrů ► „chi by eye" nebo seriozní statistické zhodnocení ► hodnotíme pomocí regularizované gamma funkce ► funkce konkrétního x2 a počtu stupňů volnosti (N -M) Q N-M x' ' 2 ► viz např. http://en.wikipedia.org/wiki/Incomplete. gamma_function ► pravděpodobnost, že zcela náhodně vybraný vzorek (Xí,yí) dá větší hodnotu x2 ► čím větší tím lepší ► Q > 0.1 je v pořádku ► Q e [0.001,0.1] je podezřelé, ale stále přijatelné, není-li distribuce chyb měření zcela normální, resp. je mírně podceněná ► Q < 0.001 znamená špatný model nebo zcela nesmyslné měření PA081: Programování numerických výpočtů A. Křenek Modelování experimentálních dat Metoda nejmenších čtverců Lineární regrese Obecný model Lineární modely SAXS Vícerozměrná data Metoda nejmenších absolutních odchylek Regularizace Lineární regrese ► data prokládáme přímkou a + bx = 0 ► obecnější než se zdá na první pohled ► data (x,y) lze předem libovolně (nelineárně) transformovat na (x',y') ► minimalizovaná funkce xHa,b) = X-- 2a- ► minimum v bodě nulových prvních derivací 0 = 0 = 3xi da db = -2l = -2Z y; - a - fox; - a ~ bxi) PA081: Programování numerických výpočtů A. Křenek Modelování experimentálních dat Metoda nejmenších čtverců Lineární regrese Obecný model Lineární modely SAXS Vícerozměrná data Metoda nejmenších absolutních odchylek Regularizace Lineární regrese ► vhodným vyjádřením faktorů S, Sx,Sy,Sxx, Sxy ► součty zlomků konstruovaných z %u Ju o~i ► vše jsou to tady konstanty ► řešíme systém lineárních rovnic Sa + Sxb = 5 Sxa + Sxxb = S XX >xy ► získáme řešení, ale nevíme nic o jeho kvalitě ► odhad podle grafu nebo výpočet Q ► uvedená lineární regrese na radioaktivní rozpad začíná být přijatelná až když připustíme 1.6 PA081: Programování numerických výpočtů A. Křenek Modelování experimentálních dat Metoda nejmenších čtverců Lineární regrese Obecný model Lineární modely SAXS Vícerozměrná data Metoda nejmenších absolutních odchylek Regularizace Obecný model ► M(x, ai.....ayi) je nelineární funkce v at ► nelinearita v x by nevadila, viz dále ► výpočet minimálního x standardními optimalizačními metodami ► existují speciální varianty právě pro tvar funkce X2(^i, ■ ■ ■, um) = X (yt -M(Xi,ai,...,aM)) 2a} ► včetně verzí s dostupnými prvními i druhými derivacemi ► díky speciálnímu použití další triky ► konkrétní metody ► Levenberg-Marquardt ► Moré ► např. NAG library PA081: Programování numerických výpočtů A. Křenek Modelování experimentálních dat Metoda nejmenších čtverců Lineární regrese Obecný model Lineární modely SAXS Vícerozměrná data Metoda nejmenších absolutních odchylek Regularizace 11/34 Lineární modely ► model M(x, ai.....ayi) je lineární kombinace M(x, ai,..., Um) = ^ajXj(x) ► linearita ve smyslu parametrů modelu a j ► základní funkce Xj mohou být jakékoli ► pro vyhodnocení modelu se použijí jen jejich konkrétní hodnoty v bodech xi ► opět minimalizujeme x2 ► definujeme matici A a vektor b Xj(Xi) (Ti (Ti PA081: Programování numerických výpočtů A. Křenek Modelování experimentálních dat Metoda nejmenších čtverců Lineární regrese Obecný model Lineární modely SAXS Vícerozměrná data Metoda nejmenších absolutních odchylek Regularizace Lineární modely ► derivováním dostáváme M rovnic pro k = 1,..., M 0 = Z ~2 \ Vi ~ ZajXjixt)] Xk(Xi) j ► po úpravách v maticovém vyjádření (ArA)a = A1 b ► řešení bývá citlivé na zaokrouhlovací chyby ► preferovaná technika je QR rozklad Ti PA081: Programování numerických výpočtů A. Křenek Modelování experimentálních dat Metoda nejmenších čtverců Lineární regrese Obecný model Lineární modely SAXS Vícerozměrná data Metoda nejmenších absolutních odchylek Regularizace 13/34 Lineární modely Rozklad na singulární hodnoty ► model nemusí být dokonalý, mohou se objevit téměř lineární závislosti ► některé základní funkce nebo jejich kombinace jsou pro danou datovou sadu irelevantní ► vede na téměř singulární matici ► standardní metody inklinují k velkým hodnotám irelevantních parametrů ► SVD dokáže tyto problémy identifikovat ► algoritmus přímo hledá nejbližší řešení, tj. minimalizuje |Aa-b|2 ► zároveň detekuje problematické funkce PA081: Programování numerických výpočtů A. Křenek Modelování experimentálních dat Metoda nejmenších čtverců Lineární regrese Obecný model Lineární modely SAXS Vícerozměrná data Metoda nejmenších absolutních odchylek Regularizace Lineární modely Rozklad na singulární hodnoty ► stejná data, kvadratické a kubická funkce 25 20 15 10 "exp.dat" + 18.100248-0.426257*x+0.002662*x*x 19.418850-0.593832*x+0.007043*x*x-0.000031*x*x*x ...... + -K --,---Ť" 10 20 30 40 50 60 70 80 90 100 PA081: Programování numerických výpočtů A. Křenek Modelování experimentálních dat Metoda nejmenších čtverců Lineární regrese Obecný model Lineární modely SAXS Vícerozměrná data Metoda nejmenších absolutních odchylek Regularizace Lineární modely ► tabulka singulárních hodnot pro různé řády polynomu PA081: Programování numerických výpočtů A. Křenek Modelování experimentálních dat Metoda nejmenších čtverců řád (Ti Lineární 1 5.18 537.24 regrese 2 39654.25 3.51 134.10 Obecný model 3 31922200.00 66338.90 564.94 25.73 Lineární modely 4 2685350000.00 4133810.00 19913.70 272.78 99.50 saxs ► koeficienty u xn pro n > 4 jsou téměř nulové Vícerozměrná data Metoda nejmenších absolutních odchylek Regularizace 16/34 Komplexní příklad - SAXS ► Small Angle X-ray Scattering ► analytická metoda přinášející informace o struktuře molekuly ► ve vodním roztoku - přirozené prostředí pro biochemii ► relativně levný a nenáročný experiment ► RTG paprsek se v kontaktu s atomem rozptýlí ► je-li v dráze odchýleného paprsku další atom, dochází k interferenci ► detektor měří intenzitu záření v celém rozsahu rozptylu ► molekuly jsou v roztoku v náhodné orientaci ► měřená křivka je průměrem všech možných orientací PA081: Programování numerických výpočtů A. Křenek Modelování experimentálních dat Metoda nejmenších čtverců Lineární regrese Obecný model Lineární modely SAXS Vícerozměrná data Metoda nejmenších absolutních odchylek Regularizace Experimentální data ► sada obrázků ► rozlišení -20 Mpix ► monochromatické - jen intenzita ► křivka SAXS ► závislost log / na q = 4n sin0/A (9 - úhel odchýlení, A - vlnová délka) ► získána průměrem v kruzích ► cca. 2000 bodů (hodnot q) ► inforamace o chybě ► spolehlivá, průměr ováných pixelů je mnoho < □ ► < igi ► < Experimentální data i i 'pokus2.dat' u 1:2 'pokus2.dat' u 1:3 0.01 0.1 0.2 0.3 0.4 0.5 0.6 0.7 PA081: Programování numerických výpočtů A. Křenek Modelování experimentálních dat Metoda nejmenších čtverců Lineární regrese Obecný model Lineární modely SAXS Vícerozměrná data Metoda nejmenších absolutních odchylek Regularizace Model ► předpokládaná strukutura molekuly ► 3D model ► typy atomů ► jejich polohy => vzájemné vzdálenosti ► dodatečné parametry ► chování solventu (vody) ► anomálie na povrchu molekuly ► voda se zahušťuje nebo ředí podle povrchových vlastností ► měřítko křivky ► jde o tvar, absolutní hodnota nemá přímý smysl ► funkce cI(q,ci,C2) PA081: Programování numerických výpočtů A. Křenek Modelování experimentálních dat Metoda nejmenších čtverců Lineární regrese Obecný model Lineární modely SAXS Vícerozměrná data Metoda nejmenších absolutních odchylek Regularizace 20/34 Regrese modelu ► lineární model pouze v c, nelineární v C\, c2 ► minimalizujeme funkci x2 = -l \—■ / ^exp (at) - a (au ci, c2) 2(JV-1) J V ctí ► v rozumném rozsahu c, C\,c2 pouze 1 lokální minimum ► bylo by možné přímo optimalizovat ► neděláme kvůli složitějšímu případu (viz dále) PA08i: Programování numerických výpočtů A. Křenek Modelování experimentálních dat Metoda nejmenších čtverců Lineární regrese Obecný model Lineární modely SAXS Vícerozměrná data Metoda nejmenších absolutních odchylek Regularizace Regrese modelu ► lineární model pouze v c, nelineární v c±, c2 ► minimalizujeme funkci x2 = 1 ^ flexp((li) -cl(quci,c2) 2(JV-1) J V (Ji ► v rozumném rozsahu c, c±,c2 pouze 1 lokální minimum ► bylo by možné přímo optimalizovat ► neděláme kvůli složitějšímu případu (viz dále) ► zafixujeme c\,c2 a derivujeme dc 1 y 2(Jexp(^i) - cl(quci,c2)) (-I(cn,ci,c2)) 2(N-1) f-t i=l PA081: Programování numerických výpočtů A. Křenek Modelování experimentálních dat Metoda nejmenších čtverců Lineární regrese Obecný model Lineární modely SAXS Vícerozměrná data Metoda nejmenších absolutních odchylek Regularizace 21/34 Regrese modelu ► optimální hodnota N Jexp c = N I(quci,c2)2 ► počítáme opakovaně pro razné hodnoty C\,C2 ► držíme se smysluplných intervalů ► stačí 20 x 80 výpočtů ► ze všech vybereme tu s nejmenším x2 ► celkem 1600AT vyhodnocení I(q3ci,C2) ► složitost I(q, Ci, C2) je O (M2) v počtu atomů PA081: Programování numerických výpočtů A. Křenek Modelování experimentálních dat Metoda nejmenších čtverců Lineární regrese Obecný model Lineární modely SAXS Vícerozměrná data Metoda nejmenších absolutních odchylek Regularizace Regrese modelu Foxs chi = 1.07 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 FOXS Chi = 1.04 O 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 PA081: Programování numerických výpočtů A. Křenek Modelování experimentálních dat Metoda nejmenších čtverců Lineární regrese Obecný model Lineární modely SAXS Vícerozměrná data Metoda nejmenších absolutních odchylek Regularizace <<□► 4 ^ >■ 4 1 ► <=► 23/34 Flexibilní molekuly ► molekula může existovat ve více tvarech ► spontánně přechází mezi nimi ► v roztoku existují všechny současně ► SAXS křivka je jejich průměrem ► počítáme s více modely //(g, Ci, £2), J = 1, ■ ■ ■ ,K ► každému přiřadíme váhu Wj a kombinujeme je lineárně ► omezení Xw, = 1 PA081: Programování numerických výpočtů A. Křenek Modelování experimentálních dat Metoda nejmenších čtverců Lineární regrese Obecný model Lineární modely SAXS Vícerozměrná data Metoda nejmenších absolutních odchylek Regularizace 24/34 Flexibilní molekuly ► různá citlivost na změny c\, cz 15 20 0 10 20 30 40 50 60 70 8Q ► vzniká mnoho lokálních minim PA081: Programování numerických výpočtů A. Křenek Modelování experimentálních dat Metoda nejmenších čtverců Lineární regrese Obecný model Lineární modely SAXS Vícerozměrná data Metoda nejmenších absolutních odchylek Regularizace 25/34 Flexibilní molekuly ► opět iterujeme různé hodnoty ci, C2 ► zavedeme substituci Wj = cwj ► parametr c následně zrekonstuujeme c = ► můžeme si dovolit - regrese netáhne k falešnému řešení ► řešíme globální optimalizací pro 01^2,1^1.....Wk s omezením 0 < wl < 1 PA08i: Programování numerických výpočtů A. Křenek Modelování experimentálních dat Metoda nejmenších čtverců Lineární regrese Obecný model Lineární modely SAXS Vícerozměrná data Metoda nejmenších absolutních odchylek Regularizace Flexibilní molekuly ► opět iterujeme různé hodnoty ci, C2 ► zavedeme substituci Wj = cwj ► parametr c následně zrekonstuujeme c = ► můžeme si dovolit - regrese netáhne k falešnému řešení ► řešíme globální optimalizací pro 01^2,1^1.....Wk s omezením 0 < wl < 1 ► lze čekat, že některá Wj budou irelevantní ► Powellova metoda BOBYQA ► startujeme z více bodů PA08i: Programování numerických výpočtů A. Křenek Modelování experimentálních dat Metoda nejmenších čtverců Lineární regrese Obecný model Lineární modely SAXS Vícerozměrná data Metoda nejmenších absolutních odchylek Regularizace Flexibilní molekuly PA081: Programování numerických výpočtů A. Křenek Modelování experimentálních dat Metoda nejmenších čtverců Lineární regrese Obecný model Lineární modely SAXS Vícerozměrná data Metoda nejmenších absolutních odchylek Regularizace 27/34 Vícerozměrná data ► místo dvojic (Xuji) máme (xi,yt) ► x je fe-rozměrný vektor ► model M je funkce Rk — R ► jinak se nic nemění ► minimalizujeme vůči parametrům ► základní funkce se pouze vyhodnocují v ► není třeba derivovat podle složek x* PA081: Programování numerických výpočtů A. Křenek Modelování experimentálních dat Metoda nejmenších čtverců Lineární regrese Obecný model Lineární modely SAXS Vícerozměrná data Metoda nejmenších absolutních odchylek Regularizace Vícerozměrná data ► místo dvojic (xuyi) máme (xi,yt) ► x je fe-rozměrný vektor ► model M je funkce Rk — R ► jinak se nic nemění ► minimalizujeme vůči parametrům ► základní funkce se pouze vyhodnocují v ► není třeba derivovat podle složek x* ► obecněji M: Rk - Rl ► potřebujeme vhodnou metriku na Kl ► počítáme se vzdáleností \yit M(x*) | podle ní Metoda nejmenších absolutních odchylek ► metoda nejmenších čtverců minimalizovala YJ(yi-M(xuall...1aM))2 tj. normu řádu 2 ► metoda nejmenších absolutních odchylek minimalizuje normu řádu 1 X \yt -M(xuai,...,aM) PA081: Programování numerických výpočtů A. Křenek Modelování experimentálních dat Metoda nejmenších čtverců Lineární regrese Obecný model Lineární modely SAXS Vícerozměrná data Metoda nejmenších absolutních odchylek Regularizace 29/34 Metoda nejmenších absolutních odchylek Výhody a nevýhody ► větší robustnost ► méně citlivá na odlehlé případy ► nemají kvadratickou váhu, snáz se „přebijí" ► menší stabilita ► malý posun v x může mít velký vliv na výsledné řešení ► nejednoznačné řešení ► lineární členy se vzájemně kompenzují proti posunu v ose y ► interaktivní srovnání na lineární regresi: ► http://www.math.wpi.edu/Course_Materials/SAS/ 1ablets/7.3/73_choi ces.html PA081: Programování numerických výpočtů A. Křenek Modelování experimentálních dat Metoda nejmenších čtverců Lineární regrese Obecný model Lineární modely SAXS Vícerozměrná data Metoda nejmenších absolutních odchylek Regularizace Metoda nejmenších absolutních odchylek Výpočet ► nejmenší čtverce ► derivace kvadratické funkce ► vede na systém lineárních rovnic ► absolutní odchylky ► úloha lineárního programování ► např. Barrodale-Robertsův algoritmus PA081: Programování numerických výpočtů A. Křenek Modelování experimentálních dat Metoda nejmenších čtverců Lineární regrese Obecný model Lineární modely SAXS Vícerozměrná data Metoda nejmenších absolutních odchylek Regularizace Regularizace ► základní regresní metody fungují dobře pro ideální případy ► mohou selhávat na reálných datech ► snaží se modelovat šum více než vlastní chování systému PA081: Programování numerických výpočtů A. Křenek Modelování experimentálních dat Metoda nejmenších čtverců Lineární regrese Obecný model Lineární modely SAXS Vícerozměrná data Metoda nejmenších absolutních odchylek Regularizace 32/34 Regularizace ► místo minimalizace výrazu ||Aa - b||2 minimalizujeme ||Aa-b||2 + IITall2 s vhodně volenou maticí ľ (Tikhonovova regularizace) ► podobně jako u nejmenších čtverců je řešením (ÁrA-ľ1ľ)-1A1b ► volbou T zvýhodňujeme nějaké řešení ► např. ľ = al preferuje menší normu T 1 aTi PA081: Programování numerických výpočtů A. Křenek Modelování experimentálních dat Metoda nejmenších čtverců Lineární regrese Obecný model Lineární modely SAXS Vícerozměrná data Metoda nejmenších absolutních odchylek Regularizace Regularizace Vztah k SVD ► rozklad A na signgulární hodnoty A = UXVT ► potom řešení regularizováného problému s r = al je VI/ UTb kde prvky Z jsou (Ji tj. a významněji ovlivní právě „skoro nulové" singulární hodnoty PA081: Programování numerických výpočtů A. Křenek Modelování experimentálních dat Metoda nejmenších čtverců Lineární regrese Obecný model Lineární modely SAXS Vícerozměrná data Metoda nejmenších absolutních odchylek Regularizace 34/34