PA081: Programování numerických výpočtů
10. Modelování experimentálních dat, metoda nejmenších čtverců
Aleš Křenek
jaro 2011
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ů ai, tj.
y = M(x,ai,...,aM)
► hledáme takové hodnoty au pro něž model odpovíd nejlépe experimentu
Modelování experimentálních dat
Příklad
t\n2
► radioaktivní rozpad N = Noe t
► N je počet atomů ve vzorku (JV0 v čase ř = 0), T poločas rozpadu
+		" príklady/ex p. dat" + 20*exp(-x/20*log(2)) -------
A +		-
		-
-	+ +"+*--	
0 10 20 30 40 50 60 70 80 90 100
3/17
Metoda nejmenších čtverců
Odvození
PA081: Programování numerických
výpočtů
A. Křenek
► „Jaká je pravděpodobnost, že konkrétní sada parametrů cii 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 ai, jaká je pravděpodobnost měření (xuyt)?"
Metoda nejmenších čtverců
Odvození
PA081: Programování numerických
výpočtů
A. Křenek
► „Jaká je pravděpodobnost, že konkrétní sada parametrů cii 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 ai, 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
Metoda nejmenších čtverců
Odvození
► předpokládáme normální rozložení chyby měření
► pravděpodobnost výskytu dané sady měření
Y\e H   -   ) Ay
► maximalizace odpovídá minimalizaci logaritmu, tj.
► N, (J, Ay jsou konstanty
Metoda nejmenších čtverců
Poznámky
► rozložení chyby všech měření nemusí být stejné ► používá se modifikovaná funkce
2      v (Jj -M(Xj))2
Metoda nejmenších čtverců
Poznámky
► rozložení chyby všech měření nemusí být stejné
► používá se modifikovaná funkce
*■ 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
Metoda nejmenších čtverců
Poznámky
► rozložení chyby všech měření nemusí být stejné ► používá se modifikovaná funkce
= 1
(yí-M(Xí))2 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čině
□     s   4 :
■OQ.O
Metoda nejmenších čtverců
Zhodnocení vypočtených parametrů
► minimalizací x2 se téměř vždy hodnot ai dopočítáme
► nevypovídá to ještě nic o kvalitě modelu
► např. lineární model radioaktivního rozpadu
0 10 20 30 40 50 60 70 80 90 100
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)
► viz např. http://en.wikipedia.org/wiki/Incomplete gamma_function
► pravděpodobnost, že zcela náhodně vybraný vzorek (x;, yi) dá větší hodnotu x2
čím větší tím lepší
► Q > 0.1 je v pořádku
► Q G [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í
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)^{yi-a2^hXl)2 minimum v bodě nulových prvních derivací
o = &- = -2Yyt-a7bXt
Ba cr2
db      ^ o-,2
Lineární regrese
vhodným vyjádřením faktorů 5, Sx, Sy, Sxx, Sxy
► součty zlomků konstruovaných zxj,Jí,(7í
► vše jsou to tady konstanty
řešíme systém lineárních rovnic
Sa + Sxb = Sy
Sxa + Sxxb = S,
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 07 > 1.6
Obecný model
M(x, a\,au) je nelineární funkce v cii ► 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
(yt -M(Xi,ai,...,aM))z
X2(a1,...,aM) = X'
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
Lineární modely
► model M(x, a\,..., au) je lineární kombinace
M(x, ai,...,aM) = YJajxjM
► 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 x;
► opět minimalizujeme x2
► definujeme matici A a vektor b
Lineární modely
► derivováním dostáváme M rovnic pro k = l,...,M
o = X^2 [yt - X"jXAxi)\ xk(Xi)
*■ po úpravách v maticovém vyjádření
(ArA)a = Arb
► řešení bývá citlivé na zaokrouhlovací chyby
► preferovaná technika je QR rozklad
Lineární modely
Rozklad na singulární hodnoty
PA081: Programování numerických
výpočtů
A. Křenek
► 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
■O0.O
Lineární modely
Rozklad na singulární hodnoty
► stejná data, kvadratické a kubická funkce
"exp.dal"
18.100248-0.426257"x+0.002662*x*x 19.418850-0.593832*x+0.007043-x*x-0.000031 *x*x*x
_1_I_I_I-1-1-1-1-1-
10 20 30 40 50 60 70 80 90 100
15/17
PA081:
Lineární modely
numerických výpočtů
A. Křenek
► tabulka singulárních hodnot pro různé řády polynomu
řád	07				
1	5.18	537	24		
2	39654.25	3	51	134.10	
3	31922200.00	66338	90	564.94	25.73
4	2685350000.00	4133810	00	19913.70	272.78 99.50
► koeficienty u xn pro n > 4 jsou téměř nulové
Vícerozměrná data
PA081: Programování numerických
výpočtů
A. Křenek
► místo dvojic (Xi,yi) máme (xuyi)
► x je k-rozměrný vektor
► model M je funkce Kfc -* K
► jinak se nic nemění
► minimalizujeme vůči parametrům
► základní funkce se pouze vyhodnocují v X;
► není třeba derivovat podle složek X;
■O0.O