PLÁNOVANIE REGRESNÉHO EXPERIMENTU Predložený text bol spracovaný hlavne podlá [5], [6], [7]. 1. ZÁKLADNÉ pojmy Prvá fáza prípravy experimentu spočíva v stanovení cielových parametrov a pokial tieto nie su priamo observovatelné (meratelné), v stanovení dostatočného počtu takých priamo observovatelných parametrov, ktoré sú vo vhodnom (vysvetlíme neskôr) známom funkčnom vztahu s cielovými parametrami. Vektor cielových parametrov označme (3. Budeme předpokládat, že (3 G lZk. Priamo observovatelné (meratelné) teoreticky bezchybné parametre (veličiny) označme /i2, /xjv0, teda počet priamo meratelných veličín je Nq. Ďalej predpokladáme, že poznáme funkciu f(-) : TZk ->• KN° (vyjadruj úcu meratelné parametre ako funkcie cielových parametrov). Túto situáciu popisuje teoretický model merania fi = f (P) = ( A(/3) \ h(P) \fN0((3)J Príklad 1.1. Úlohou je stanovit súradnice /?i,/?2 bodu A = (/?i,/?2), ked meriame vzdialenosti /ii,/i2,M3 (daných) bodov B = (xi,yi), C = (002,1/2) D = (x3,y3) od bodu A. Teoretický model merania je ^(h-X2Y + (fc-V2Ý V(ä-X3)2 + (/32-y3)2 Príklad 1.2. Treba určit váhu troch predmetov /3i,/?2,/?3 na ručičkových váhach (majú jednu misku). Vážit môžeme veličiny /ii, ...,/is, pričom teoretický model váženia je m2 a) + a M3 a) + P2 fi4 = a> + /33 M5 Po + P\ + h M6 /30 + /?1 + /?3 \/i8J Vä + ä + ä + ä/ kde napr. /íq znamená, že vážime prázdne váhy, /X5 znamená, že vážime spolu prvé a druhé závažie, atd. Stretli sme sa tu aj s novým fenoménom, a síce parametrom 1 2 pV Je to nulový údaj váh (prázdne váhy). Voláme ho tzv. rušivým parametrom. Z našich úvah ho vylúčime, hoci z modelu merania ho vy lučit nemôžeme. V nasledujúcom budeme předpokládat, že teoretický model merania je lineárny (alebo linearizovaný) v cielových parametroch, t.j. (1.1) fi = f (P) = F/3 = V r / kde F je známa N0 x k matica plánu, je jej z—ty riadok. Model z príkladu 1.2 je lineárny. Model z príkladu 1.1 vieme linearizovat. Ako ? Príklad 1.1 - pokračovanie. Nech = 0, 00m, yi = 0, 00m; x2 = 2365, 22m, í)2 = 0,00m; x3 = 3603,67m, y% = 823,35m. Približné (namerané) hodnoty IH « 1980,102m, h2 ~ 2040, 243m, /z3 « 2598, 878m. Z rovníc 1980,102 = y/(p10 - 0,00)2 + (/?20 - 0,00)2 2040, 243 = V(Äo - 2365,22)2 + (/320 - 0,00)2 spočítame hodnoty p\o a /320, z čoho dostávame /?io = 1131,5m a /?2o = 1625, 0m. Funkciu f linearizujeme (rozvinieme do Taylorovho radu a zanedbáme členy rádu druhého a vyšších) okolo hodnôt /?io a /?2o, teda V(Pio - ^) + (P20 - Vi)2 \/(/3io - ^)2 + (äo - y*)2 Dostávame linearizovaný model (1.1) Mi- 1980,131 \ / 0, 571<% + 0,8215/32 \ / 0, 571 0, 821 \ /z2 - 2040,267 = -0, 6055/?i + 0,796<% = 0, 605 0, 796 m - 2598,897/ \-0, 951<% + 0,308<% / \-0,951 0,308/ s (novými) parametrami 5j3\ a ôfo. Observovatelný parameter /Zj samozrejme nepoznáme presne. Výsledok jeho zmerania je číslo yi, ktoré považujeme za realizáciu náhodnej veličiny Yi. Váha tohto merania je Aj a je nepriamo úmerná disperzii 2?(YÍ). Všetky merania považujeme v nasledujúcom za neskorelované. Dostávame sa k stochastickému modelu merania. V prípade, že práve jedenkrát (nezávisle) meriame každý priamo observovatelný parameter /Zj, i = 1,2,...,/Vo a váhy jednotlivých meraní sú Aj, í = 1,2,...,/Vo, stochastický model merania je lineárny regresný model (Yjv0,i, Fjv0,fe/3fe,i, cr2A_1). 3 Observačný vektor (vektor meraní) Y má vektor stredných hodnôt F/3 a kova-riančnú maticu a2 Ar1, kde A /Ai 0 ... 0 \ ' 0 A2 ... 0 ' V 0 0 AjV„ / Ďalšia fáza prípravy experimentu spočíva v respektovaní pravidiel optimálneho návrhu (regresného) experimentu, kedy odpovedáme na otázku, kolkokrát ktorý priamo observovatelný parameter treba merat, aby výsledok spracovania respektoval určité dopredu zadané kritéria optimality, napr. aby experiment pri zadanej presnosti cielových parametrov bol čo najlacnejší, aby určité vybrané parametre boli odhadnuté s najväčšou možnou presnostou, t.j. s minimálnou možnou disperziou ich odhadov, atd. Príklad 1.3. Každý priamo observovatelný parameter meriame práve jedenkrát. Lineárny regresný model merania je (Yjv0ji, Fjv0,fe/3fe,i, c2A 1), teda stredná hodnota observačného vektora je V f 7 13 = F/3 a jeho kovariančná matica je a2 A 1. Najlepším lineárnym nevychýleným odhadom (NNLO) parametra (3 je /3(Y) = (F'AF)_1F'AY. Kovariančná matica NNLO $ je T,0 = cr2(F'AF)_1 = a2 < (fi,fjv0) /Ai 0 ' 0 A2 VO 0 AjVn / V f 7 Predpokladáme, že experiment je navrhnutý tak, že matica F'AF je regulárna (váhy Aj, i = 1, 2,No sú kladné a hodnost matice F, teda h(F) = k ^ No). Príklad 1.4. Vyberme r rôznych priamo observovatelných veličin /z^, Ve-ličinu híj merajme krát, pričom všetkých meraní nech je opát No, čiže ELi^ = No- 4 Observačný vektor a matica plánu v tomto prípade sú ■ N0,l — ' N0,k Y Y2,i Y 1 í2-™2 (riadok f4' je práve krát, j = 1,2, ...,r). Matica váh observačného vektora je * A = diaglX^,Aj1, \i2,Ai2,\ir,A^}, pričom na diagonále je Aj. práve nj krát, j = 1, 2,r. Ľahko vidíme, že NNLO /3 je v tomto prípade /3(*Y) = (*F' *A *F)~1 *F' *A *Y = („F' *A „F)"1 „F' *A „Y = /3(*Y), kde f f Í2 f Í2 f íl f,' 0 0 «í2^Í2 ° \ (Y ' Y i. \0 0 ... nir\ij pričom Y i = ^- Ya=i Yij,t, í = 1, 2,r, Y. !,Yi 2, Y. ni sú nezávislé mera-nia veličiny /Zj . Kovariančná matica odhadu /3(*Y) je = cr2(,F' „A „F)"1 = = v2 {({tl,...,{lr) nllXll 0 0 nl2 Xl2 V 0 0 \ x J f Í2 > = = íj2 ^ "•A- '• Zadefinujme si niektoré základné pojmy. Definícia 1.5. Funkcia 5 : {l,2,..,iVo} -> <0,1 > pre ktorú platí ô (i) = 1 sa nazýva návrh, plán alebo projekt experimentu (design of experiment). Ak celkový počet meraní je N, potom číslo rii = N ô (i) udáva počet opakovaných meraní hodnoty /íj. Číslo ô (i) je relatívny počet replikácií (opakovaní) merania hodnoty /íj. 5 Definícia 1.6. Množinu indexov S p (ô) = {i : ô (i) > 0} nazývame spektrom (suportom, nosičom) návrhu S. Samozrejme S p (ô) c {1,2, ...,Nq}. Meráme tie veličiny /zj1,/zj2, ...,Hir spomedzi všetkých experimentálnych bodov (observovatelných veličín, parametrov), teda z množiny £ = {/zi, /íjv0}, ktorých indexom plán S priradil nenulovú hodnotu, čiže pre ktoré S (i j) > 0. Definícia 1.7. Matica kde f/ je daný vektor, pre ktorý platí í[(3 = /Zj, sa nazýva informačná matica experimentu pri návrhu S. Definícia 1.8. Majme £ = {/zi, /zjv0} (množinu priamo observovatelných parametrov) a návrh S. Celkový počet meraní je N. Cielové parametre sú (3 G lZk. Poznáme ii - vektor, pre ktorý platí i^fS = /Zj a Aj - váhu merania veličiny /Zj, i = 1,2,Nq. Rešpektujeme plán S, t.j. opakujeme rii = N ô (i) krát meranie veličiny Hi, ak prirodzené číslo i G Sp(ô). Potom lineárny regresný model experimentu s (presným) návrhom S je (1.2) (Yä.FäA^A^1). Ak {ii, ...,ír} = Sp(S), tak matica F$ je vytvorená riadkami f'[ , j = 1,2, ...,r a A$ je diagonálna matica /Annn 0 ... 0 \ /KS(ii) 0 ... 0 \ ' 0 \l2nl2 ... 0 0 Aia5(i2) ... 0 » Vo 0 ... Avnv/ VO 0 ... \lr5{ir)J Yg je r rozmerný náhodný vektor, ktorého j - ta súradnica je {Y g} j = ^-(Yi i + Ytj,2 + - + Yh^), j = 1,2,...,r. Veta 1.9. Majme lineárny regresný model (1.2) z definície 1.8. Platí 1. Lineárny funkcionál g(-) vektora parametrov (3 je (lineárne a nevychýlené) odhadnutelný práve vtedy ak g {(S) = g'/3, pričom g G M(M(Ô)) = {M(S)u : u G Kk}. 2. Vektor parametrov (3 je (lineárne a nevychýlené) odhadnutelný práve vtedy ak M(t>) je regulárna matica. 3. Najlepší nevychýlený lineárny odhad (NNLO) parametra (3 je P(YS) = (FjAíFír^AíYí. 4- Kovariančná matica NNLO (3(Ys) je 2 cov(f3(Ys),N) =v2(F'sAsFs)-1 = ^-M_1(5). 6 Dôkaz. Spravte ako cvičenie. Pre experiment s cielovými parametrami Pi,..., a observovatelnými parametrami ni,/íjv0 možno určit tolko návrhov, kolko je funkcií S : {í, 2,Nq} —> < 0,1 > spĺňajúcich podmienku X)i=i^(*) = 1- -Pre = 2 je týchto funkcií nekonečne vela. Triedu všetkých návrhov označme A. Návrh t) e A nazveme regulárnym, ak det(M.(S)) =/= 0. Triedu všetkých regulárnych návrhov označíme Areg. Príklad 1.10. (podlá [7], str. 13) Majme tri predmety Ai,A2,A3, ktorých hmotnosti sú /3i,/?2,/?3 (neznáme). Pomocou N = 4 vážení treba určit (odhadnut), na ručičkových váhach tieto hmotnosti. Vytvorme si teoretický model váženia podlá príkladu 1.2 a stochastický model podlá príkladu 1.3. Nq = 8, £ = {fii, ■ 8,1 (Yl\ 8,4 f1 0 0 1 i 0 0 1 0 1 0 1 0 0 1 1 i 1 0 1 i 0 1 1 0 1 1 Vi i 1 1/ L8,8- I. organizácia váženia (nazvime ju bežnou): 1. váženie - zistenie nulovej výchylky váh n 2. váženie - váženie predmetu Ai 3. váženie - váženie predmetu A2 4. váženie - váženie predmetu A% Podlá príkladu 1.2 sme si vybrali 4 observatelné veličiny, a síce ni, m2; M3j M4- Počet všetkých vážení je N = 4, plán tohto (bežného) experimentu je Sb, pre ktorý platí ôb(í) = ôb(2) = 56(3) = ôb(4) = I 56(5) = 56(6) = ôb(7) = ôb(8) = 0, Sp(6b) = {1,2,3,4}. stredná hodnota je Observačný vektor je YSb = (Yia,Y"2,i,Y3A,>4,i)'- FsbJ = teda Jeho pričom 7 = (k, Pi, p2, P3)'■ Matica váh meraní je A^ = i44. NNLO parametra 7 je (F'SbAsbFgb)^1F'Sb AsbYgb a kovariančná matica tohto odhadu je ^(F'SbAsbFSb)-1 = -M-1(Sb) = ^ II. organizácia váženia (nazvime ju premyslenou): 1. váženie - váženie všetkých troch predmetov 2. váženie - váženie predmetu Ai 7 3. váženie - váženie predmetu A2 4. váženie - váženie predmetu A3 Tentokrát sme si vybrali 4 observatelné veličiny /i8, /i2, /i 3, /X4. Počet všetkých vážení je opát N = 4, plán (premysleného) experimentu je Sp, pre ktorý platí 5P(8) = Sp(2) = ôp(3) = Sp(4) = \, Sp(5) = Sp(6) = Sp(7) = Sp(í) = 0, teda Sp(óp) = {8,2,3,4}. Observačný vektor je Yáp = (y8,i>*2,i>*3,i>y4>1)'. Jeho stredná hodnota je Fáp7 = 7 = {n, P2, Ps)'■ matica váh meraní je opát A$p = I44. NNLO parametra 7 je (F^ AspFgp)^1F's AspYgp a kovariančná matica tohto odhadu je (T2(FLAápFáp)-1 = -M-1(5p) = (r2 1 -0,5 -0,5 -0,5 0,5 1 0 0 0,5 0 1 0 0,5 0 0 1 Ked porovnávame oba organizácie váženia (oba plány), vidíme, že pri oboch dostávame nevychýlené odhady neznámych hmotností P2, P3. Pri bežnom pláne majú tieto odhady disperzie 2a2, pokial pri premyslenom pláne sú disperzie odhadov a2, teda menšie. Navyše pri premyslenom pláne sú odhady neskorelované. 2. Najdôležitejšie kritériá optimality Návrh S treba vybrat tak, aby spĺňal nejaké kritéruim. Máme napr. určit čo najpresnejšie hodnotu h'/3, (3 e lZk (h je daný vektor). Merat môžeme N- krát. Za optimálny budeme považovat ten návrh S*, pre ktorý platí h'M^Onh = mMh'M-'fíjh : ô e Areg}. Návrh S* v tomto prípade minimalizuje disperziu odhadu h'/3. V praxi sa najčastejšie vyskytujú tie kriteriálne funkcie, ktoré sú uvedené v nasledujúcom. Používa sa pre ne označenie D—optimalita (podlá slova dispersion), A—optimalita (average), L—optimalita, reštringovaná A—optimalita a reštringova-ná D—optimalita. Iné kritériá pozri napr. v [7]. Definícia 2.1. Zobrazenie L(-) : Sk —> je lineárny pozitívny funkcionál definovaný na priestore Sk symetrických matíc typu k x k, pre ktorý platí (i)V{A,BeSk} L(A + B) = L(A) + L(B), (n) V{a e ft}V{A e Sk} L (a A) = aL(A), (iii) V{A G Sk : A je pozitívne definitná matica } L (A) > 0. Definícia 2.2. Návrh 5*D G Areg je D—optimálny ak deŕpVT1^)] = minidetlM.-1 (6)] : ó G Areg}. 8 Definícia 2.3. Návrh 5*A G Areg je A—optimálny ak Tr[M_1(^)] = mm{Tr[M_1(5)] : ó G Areg}. Definícia 2.4. Návrh 5*L G Ares je L—optimálny ak L\MT1(5*L)} = mm{L[M_1(5)] : 5 G Ares}. Nech je vektorový parameter (3 G 7?.fe vyjadrený v tvare či kde /3i (užitočný parameter) je k\ rozmerný a /32 (rušivý parameter) je kí rozmerný, pričom k\ + kí = k. V súlade s týmto rozkladom je rozložená aj informačná matica a jej inverzia. Teda platí, že pri použití plánu S a celkovom počte meraní N je (T2 kovariančná matica cov(0i) = —M.1'1 (ô), kde {lA) M ^ ' ~ \ M2'1(S) M2>2(5))- Definícia 2.5. Návrh 5*D G Areg je reštringovane D—optimálny ak detlM1'1 (ó^r)} = minidetlM1'1^)] : ó G Areg}. Definícia 2.6. Návrh 5*A G Areg je reštringovane A—optimálny ak TrlM1'1 (ô*Ar)] = mmjTrtM1'^)] : ó G Areg}. Trochu odlišné je kritérium S—optimality. Definícia 2.7. Návrh G Areg je S—optimálny ak T2 S--M-VSf) £- — M'1 (6) S e Areg,N = í,2, kde S je dopredu zadaná cielová kovariančná matica výsledného odhadu vektorového parametra (3, pričom ||A|| = ^/TV(AA'). V prípade S—optimality hladáme nielen plán 5^, ale aj optimálny počet meraní N*. Kritérium D—optimality má nasledovnú interpretáciu: Ak Ys ~ Nn(Fsf3, ct2AJ1), tak (1 — a)— konfidenčný elipsoid pre vektor j3 G 7čfe pri iV meraniach je £^a((3) = (u : (u - 0yN^fl{u -p) í x2(0; 1 - a) 9 a jeho objem je v{5)_ ^ ixim-<*)]' r(i + f)^eí[^(nAáFá)] (pozri [1], kapitola 11.12). Pretože 1 ^ VdetpVl-1^)], det[£ (FJAÄFÄ)] N' D—optimalita návrhu zaručuje minimálny objem konfidenčného elipsoidu. Pri použití tohoto kritéria je niekedy potrebné kontrolovat približnú gulatost konfidenčného elipsoidu. Príliš velké rozdiely medzi velkostami jeho hlavných poloosí môžu niekedy signalizovat nežiadúce vlastnosti návrhu. Na druhej strane D—optimalita má tzv. minimaxnú vlastnost (pozri vetu 3.1 v kapitole 3). Táto vlastnost návrhu S môže byt v niektorých prípadoch velmi dôležitá. Preto sa D—optimalita v praxi pomerne často používa. Pretože 2 cov[f3s] = ^rM-1^), A—optimálny plán minimalizuje súčet disperzií odhadov zložiek vektora (3. Kritérium A (A—optimalita) je špeciálnym prípadom L—optimality, lebo TV(-) je lineárny a pozitívny funkcionál. Pri riešení odhadu lineárnej funkcie h{(3) = h'/3 s minimálnou disperziou odhadu (pozri začiatok tejto kapitoly) ide zase o L—optimálny plán, lebo funkcionál L(-) definovaný vztahom L (A.) = h'Ah (h je daný pevný vektor) je opát pozitívny lineárny funkcionál. Špeciálnym prípadom L—optimality je aj reštringovaná A—optimalita, ked minimalizujeme TVpVl1'1^)] (pozri (1.2)) vhodnou volbou S G Areg. Matica M1'1 (5) patrí parametrom, pre ktoré chceme minimalizovat súčet disperzií ich odhadov. V prípade S—optimality ide o maximálne priblíženie (v danej norme) matice cov(J3) = — M_1(5) k cielovej matici S. Okrem uvedených kritérií sa objavujú kritériá motivované špeciálnymi požiadavkami užívatela. Obyčajne majú konvexnú (alebo konkávnu) kriteriálnu funkciu. Niekedy sa použije kritérium, ktoré je konvexnou kombináciou uvedených kritérií. Jedná sa o snahu udržat dobré vlastnosti oboch kritérií, alebo potlačit nežiadúcu vlastnost návrhu optimálneho podlá jedného kritéria priblížením k návrhu optimálneho podlá iného kritéria. Podrobná teória o kritériách optimality je napr. v [7], kde je uvedená aj bohatá literatúra o optimálnom navrhovaní regresného experimenta. 3. Vety o ekvivalencii pre niektoré kritériá optimality Veta 3.1. Nasledujúce tvrdenia sú ekvivalentné. 1. Návrh 5*D G Areg je D—optimálny, teda cfeípVT1^)] = minidetlM.-1 (ô)] : S G Areg}. 2. Návrh 5*D G Areg minimalizuje d(ô) = max{\iťiM.^1{5)íi : i = 1,2, ...,Nq}, teda d(ôjy) = min{d(5) : S G Areg}. 3. d(S*D) = k (dimenzia vektora (3). 10 Dôkaz. Dôkaz vety realizujeme nasledovným spôsobom: 1. => 2. a súčasne 2. •<=> 3. a 2. 1. 1. =>• 2. a súčasne 2. •<=> 3. Nech S*D G Ares je D-optimálny, teda deŕ[M_1 = mm{deŕ[M_1(5)] : S G Ares}. Vezmime lubovolný návrh S a a G< 0,1). Návrh S = (1 — + at) je podia lemy 8.11 regulárny, pričom M(6) = [(í-a)ô*D(i) + aô(i)]Xiíiíl = (l-a)M(ô*D) + aM(S). i£Sp(Š) Podlá lemy 8.12 je pre a G< 0, 1) funkcia g (a) = lndet~M.(S) spojitá diferencovatelná, pričom g(0) = max{g(a) : a G< 0,1)}. Teda g(0) ^ g(a) pre a G< 0,1) a musí byt d lim -^q(a) S 0 a^o da y ' - (ide o deriváciu sprava). Dostávame (pomocou lemy 8.6) lim —IndetUí - a)M(ó*n) + aM(ó)} = a^o da = lim Tr{[(l - a)M(ó*D) + aMfä)]"1 [-^-((1 - a)M(ó*D) + aM(ó))}} = a^>o da = Tr{M-1(S*D)[-M(S*D) +M(S)]} = -Trlk,k + TrM'1 (5*D)M(5) = = Tr[M-1(S*D) S(i)Xiíiíi] - k = SiWlM-^ô^íi-k^O, teda (3.i) J2 m^M-^s^^k. Ak vezmeme jednobodový návrh S so S p (ô) = {io}, tak z (3.1) pre íq = 1, 2,Nq je A^ťM-1^)^ S A;, (3.2) d(<^) =max{Aíf;M"1(^)fí : z = 1, 2,iV0} 5; k. Pravda pre každý regulárny návrh S G Areg platí N0 No (3.3) jfc = TrM-1(S)M(S) = TrpVT1^) ^ A^M] = '£iô(i)ÁitlM-1(ô)ti. i=i i=i Okrem toho pre každý S G Areg je (3.4) d(ô) = maxiXďJs/r1 {5% : í = 1, 2,N0} ^ k. 11 Dôkaz tvrdenia (3.4) vykonáme sporom. Ak by pre nejaký návrh r) G Areg bolo max{Xlí'iM^1(7f)íl : í = 1,2, ...,N0} < k, tak z (3.3) N0 N0 k = ^iti(Í)Xí{ÍM-1(ti){í S ^v^maxiXiíiM-1^ : i = 1, 2,JV0} = i=i i=i No = max{\£M.-1(r])tl : í = 1, 2, N0} ^ 7?(i) = i=l = maxíAjf/M"1^)^ : i = 1,2, ...,-/V0} < fc, čo je v spore s (3.3). Z (3.2) a (3.4) pre 5*D (pretože 5*D G Areg) dostávame k ^ d(^) ^ /j, čiže d(5p) = k = min{d(5) : S G Areg}. Dokázali sme 1. 2. a tiež 2. •<=> 3. 2. 1. Ak teda máme regulárny návrh S G Areg, tak = maxiXií^M-1 (ô)íi : i = 1,2, W0} ^ d(<%) = Ä:, kde je (lubovolný) D—optimálny návrh. Nech S G Ares minimalizuje d(ô) = max{Xlí'iMr1 {5% : i =1,2, ...,iV0} na množine Ares. Musí byt (3.5) = k, (lebo podlá (3.4) pre každý S G Ares je d(S) k a, pre (lubovolný) D—optimálny návrh 5*D je d{5*D) = k). Platí (pretože t)^ je D—optimálny) 0 < deŕM_1(^) S detM'1^), teda (podlá lemy 8.9) k (3.6) 1 ^ deťM(5*D)deťMT1 (5) = det^ô^MT1 (5)] = J|7í i=l (71,, 7fe sú vlastné hodnoty matice M(#£,)M_1(5), ktoré sú podlá lemy 8.8 reálne a kladné), čiže 1 (k \ k 12 Na druhej strane podia lemy 8.10, lemy 8.9 a (3.5) je i = -ttm(ô*d)m-i(ô) = -ttm-i(ô) s*Dmw = ieSP(S'D) iesP(S'D) = \ J2 S^maxiXitiM-1^: i = 1, 2, N0} = iesP(S'D) (3.8) = - = ieSp(s*D) z čoho dostávame k (3.9) n'k 1=1 Zo vztahov (3.6) a (3.9) máme 1 S detMió^detM.-1 (ó) ^ 1, teda detM(ô*D) = detM(ô), čo znamená, že deŕM(5) = max{det~M.(ii) : /i G Ares}, alebo ekvivalentne deťMT1(5) = minideťMT1 (/i) : /i G Ares}. Dokázali sme 2. =>• 1. celú vetu. □ Veta 3.2. Nasledujúce tvrdenia sú ekvivalentné. 1. Návrh 5*A G Areg je A—optimálny, teda Tr[M_1 (SA)] = min{Tr\MT1{5)] : S G Areg}. 2. Návrh 5*A G Areg minimalizuje A{5)= max {\iťi\$A~1 {ô)]2 íi : i =1,2, ...,N0}, teda A(S*A) = min{A(S) : S G Areg}. 3. A(ô*a) = TrM"1^). Dôkaz. Dôkaz vety realizujeme nasledovným spôsobom: 1. => 2., 1. => 3., 2. => í., 3. 1. 1. => 2. -kk y. %(o = i, i£Sp(S*D) = 1. 13 Nech Ô*A e Areg je A-optimálny, teda TrM_1(^) = mín{TrM.-1(5) : S e Areg}. Vezmime lubovolný návrh S a a G< 0, 1). Návrh S = (1 — + aS je podia lemy 8.11 regulárny, pričom M(6) = J2 [(í-a)ô*A(i) + aô(i)]Xiíiíl = (l-a)M(ô*A) + aM(S). i£Sp(Š) Pre a G< 0, 1) je funkcia h(a) = TrM_1(5) spojitá diferencovatelná, pričom h(0) = min{h(a) : a G< 0,1)}. Teda h(0) 5; h(a) pre a G< 0, 1) a musí byt lim -^-h(a) > 0 a^o da (ide o deriváciu sprava). Pomocou dôsledku 8.4 (položíme tam C = I) dostávame lim -^-TrM_1((l - a)5*A + ctó) = a^o da = - lim TrM-^íl-o;)^-!-^) -^-M((l-a)^ + a5) sa M-^íl-o;)^-!-^) = = - lim TrM_1((l - a)5*A + ač) a—^0 |[(l-a)M(ä*A) + aM(í)] M_1((l - a)<% + ač) = = —TVM-1 (íľ^) [—M(t>^) + M^M"1^) = (3.10) = TrM-1(^)[M(^) - M^lVr1^) ^ 0. Ak si vezmeme jednobodový návrh S so spektrom S p (ô) = {io}, tak z (3.10) pre z0 = 1,2,A0 je TrM-!(^)[M(^) - A^JM-1^) = = TrM-1^) - TrM-^AiJ^M-1^) = = TrM-1^) - Aíof/0M-2(^)fío ^ 0, čiže pre i = 1,2,Nq TVM_1(<ľ^) ^ Xlí-M^2(S*A)íl, teda (3.11) TrM-^D^mfliíA^M-2^)^ í = 1, 2,A0} = A{5*A). Pre lubovolný regulárny návrh S G Ares platí N0 N0 i=i i=i 14 (3.12) N0 i=l teda pre lubovolný regulárny návrh S G Areg je N0 i=i No i=i = max{\1ťl[M.-1{5)fí1 : i = 1, 2,N0} = A(ó), Z predpokladu TrM^1(S*A) = miniTrMr1 (S) : S e Areg} teda TrM'1^) S A(5). ■jA) = min{Tr'Mr1{5) : ó t t-±reg) TrMT1 (5*A) = míniTrMT1 {5) : 5 e Areg} ^ min{A{5) : 5 e Areg}. Ale z (3.11) máme Dostávame A(Ô*A) S TrM-1^) S mira{A(5) : 5 G Ares} ^ A(^), čím sme dokázali 1 2. Súčasne musí byt (3.13) A(<%) ^ TrM_1(^) = min{A(5) : S e Areg}, teda (3.14) A(^) = TrM-1(^) a dokázali sme 1. => 3. Teraz 2. 1. Nech S* G Ares minimalizuje -A(á) na množine Areg.Z predpokladu (3.15) TrM.-1^*) > míniTrMT1 {5) : ó e Areg} = TrM^1 (óA) vyplynie spor. Potom totiž musí existovat návrh ô, že pre funkciu s (a) = TrM_1((l - a) ô* +aô), a£<0,l) platí (3.16) lim -^-TrM_1((l - a)ô* + aô) < 0. 15 Keďže S* minimalizuje A(ô), je maxiX^iM.-1^*)}2^ : 2=1,2,N0} S SmaxiX^lM-1^)}2^ : i = 1, 2, N0} = = TrM-1^) = míniTrMT1 {5) : ó G Areg} a súčasne z platnosti (3.15) vyplýva (analogickou cestou ako odvodzovanie (3.10)), že lim -^-TrM_1((l - a)(T + aô) = a^>o da = TrM"1^*)^^*) - M^M"1^*) = N0 i=i (3.17) = TrM'1^*) - A(ô*) ^ TrM.^1 (ô*A) - A(ô*) ^ 0, čiže (3.17) je v spore s (3.16). Z toho vyplýva, že nemôže platit TrM_1(5*) > TrM"1^), čiže (3.18) TrM'1^*) S míniTrMT1 {5) : ó e Areg} = TrM^1 (ó*A). Samozrejme (3.19) míniTrMT1 {5) : ó G Areg} ^ TrM_1(5*) a dostávame, že TrM_1(5*) = mm{TrM_1(5) : 5 G Ares}. Dokázali sme 2. ^> 1. Konečne dokážeme sporom 3. => 1. 63 Nech 5** G Ares je taký návrh, ktorý neminimalizuje TrM 1(S), čiže TrM-1^**) > TrM-1^) = mm{TrM_1(5) : 5 G Ares}. Potom ale musí existovat návrh ô, že 5 = (1 — a)5** + at) má vlastnost, že (3.20) lim -^-TrM_1(5) < 0. Položme (3.21) lim -^-TrM^ÍÔ) = d < 0. 16 Podobnou cestou ako pri odvodzovaní (3.17) a s využitím (3.12) a (3.21) dostávame d N° d = lim —TVM_1(áh = TVM_1( 0, teda A( TrM"1^**). Dokázali sme 3. =^ 1. aj celú vetu. □ A—optimalita je špeciálnym prípadom L—optimality, ked funkcionál L(-) je definovaný pre daný pevný vektor h G lZk ako L(A) = h'Ah (pozri záver 2. kapitoly). Vetu o ekvivalencii pre L—optimálny návrh dokazujeme analogicky ako pre A—optimálitu (pozri [5], str. 68). Tu šiju len sformulujeme. Veta 3.3. Nasledujúce tvrdenia sú ekvivalentné. 1. Návrh S*L G Areg je L-optimálny, teda L[M_1(5£)] = mm{L[M_1(í)] : S G Areg } ■ 2. Návrh S*L G Areg minimalizuje l(S) = max{AiL[M_1(5)fjf/M_1(5)] : i = 1,2, ...,N0}, teda l(5*L) = min{l(ô) : S G Areg}. 3. l{5l) = L[M.-\5l)]. Veta o ekvivalencii pre reštringovaný A—optimálny plán vyplýva z vety 3.3, ked minimalizujeme L[M_1(5)] = Tr\M1'1(5)] (pozri (2.1). Teraz si ešte uvedieme vetu o ekvivalencii pre reštringovaný D—optimálny plán. Veta 3.4. Nasledujúce tvrdenia sú ekvivalentné. 1. Návrh 5*D G Areg je reštringovane D—optimálny, teda det\$A1,1{5*D )] = mmídeŕpVE1'1^)] : 5 G Areg}. 2. Návrh 5*D G Areg minimalizuje dr(ô) = max{XiťiM.^1{5)íi — \i{í\2))'M.^{5)í\2) : i = 1,2, N0}, teda dr(ó*Dr) = min{dr(ó) : ó G Areg}. Tu (3 = (/3^, /32)' a rozklad M(t>) a zodpovedá rozkladu vektora (3 na jednotlivé subvektory, 3. dr{5*D ) = k\ (dimenzia vektora f3\). Dôkaz, je zložitejší a vynechávame ho (pozri [2], str. 115). Skúsený uživatel môže na základe intuície alebo predchádzajúcich skúseností vypracovat taký návrh, že je blízko optimálneho (podlá zvoleného kritéria optimality). Napr. ak max {Xlí'iM.-1(5*)íl : i = 1,2,...,Nq} je len nepodstatne 17 väčší ako dimenzia k vektora (3, potom návrh S* je z hladiska praktického použitia D—optimálny a nemusíme ho už nijako vylepšovat. Podobne uvažujeme i pri ostatných vyššieuvedených kritériách. Používame pritom 3. tvrdenia viet o ekvivalencii. Ak predbežný návrh výrazne nespĺňa požiadavku 3. tvrdenia príslušnej vety o ekvivalencii, potom tento návrh iteračne vylepšíme postupom uvedeným v dalšej kapitole. 4. Iteračne určenie optimálneho návrhu Lema 4.1. Nech d(i,5) = f/M 1(ô)ti, pričom ô G Areg a i G {1, 2,N0}. Pre ľubovolný plán s G Areg, i G {1, 2,Nq} a a G (0, 1) platí detUí - a)M(S) + aXlílí'} = (1 - a)k i 1 + . a ,Xld(i, S) \ detM(S). I (1-a) J Dôkaz. Ak si zvolíme v leme 8.14 A = (1 — a)~M.(S), B = Va^Mi, C = — \fa\ií[, D = 1, dostávame z rovnosti detAdet(D - CA *B) = deťDdet(A - BD^C) tvrdenie lemy. □ Lema 4.2. Nech ô G Areg, 5 ^ 5*D. Potom max{det[(í - a)M(S) + a\&{-] : a G (0, í), i = 1, 2, N0} = x^d(z,,ó)\k ( k-1 Ý 1 detM{s) >detM{s)j k J yA^d^, 5) — 1 pričom i* je také číslo z množiny {1, 2,Nq}, že \i*d(i*, s) = max{\id(i, s) : i = 1,2, Nq}. Dôkaz. Z lemy 4.1 je vidiet, že det[{\ — a)~M.(S) + aXii^] je rastúca funkcia veličiny \id{i,5) (táto veličina nezávisí od a). Aby sme (pre daný plán S G Areg) dosiahli maximum det[{\ — a)~M.(5) + aXiiiil], musíme použit hodnotu Xi*d{i*,5) a maximalizovat (1 — a)k < 1 +---Xi*d(i*, S) > det~M.(S) vzhladom na a. Budeme maximalizovat Zndeí[(l — a)M(5) + aA^f^fj',] vzhladom na a. -^-lndet[(í - a)M(ó) + aX^i^i-,} = sa = -^-{fcžn(l-Q:) + žn{l+ ^ " . A,», ô) \ + žnrfeŕM(^) i = sa [ L (1 — aj J J /j 1 Xi,d(i*,ô) í — a l-al-ct+ a\i-d(i*, ô) Xl,d(i* ,S) — k oĺ = - [Ái,d(i*,S) - í]k 18 Pretože z predpokladu návrh S nie je D—optimálny, musí byt Xi»d{i*, S) — k > 0, teda a* > 0. Ľahko sa presvedčíme, že d2 -—lndet[(í - a)M(ô) + aXi.íi.í^,] da2 <0, čiže takto určený extrém je maximum. Teda max{det[(í - a)M(S) + aX&i-] : a G (0, í), i = 1,2, ...,N0} = = max{(l - a)k \ 1 + . " ,Xld(i, S) \ detM(S) : a G (0, í), i = 1, 2, N0} l (1-a) J = (1 - a*)k |l + ^yAi-d(i*, 5)| deŕM(5) = Xi.d(i*,S)\h f k-í Ý 1 detM{s) > detM{sy D k J \Xi,d(i*,ô)-l Veta 4.3. Nech pre postupnost plánov Sq, Si, ... platí Ss+i = (1 - a>l+1)6s + a*s+16i*+i, M(SS+1) = (1 - a*s+1)M(Ss) + a*s+1Xi:+iílUin:+i, (5i*+i je "jednobodový" plán so spektrom Sp(ôi* ) = {i*s+i})> pričom je určené z rovnice ^'s+1d(i*s+ii $ s) = max{Xld(i, Ss) : i = 1, 2, ...,N0} = = max{X£M-1(ós){l : i = 1, 2,N0} a * _ At;+irf(^+1A) - k as+1 ~ [Ai*+id(i*+1,5s) - í]k' Nech Si =/= S*D, i = 1,2,potom lim deŕM(5s) = detM(ó*D). Dôkaz, je zložitý a potrebuje hlbšie vniknut do teórie, pozri napr. [7],[2]. Preto ho tu vynecháme. V predchádzjúcich tvrdeniach opísaná optimálna volba čísel o;*, s = 1, 2, je dost zložitá. Jednoduchší postup pre volbu čísel as a tiež iteračný postup, ktorý sa osvedčil v praxi je nasledovný. Nech Sq je štartovací návrh. Prvé zlepšenie, ktoré vedie k návrhu Si je konvexná kombinácia plánu S0 a vhodne zvoleného "jednobodového" návrhu S^ so spektrom Sp(S$) = {it}, teda Si = (1 - a0)ó0 + a0óq , a0e(0,l). 19 Návrh 62 je opát konvexná kombinácia plánu t>i a vhodne zvoleného "jednobodového" návrhu 5i» so spektrom Sp(5i») = teda h = (1 -u1)51 +axSq, e (0,1). Takto postupujeme, až získame návrh Sopt, ktorý spĺňa zvolené kritérium optimality dostatočne presne. Volba štartovacieho návrhu S0, volba čísel a0, ai,... a volba postupnosti je daná nasledujúcimi pravidlami. Nech Sp(So) = «2, ífe} a hodnoty návrhu Sq v bodoch spektra nech sú ^o(i) = 7j i g Sp(S0). Je vhodné, aby štartovací plán mal v spektre práve k indexov, alebo len o málo väčší počet indexov ako k (k je dimenzia vektora parametrov (3). Volba štartovacieho plánu Sq musí byt taká aby jeho informačná matica M(50) bola regulárna. Regularita tejto matice je ekvivalentná nevychýlenej (nestrannej) odhadnutelnosti vektora parametrov (3 (pozri vetu 1.9). V dalšom budeme pokračovat tak, že počet bodov spektra štartovacieho plánu bude rovný k. lahko sa dajú upravit nižšie uvedené vztahy pre prípad inej volby (väčšieho počtu) bodov spektra štartovacieho plánu. Ak i\ g Sp(So), potom zrejme Sp(Si) = Sp(So). Hodnoty 5i(j) zvolíme podlá nasledujúceho pravidla o ak j = i\, k+V k+í 0, inak. Ak i\ £ Sp(50), potom Sp(S1) = Sp(S0) U {i^} a 1 ak j =/= i\ a súčasne j g Sp(S0), k+í SÁJ) = < 1 k+í 0, inak ak j = i\, ak j 7^ i\ a súčasne j g Sp(So), Úplne analogicky postupujeme v dalších iteráciach. Ak Sp(Sp) = {i1,i2,...,ikp}a, í*+1 g Sp(Sp), potom Sp(Sp+i) = Sp(ôp), Sp+1(í*+1) (k+p)5p(i;+l) + l k+p+1 a pre ostatné indexy j g Sp(Sp) je Sp+i(j) = (k + p)Sp(j) k + p+1 Ak í*+1 £ Sp(Sp), potom Sp(Sp+i) = Sp(Sp) U {í*p+1} 20 a pre j e Sp(Sp) je (k + p)Sp(j) k + p+1 a pre í*+1 je = k + p+1' Takto popísaný postup, ktorý určuje 5p+i = (1 — ap)Sp + apSi'+i, dáva pre čísla ap vztahy ap =---, p = 0,1,2,... . V literatúre [7],[2] sú i iné volby čísel k + p + 1 a0, ati,... . Určenie postupnosti i\,i*2, pri jednotlivých kritériách optimality: D—op tinmlita: Ai;+1fí;+1M-1(5fl)fi;+i = maxiXjťjM-^ôsVj : j = 1,2, ...,JV0}, S = 0,l,2,... . Reštringovaná D—optimalita pre určenie prvých k\ súradníc vektora (3: = maxlX^M-1^ - (ff ))'M^(5fl)ff>] : j = 1, 2, ...,N0}, S = 0,l,2,... . A—op timali ta: Ai;+1fí;+1[M-1(5fl)]2fi;+i = moaríAjfJpVI-1^)]2^ : j = 1,2, ...,JV0}, S = 0,l,2,... . L—optimalita: = mcuLiXjLlM^ios^jťjM-^Ss)] : j = l, 2,N0}, S = 0,l,2,... . Pri uvedenom sekvenčnom vylepšovaní štartovacieho návrhu S0 je potrebné na každom dalšom kroku znovu určit inverziu príslušnej informačnej matice plánu. Pri väčšom počte parametrov k a väčšom počte iterácií (niekedy rádovo 100 — 1000 iterácií) je iterovanie informačných matíc náročnou numerickou úlohou. Poznámka 4.4. Pri iteratívnom vylepšovaní štartovacieho návrhu Sq s výhodou používame vztah 21 (pozri Lemu 8.15). V našom prípade (ak v štartovacom návrhu bolo práve k experimentálnych bodov) je M(5fl+1)= W^fl+i(i) = XA^ôs^k++sS+í + ^+1^+1^:+! ^+1(^+1) TM(5fl) Preto k + s + 1 fc + a + Ai:+if/. M-i(5fl)fi:+i 5. Pravidlá pre zastavenie iterácií Pravidlá z predchádzjúcej kapitoly zaručujú konvergenciu iteračne vylepšovaných návrhov k návrhu optimálnemu. Konvergencia nemusí vždy postupovat tak rýchle, ako by sme si to priali, ale na druhej strane ani nepotrebujeme, aby proces dokonvergoval. Postačí nám vyhovujúce priblíženie k optimálnemu plánu. Toto priblíženie si stanovíme pravidlom zastavenia pomocou dostatočne malého kladného čísla s > 0. Iterácie zastavíme pri návrhu, ktorý označíme Sposi (posledný). D—optinmlita: Posledný návrh Sposi bude ten, pre ktorý platí max^X^M-^ópo^ : j = 1,2, ...,N0} <í+e. D cl Scl dokázat, že v tomto prípade platí k <\ maxi^Xj^M-^ôpos^fj : 3 = 1,2,N0}. A—optinmlita: Posledný návrh Sposi bude ten, pre ktorý platí mcu:{Ájq\M.-1(ôpoal)]2tj : j = 1,2, ...,N0} - TrM.-1^*) < s. L—optinmlita: Posledný návrh ôposi bude ten, pre ktorý platí maxiXjLiM-^Spo^i^M-1^)] : j = 1,2,N0} — L[M_1(5pos;)] < s. detM-^ópost) detM-1^*) 22 Reštringovaná D—optimalita (určujeme prvých k\ súradníc vektora (3): Posledný návrh ôposi bude ten, pre ktorý platí 1 max{~r~^i[fjM~1 (Sposi)íj - (í^YM^iôposi)^] : j = 1,2, ...,N0} <í + e. D cl Sel dokázat, že potom platí rietM1-1^) _ detM1'1(S*Dr) 1 fei < maxi^-X^M.-1^)^ - (ff tyM^(5poflí)ff >] : j = 1, 2,N0}. 6. S — optimalita Podia definície 2.7 návrh t)^ G Ares je S—optimálny ak 5e Areff,iV= 1,2, kde S je dopredu zadaná cielová kovariančná matica výsledného odhadu vektorového parametra (3. Hladáme 5^, a N*. Označme £, i = 1,2.....JV, G' = (í ÍJVr o)- Pri nejakom /V a návrhu t)^ má NNLO (3 kovariančnú maticu "^"M 1(S^), ktorej inverzia je N í M(<^) = G'N 0 5f,(2) V o G. 0 Ak existujú nezáporné čísla 5^(1), 5^(2),S^(N0), (pre ktoré = 1) že platí (6.1) G'jV 0 5|,(2) V 0 o 0 \ o 5*(/V0)/ G = S 23 tak sme našli aj N* = N. Podia lemy 8.18 sa dá systém (6.1) prepísat ako (G' G>ec 0 ... 0 \ ' 0 5* (2) ... 0 » V 0 0 S^(N0)J ktorý sa dá ešte zjednodušit tým, že sa vynechajú rovnice pre {E-1}^-, pre ktoré k(k — 1) i < j (dostaneme--- rovníc). Ak označíme y taký vektor, ktorý dostaneme keď vo vektore vecS^1 vynecháme súradnice {S-1}^j, pre i < j a A takú maticu, ktorú dostaneme ked v matici G' ® G' vynecháme práve tie isté riadky, ktoré zodpovedajú súradniciam vynechaným vo vektore iiecE-1, tak nájst S—optimálny návrh (pri nejakom N) je ekvivalentné hladaniu vektora x G 1ZN° íl M, kde N0 M = {x G KN° : x1^0,x2^ 0,...,xNo ^ 0,^%j = N}> i=i ktorý (vektor) minimalizuje veličinu ŕíľ(x) = x'A'Ax - 2y'Ax. Problém riešime ako úlohu kvadratického programovania (pozri napr. [3], str. 118). Riešenie x0 = N(6^,(1), 5^,(2),S^(Nq))' nám dáva optimálny návrh 5^ pri danom počte meraní N. Uvažujme teraz dva návrhy s rôznym počtom meraní. Pri prvom experimente nech je tento počet Ni, pri druhom N2- V obidvoch experimentoch nech sú matica A a vektor y rovnaké. Nech xi je riešenie našej úlohy pre prvý experiment a x2 je riešenie úlohy pre druhý experiment. Môže nastat situácia, že Ä"(xi) < k(tí2), teda s menším počtom meraní sa v prvom experimente lepšie priblížime k predpísanej matici S-1. Samozrejme dáme přednost návrhu prvého experimentu. Táto úvaha nás vedie k nasledovnej modifikácii definície S—optimálneho návrhu, s ktorou v praxi vystačíme. Vektor x0 = Nopt8^ nazveme približne S—optimálnym, ak minimalizuje veličinu Ä"(x) = x'A'Ax — 2y'Ax, jeho komponenty xo,i, xot2, xo,n0 su nezáporné a pre komponenty vektora á|. platí (*) = 1 ■ číslo Nopt je optimálny počet meraní, pričom Nopt = Nmax, kde Nmax je maximálny počet meraní, ktorý v danom experimente ešte pripúštame. Zrejme No ôh(i) = =ň02—, Nopt = ^2x0l. l^i=l x0i Ak použijeme označenie K (x) = x'A'Ax - 2y'Ax, B = i=l 1/ 1,1,..., 1 / \ Nmax No,No \ u — í UNo,l tak xo minimalizuje Ä"(x), pričom spĺňa podmienky Bxo íí b (toto označenie chápeme tak, že pre zvolený komponent vektora na lavej strane zodpovedajúci komponent vektora na pravej strane nie je menší). Problém riešime opát metódami kvadratického programovania. 24 7. Určenie optimálneho návrhu experimentu v niektorých špeciálnych prípadoch Niekedy máme za úlohu určit optimálny návrh experimentu, ked množina priamo observovatelných parametrov nie je konečná. Predpokladajme, že množina priamo observovatelných parametrov je {/j,x = í'x(3 : x g< a,b >}, kde íx = (1, x, x2,aľfe_1)'. Majme tiež dané váhy A^. ako A^. = , pričom aA(x) cr2(') Je polynóm, ktorý nemá korene v intervale < a, b > a jeho stupeň je menší alebo sa rovná 2(k — 1). V tomto prípade existuje D—optimálny návrh 5*D, že Sp(S*D) = {x(1),x(2),-;x^}. Pre návrh Ó*D platí S*D(x(i)) = l, * = 1,2,..., k. Dôkaz tvrdenia nájdete v [7]. Ako určit body x^\ x^2\ x^ spektra tohto D—optimálneho návrhu ? Dostaneme ich riešením maximalizačnej úlohy max 11 G'2 11 (X&-Xti)) x{1) e,...,x{k) e . - _, i>j za podmienky x^ < x^ < ... < x(k\ Ide o úlohu maximalizácie reálnej funkcie k reálnych premenných. Iný prípad, ked vieme napísat optimálny návrh je ak {/ix = í'x(3 : x g< a, b >}, kde íx = (í,x,x2, ...,xk^1)' a váhy A^ sú konštantné. Nech y g 1Z— < a, b >. Odhadujeme hodnotu íý(3 = Ej=1 VJ(predikcia, alebo extrapolácia). Optimálny návrh experimentu v tomto prípade je ten, ktorý vedie k minimalizácii disperzie odhadu čiže je to L—optimálny návrh, kde lineárny funkcionál L(-) je definovaný ako L(A) = iýAiy. Označme = COS (---7T I , i = 1,2, ...,k. \k- 1 J Hladaný L— optimálny návrh pre extrapoláciu s*l má spektrum pozostávajúce z bodov xi, = - + b+(b-a)z(^ L.,.....;. takých, že platí n {x(t)~v) Dôkaz pozrite v [7], kde sú dokázané aj iné podobné tvrdenia pre nájdenie optimálnych návrhov v niektorých špeciálnych situáciách. 25 8. Pomocné tvrdenia Nech B je regulárna n x n matica, ktorej prvky sú diferencovatelnými funkciami premennej t, čiže {B}i j = = bi j {ť), i,j = 1,2, ...,n, dB ~dt ■ , • , , dbuCt) je n x n matica, ktorej prvky su —^-, i, j = 1, 2,n dt ddetB . . ddetB 1 n je n x n matica, ktorej prvky su ——-, i, j = 1,2, ...,n, dB diagB /{B}i,i 0 ' 0 {B}2,2 0 0 Lema 8.1. Platí V 0 dt {B}ntn / B l9BB 1 Dôkaz. Prvky matice B označme 6^- , i,j = 1,2,...,n. Tiež sú diferencovatelnými funkciami premennej t, čiže b\- ^ = b\- ~^{ť),i,j = 1,2,..., n. Pre i, j = 1,2,..., n je {BB-1}ÍJ=^6ífe(í)6Í71)(í) = 5ÍJ fe=i (Sij je tzv. Kroneckerovo delta, čiže Síj = 0 pre i =/= j a t)jj = 1 pre z = j.) Preto 0 = §-{BB-%,=J2ÍWk(t)bk-J1)(t) k=l VÔMt)6<-i>m + v6 m^1)(ŕ) i, j = 1,2, ...,n, fe=i čo v maticovom zápise je dB ~dt k=l B 1 + B dB- = 0, Lema 8.2. iVec/i C je n x n matica konštánt. Platí ÔTrBC dB Dôkaz. ÔTrBC d dt dt = 757 ZZ M*)** = Z! Z! 96ý-(í) dB ,<9B 2—1 j — l i — l j — l Z predchádzajúcich dvoch liem priamo dostávame dt -c- = TrHtc = TrCl>t □ 26 Dôsledok 8.3. Platí dTrB ^ dB dt dt Dôsledok 8.4. Platí 9TrB-Ht)C_ _TrCB_1ÔB(t)B_1_ dt dt Lema 8.5. Platí ddetB dB (detB)(B ak je B nesymetrická (detB)(2B_1 — diagB^1), ak je B symetrická. Dôkaz. Determinant regulárnej n x n matice B sa dá písat ako detB = {B^.iSi.i + {B}ij2Bij2 + ... + {B}i>nSijn pre i G {1, 2, n}, pričom Bst je doplnok (n — l)-ho stupňa determinantu detB patriaci k prvku {B}st (pozri napr. [4], str. 270). Preto ddetB d »{»!.,, n{B\,., Dostávame ({Bj^iS^i + {B}ij2Sij2 {^}i,nBi,n)- ddetB = B ddetB dB ■82,1 S2,2 = (detB)(B-1)' lebo B ■82,1 Bn,i \ detB detB detB -Bl,2 -82,2 -B„,2 detB detB detB B2,n ^ 71,71 det B det B (pozri napr. [4], str. 320). Toto platí o nesymetrickej matici B. V prípade, že B je symetrická, platí {B}rjS = {B(6n, 612, ...,bin, b22,b: '23 , ■ ■ ■ , ^2n, ■ ■ ■ , , ^n —1 n—l,^n —1 íi,^íiíi)}r,s 6rs, ak r íí s, 6sr, ak r > s. 27 Pre symetrickú maticu teda B Preto /{B}M {B}1>2 ... {B}1>n\ {B}2jl {B}2,2 ... {B}2,„ V{B}nji {B}„j2 ... {B}nj„/ ddetB _ ™ ™ ddetB d{B}Kl dbu -2^2^d{B}k>l dbu /hi b12 bi2 622 bin \ b2n V bin b2n ■ ■ ■ bnn ) ddetB <9{B}M _ d{B}M dbu d d {B} i,i Pre i < j je dcfeíB [{B^.iSi.i + {B}íj2Síj2 + ... + {B}i>nSijn] .1 = Bij, í = 1,2,..., n. ddeíB <9{B}M ddeíB ^B}^- ddeíB 9{B}j-í d{B\ d{B}. [{Bj^iS^i + {B}ij2Sij2 + ... + {lS}i,nBi,n] -l" 9 9{B} {Bjj^Sj^] .1 Úplne rovnako pre z > j dostaneme ddetB díieíB dB ddetB ddetB ddetB dbu ddetB dbi2 ddetB dbin ddetB dbi2 db22 db2n ddetB ddetB ddetB dbin db2n dbnn (det B)(2B_1 - díagB-1). □ Lema 8.6. Pre symetrickú regulárnu n x n maticu B platí dlndetB(t) R_i<9B(ŕ) dt ~ r dt ' Dôkaz. Ak si uvedomíme, že B aj B-1 sú symetrické matice, teda pre í > j platí JdB 1 f (9B 1 —— > = < > a tvrdenie predchádzjúcej lemy, čiže {B-Hij = {B-1},, dostávame dlndetB(ť) ďt ddetB _ í 2{B-1}íJdeíB, ak i < j, dbtj \ {B_1}jjí(ieíB, ak i = j, 1 ddetB _ 1 ddetB dbtj _ T ß-i dB detB dt ~ detB ^ ^ dbtj dT ~ T ~dt' Í — 1 j—Í J □ 28 Lema 8.7. A, B nech sú symetrické m x m matice, A je pozitívne definitná. Potom existuje nesingulárna m x m matica U taká, že platí UAU' = I, UBU' = A, pričom A je diagonálna. Dôkaz. Označme w1,...,wm ortonormálně charakteristické vektory matice A a di,dm jej charakteristické čísla prislúchajúce týmto charakteristickým vektorom. Teda Awj = djWj, i = 1,2,to. Čiže existujú matice W = (wi:w2:... :wm) a D /di 0 ' 0 d2 ° \ 0 dm' alebo Označme Vo o . AW = WD, WW' = W'W = I, W'AW = D, WDW' = A, WD 1 W' = A 1 / Vď[ o ' o ^fd~2 Vo o 0 \ o nrz) o -k D 0 o 0 \ o ml C = WD5W'. Matica C je regulárna, symetrická, pričom C2 = WD i W' WD i W' = WDW' = A, teda WD5W' = CT1 a CT2 = A"1. Matica S = C^BC"1 je symetrická. Nech si,sm sú jej ortonormálně charakteristické vektory a Ai,Am charakteristické čísla prislúchajúce týmto vektorom, teda C 1BC 1s,=A,s,, i = 1,2,...,m Označme Zrejme C Si, i = 1, 2,to. A^Bz, = A^BC V = A^CC^BCTV = A^CA.s,- = C"1 A,-s,- = A Preto Zj, z = 1,2,...,to sú charakteristické vektory matice A ^aA^, i = 1,2,...,m sú im prislúchajúce charakteristické čísla, čiže A_1Bz,- = A, i = 1,2,..., TO 29 Pre maticu U' = (zi:z2:...:zm) platí UAU' = V z' / A(zi:z2:...:zm) = / ziAzi ziAz2 z2Azi z2Az2 iAzm \ VzmAzi z A Z2 lebo z^Azj = z^C'Czj = s^s,- = fc. Ďalej UBU' V z' / B(zi:z2:...:zm) = /A, 0 0 . A2 . . 0 Vo 0 . Am / = A, lebo z^Bzj = z^AA-^Zj = z^C2A"1Bzj = C ^BC 1 s j = Aj s j. = síC^CHC-^BC-^,- = □ Lema 8.8. Nech Mi, M2 sú mxm symetrické a pozitívne definitné matice, \i,...,\Tl sú charakteristické čísla matice M^M^1. Potom sú Ai, Xm reálne a kladné. Dôkaz. Ai,Am sú riešením rovnice det(M1'M.21 - AI) = 0, ktorá je ekvivalentná nasledujúcim rovniciam deŕ[(Mi - AM2)M2_1] = 0, deŕ(Mi - AM2) = 0, deŕ[M| (M^MiM^ - AI)Mf ] = 0, deí(M^MiM^ - AI) = 0. _ i _i Pretože M2 2MiM2 2 je pozitívne deŕinitná matica, všetky jej vlastné čísla, teda _ i _i riešenia rovnice deŕ(M2 2MiM2 2 — AI) = 0 sú reálne a kladné. □ Lema 8.9. Majme n x n maticu F a nech 71, 7„ sú korene jej charakteristickej rovnice, t.j. rovnice det(F — 7I) = 0. Potom detF = Yľj=i Ij a TrF = 5^™=i 7«- Dôkaz. Charakteristická rovnica det(F — 7I) = 0 sa dá písat ako (-l)n7n + M""1 + - + 6n-l7 + bn = 0, čiže (8.1) (-1)"(7" + «17"_1 + - + On-17 + On) = 0, pričom (priamo z definície determinantu) platí (8.2) bn = (-í)nan = detF, b\ = (-l)nai = (-l)r TrF. 30 Rovnicu (8.1) môžeme písat ako (8.3) (-l)n(7 - 7i)(7 - 72)-..(7 - 7n) = 0 a roznásobením (8.3) dostávame vztahy medzi koreňmi a koeŕicientami charakteristickej rovnice, teda -ai =7i + ...+7„(=TrF) 0-1 = 7172 + ••• + 7n-l7n -a3 = 717273 + ■■• + 7»-27n-l7n (8.4) (-l)"a„ = 7i72-7n(= Ď« = detF). Z (8.2) a (8.4) dostávame tvrdenie lemy. □ Lema 8.10. Nech Ai,Xn sú nezáporné čísla. Platí (8.5) (n^V ^í> n ■ i=l Dôkaz. Ak a ^ 0, 6^0, tak (8.6) ab^^(a2 + b2), lebo (a — Ď)2 ^ 0, teda postupne a2 - 2ab + b2 ^ 0 a2 + Ď2 ^ 2aĎ i (a2 + Ď2) ^ a&. Nech Ai ^ 0, A2 ^ 0, a položme v (8.6) a = V^T, 6 = Potom z (8.6) (8.7) VAiÄ^=a&^ -(Ai+A2). Nech ďalej Ai,a2,A3, A4 sú nezáporné. Podlá (8.7) platí \/AiA2 ú ^(Ai + A2) a \a3a4 ^ ^(A3 + A4), čiže (znovu využijúc (8.7)) \VA1A2A3A4 = \f x/xJ^x/X^ S y|(Äi + A2)i(A3 + A4) < 31 < ^(Ai + A2) + i(A3 + A4) Ai + A2 + a3 + a4 Matematickou indukciou dokážeme, že (8.5) platí pre n = 2k (k je prirodzené číslo). Pre k = 2 (8.5) platí. Nech platí (8.5) pre n = 2k, ukážeme, že platí aj pre n = 2k+1. Teda nech platí ' 2" \ 2" 2" IIM ^ÄEA 2k 1-1 = 1 I i=l Vezmime Ai, A2,A2fc, A2fc+1,A2fc+i=2 2t. Naozaj „fc+i \ 2fc+i n a, (2fc \ 2fc / 2fc J_ 2fc 2fc+i < \ 2~fcEAí^fc Ea2fc+i = \ i=l i=i < 2 2k 2" 2" EA* + ÄEA 2fc \ 2" ~lJ2^ 2k /-^i "2fc+j I 2fe+l i=i / í=l Teraz ukážeme, že ak (8.5) platí pre nejaké n ^ 2, tak platí aj pre n — 1 (spätná indukcia). Ai + A2 + ... + A„_i —. Všetky äí su nezáporné. Uvažujme Ai, A2,\n-i, K = n — 1 Teda ak nerovnost (8.5) platí pre (nejaké) n ^ 2, má tvar AiA2...A„_ Ai + A2 + ... + A„_i < Ai + A2 + ... + A„_i n — 1 Ai + A2 + ... + A„_i n — 1 n — 1 (Ai+A2 + ... + A„_i). Nerovnost (8.8) sa dá písat tiež ako čiže (ak aspoň jedno Ai > 0) ŕn— 1 \ n—1 i=i n —1 n n—l ^tE^ ^Ea - 1 4^ J ~ n - 1 ^ /n — l \ n — l < n—1 Umocnením na n — 1 dostávame ''n —1 \ n —1 i=i n-1 < 1 " ^ Tým sme lemu dokázali pre každé n = 1,2, i=l □ 32 Lema 8.11. Majme dva návrhy 61,62, pričom 61 G Areg. Pre a G< 0,1) patrí 6 = (1 — a)6i + 0:62 do Areg a M(6) = (1 - a)M(5i) + aM(52). Dôkaz. Najprv ukážeme, že S je návrh. Skutočne S je zobrazenie z {1, 2,Nq} —> < 0,1 >, pre ktoré platí N0 53<^) = Z (l-a)5i(i) + a52(i) = i=i ie{s,p(5i)us,p(Ä2)} = (l-a) 53 Sx(i) + a 53 S2(i) = (1 - a)+ a = 1. Informačné matice návrhov #i,52 sú M(5!)= ^ a M(52)= ^ ^./'Wi- pričom M(t>i) je pozitívne deŕinitná a M(52) je pozitívne semideŕinitná. Informačná matica návrhu S je M(~S)= J2 [(l-a)5!(i) + a52(i)]W? = = (l-a) ^ ^a^ + a ^ 52(i)Aifif/ = (l-a)M(51) + aM(52) a pre a G< 0,1) to je pozitívne deŕinitná matica (vyplýva priamo z definície pozitívnej deŕinitnosti matice). □ Lema 8.12. Majme návrhy 61,62, pričom Si G Areg, «£< 0,1), ô=(l — a)ôi + 0162 (G Areg podlá Lemy 8.11). Reálna funkcia g (a) = lndetM(6) = lndet[(í - a)M(51) + aM(52)} je spojitá a diferencovatelná na < 0,1). Dôkaz. Z definície determinantu vyplýva, že deťM.(ô) = deŕ[(l —a)M(#i)+o:M(52)] je polynóm k—teho stupňa premennej a, teda to je spojitá a diferencovatelná funkcia. Pretože pre a G< 0,1) je det[(í — a)'M.(Si) + o;M(52)] > 0 a In(-) je všade na (0, 00) spojitá a diferencovatelná, je g(-) (ako funkcia a) na < 0,1) spojitá a diferencovatelná. □ Uvedme jednu z charakteristických vlastností funkcie lndeťM.(ô) na množine {M(ô) : ô G Areg}. 33 Lema 8.13. Funkcia lndet(.) je na množině {M(t>) : S G Areg} konkávna. Dôkaz. Vezmime lubovolné ôi, S2 G Areg. Podlá lemy 8.7 existuje regulárna matica U, že UM(5i)U' = I a UM(ä2)U' = A (diagonálna). Pretože In(-) je na (0, oo) (rýdzo)konkávna funkcia, dostávame pre každé a G< 0, 1 > lndet[(í - a)M(5i) + aM(ó2)} = IndetV'1^ - a)l + aA]U'_1 = k = ln{det\J~2det[(\ - a)l + aA]} = lndet\J~2 + lnJJ[(í - a) í + a\,] = i=i k k = lndet\J~2 + ^ M(l - a)l + a^i] = lndet\J~2 + J^[(l - ajlní + alnXl] = i=i i=i = lndet\J~2 + alndetA = (1 — a)lndet\J~2 + a{lndet\J~2 + IndetA] = = (í-a)lndetM(ó1) + alndet\J-1\J'^1A = (í-a)lndetM(ó1)+alndetlJ-1AlJ'^1 = = (1 - a)lndeťM(51) + alndetM(S2). Nerovnost je ostrá, ak a G (0,1) a ak Aj ^ 1 aspoň pre jedno i G {1, 2,k}, t.j. ak M(5i) ^ M(S2). □ Lema 8.14. Nech A, D sú regulárne matice. Potom platí A B det C D detAdet(D - CA_1B) = detY)det{A - BD_1C). Dôkaz. Zrejme A B det C D = det I 0\ í A B\ íl A *B CA 1 I M C D H 0 I Tiež det det A B C D A 0 0 D — CA_1B detAdet(D - CA_1B). det 1 BD 1 0 I A B C D I 0 -D_1C I detlA BQD 1C ° ) = detDdet(A-BD-1C). □ Lema 8.15. Nech A je regulárna n x n matica, u, v G lZn. Platí: /» »-i A_1uv'A_1 (A + uv') 1 = A — 1 +v'A-!u ' Dôkaz. Lemu lahko dokážeme vynásobením matíc (A + uv') (A + uv') . □ 34 Definícia 8.16. Nech A — / «11 «21 aln \ »2n V aml ... í Kroneckerov súčin matíc A a B je , Br A (8) B / anB ai2B a2iB &22B \amlB am2B hi \brl alnB \ «2nB 62s Vlastnosti kroneckerovho súčinu matíc pozri napr. v [9]. Definícia 8.17. Ak napíšeme "pod seba" stĺpce matice K, povieme, že sme vykonali na matici operáciu vec. Teda vecKmtn = í;ec(ki:k2:...:k„) = Vk„/ Lenia 8.18. Pre matice príslušných rozmerov platí vecABC = (C (g) A>ecB, Tr AB = (vecB')'vecA. Dôkaz. Lemu dokážte ako cvičenie. Viac o optimálnom návrhu experimentu nájdete v monografii [7], príklady sú v [8]. V [7] sa nachádza aj obšírny zoznam dalšej literatúry k danej téme. References [1] Cramer, H., Mathematical Methods of Statistics, Princeton University Press, Princeton, N.j., 1946. [2] Fedorov, V., V., Theory of Optimal Experiment, Duxbury Advanced Series, Second Edition, Belmont CA, 1990. [3] Hamala, M., Nelineárne programovanie, ALFA, Bratislava, 1972. [4] Kořínek, V., Základy algebry, Nakladatelství ČSAV, Praha, 1953. [5] Kubáček, L., Kubáčková, L., Statistika a metrologie, Univerzita Palackého v Olomouci -vydavatelství, 2000. [6] Kubáčková, L., Kubáček, L., Kukuča, j., Pravděpodobnost a štatistika v geodézii a geofyzike, VEDA, Bratislava, 1982. [7] Pázman, A., Základy optimalizácie experimentu, VEDA, Bratislava, 1980. [8] Pázman, A., Mikulecká, j., Raffaj, j., Tokošová, M., Riešené situácie z navrhovania experimentov, ALFA, Bratislava, 1986. [9] Rao, C. R., Lineární metody statistické indukce a jejich aplikace, Academia, Praha, 1978.