Obsah 1 Neurčitost a citlivost 2 1.1 Analýza neurčitosti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.1.1 Pravděpodobnostní analýza neurčitosti . . . . . . . . . . . . . . . . . 4 1.1.2 Bayesova metoda . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.2 Zobrazovací metody a lokální analýza citlivosti . . . . . . . . . . . . . . . . . 6 1.2.1 Metoda konečných diferencí . . . . . . . . . . . . . . . . . . . . . . . 7 1.2.2 Přímá metoda . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.2.3 Lokální citlivosti a analýza neurčitosti . . . . . . . . . . . . . . . . . 8 1.3 Globální analýza citlivosti . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.3.1 Standardizované regresní koeficienty . . . . . . . . . . . . . . . . . . 9 1.3.2 Spearmanův koeficient pořadové korelace . . . . . . . . . . . . . . . . 10 1.3.3 Metody založené na rozptylu . . . . . . . . . . . . . . . . . . . . . . . 11 1.4 Ilustrační příklady . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 1.4.1 Příklad 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 1.4.2 Příklad 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 1.4.3 Příklad 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 1.4.4 Příklad 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 1.5 Příklady k procvičení . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 1 1 Neurčitost a citlivost 1.1 Analýza neurčitosti S výpočetní vědou je neodmyslitelně spjatá vlastnost zvaná neurčitost. Jakékoliv řešení reálného problému z praxe v sobě zahrnuje neurčitosti a není možné je odstranit1 . To je na jednu stranu trochu frustrující, na druhou stranu nám nezbývá, než se s tím smířit a „minimalizovat škody . Každé řešení problému, které vede na nějaké „důležité rozhodnutí musí sestávat z analýzy neurčitosti. V knize [5] jsme neurčitosti rozdělili do tří hlavních skupin na neurčitosti v matematickém popisu, neurčitosti v datech a neurčitosti v aplikaci modelu. V literatuře je možné najít různá dělení neurčitostí, vždy se však týkají tří zmíněných oblastí. Každý popis reality je nepřesný, podobně nepřesná jsou data (závislá jednak na přesnosti měření, ale také na jejich „reprezentativnosti – viz dále) a zejména v případě vědeckých výpočtů se neobejdeme bez použití ICT (a případných dalších nepřesných metod hledání řešení problému). Pro lepší pochopení a uvědomění si důležitosti analýzy neurčitosti uveďme tři jednoduché příklady. Příklad 1.1: Na vymezeném území žije populace nějakých živých organismů. Určete, jak se bude vyvíjet její velikost v následujícím období (dnech, měsících, ...), jestliže počáteční velikost populace v čase t = 0 je N(0) = 100 a vnitřní koeficient růstu populace je r = 0, 1. Řešení: V knize [5] jsme uvedli několik matematických modelů růstu populace živých organismů. Předpokládejme, že naše populace splňuje kvantovací a vzorkovací podmínku (viz [5]). Jedním z možných popisů jejího růstu je rovnice d dt N(t) = r · N(t), (1.1) kde N(t) je velikost populace v čase t. Řešení rovnice (1.1) má tvar: N(t) = N(0) · er·t . Pokud nemáme k dispozici žádnou další informaci o tom, na čem závisí růst populace (a jak), odpovědí na zadaný problém je, že se velikost populace bude v čase zvětšovat podle vztahu N(t) = 100 · e0,1·t . Pokud ovšem zjistíme, že na vymezeném území je omezená kapacita prostředí hodnoty K, je třeba pro predikci růstu populace použít jiný model, např. d dt N(t) = r · N(t) · 1 − N(t) K . (1.2) 1 Je však většinou možné je snížit. 2 V takovém případě by už však naše odpověď na zadaný problém byla, že velikost populace se „po určitém čase ustálí na hodnotě K. K této hodnotě se bude velikost populace přibližovat buď neustálým růstem, či neustálým poklesem v závislosti na hodnotě K v porovnání s hodnotou N(0) (viz [5]). Další změnu přístupu k modelování růstu populace může přinést doplňující informace například o výskytu nějakého predátora na vymezeném území či jiných vlivů ovlivňujících počet jedinců sledované populace. To vše zahrnujeme do neurčitosti v matematickém popisu modelu. Příklad 1.2: Uvažujme opět, že na vymezeném území žije populace nějakých živých organismů. V určitý den byla změřena její velikost, která činila N(0) = 100 jedinců. Nechť dále víme, že prostředí, v němž populace žije, uživí (dlouhodobě) maximálně 500 jedinců (tj. kapacita prostředí K = 500) a že pro vnitřní koeficient růstu populace r platí r ∈ [−0, 1; 0, 1]. Určete vývoj velikosti populace v následujícím obodbí (dnech, měsících). Řešení: Předpokládejme nyní, že růst sledované populace odpovídá rovnici (1.2). V tomto případě se nám však objevuje neurčitost v hodnotě vnitřního koeficientu růstu populace r, která způsobuje, že odpověď je opět velmi neurčitá. Velikost populace se bude buď postupně zvyšovat až k hodnotě K = 500 (pro r > 0) nebo naopak klesat až po úplné vymření populace (pro r < 0), či může zůstat konstantní v závislosti na hodnotě (pro r = 0). Příklad 1.3: Nalezněte pevný bod funkce f(x) = 11 · x − 2. Řešení: Nejprve si připomeňme, že pevným bodem x funkce f : M → M je prvek x ∈ M, pro nějž platí x = f(x). Nejjednodušší způsob nalezení správné odpovědi je zřejmě vyřešení rovnice x = 11 · x − 2, (1.3) kdy okamžitě získáme x = 0, 2. Zkusme však hledat řešení metodou prosté iterace (viz [4]) pomocí ICT s počáteční iterací rovnou správnému řešení, tj. x1 = 0, 22 . V systémech Maple a Matlab zkusme vypočítat několik prvních iterací. Nehledě na jejich počet bychom očekávali, že každá z iterací by měla mít hodnotu 0, 2. V systému Maple bychom pro výpočet a výpis prvních 15 iterací v prostředí Standard Worksheet mohli použít následující kód: x [ 1 ] : = 0 . 2 for i from 2 to 16 do x [ i ] := 11∗x [ i −1]−2 end do V Matlabu bychom pro totéž použili kód: x (1)=0.2; for i =2:16 x( i ) = 11∗x( i −1)−2; end x Zatímco v systému Maple dostaneme všechny iterace rovné správné hodnotě 0, 2, v Matlabu se dvanáctá iterace (při standardním zobrazení na 4 desetinná místa) již liší od správné 2 Počáteční iteraci značíme jako x1 kvůli Matlabu, v němž je pole indexováno od čísla 1. 3 hodnoty a tento rozdíl se s dalšími iteracemi stále zvyšuje. Například dvacátá iterace (pokud bychom si ji vypočetli) již má hodnotu řádově 104 ! Když si nastavíme příkazem format long zobrazení čísel s vyšší přesností, zjistíme, že již první počítané iterace nebyly přesně rovné hodnotě 0, 2. Na vině je zobrazení čísel v počítači a fakt, že pevný bod funkce f je odpuzujícím pevným bodem (viz [4]). Čísla jsou v počítači ukládána v binární soustavě, v níž má většina desetinných čísel nekonečný rozvoj. Již při uložení čísla tak ztrácíme přesnost. Relativně špatná podmíněnost úlohy a odpuzující vlastnost pevného bodu tuto nepřesnost při každé další iteraci zvětšují. V systému Maple se chyba neprojeví, neboť ten má implementováno ukládání čísel v desítkové soustavě, takže číslo 0, 2 je uloženo přesně. Předchozí tři příklady byly vybrány záměrně. V praxi nemusí (ale samozřejmě mohou) neurčitosti způsobovat takové rozdíly (neurčitosti, chyby) v řešení zadaného problému, navíc příklad 1.3 je trochu „umělý . Ovšem i v praxi se může stát, že řešení problému povede na nějakou iterativní metodu, či že bude provedeno několik výpočtů, které zpočátku zanedbatelnou chybu zmnohonásobí. Z uvedeného vyplývá, že analýza neurčitosti je velmi důležitá k ověření „správnosti (resp. „věrohodnosti ) nalezeného řešení. Neurčitosti můžeme analyzovat postupně po skupinách, jak jsme si je dříve rozdělili. Neurčitost v matematickém popisu souvisí s nedostatkem informací (resp. provedenými předpoklady) při sestavování modelu. Jejich analýza pak sestává z porovnání různých struktur matematického modelu s měřenými daty. Jedná se o nejsložitější část celé analýzy. Datová neurčitost zahrnuje nepřesnost měření parametrů vystupujících v modelu (způsobenou jak chybou měřícího přístroje, tak lidským faktorem) a nepřesnost vzniklou omezeným množstvím vzorků ze sledovaného systému. V závislosti na informacích o hodnotách parametrů reprezentujeme datovou neurčitost nejčastěji za pomoci intervalové aritmetiky, fuzzy teorie, či pravděpodobnostní analýzy [5]. Poslední typ neurčitosti je předmětem numerické analýzy. Zkoumá se, jaký vliv na hodnotu řešení má použití ICT (uložení čísel v počítači/softwarovém systému) a metod hledajících přibližné řešení. Více k této problematice je možné najít např. v [4]. 1.1.1 Pravděpodobnostní analýza neurčitosti Parametry matematického modelu obsahují neurčitost, kterou je sice většinou možné popsat intervalem, typicky má však parametr nějakou střední hodnotu, jež je v jistém smyslu „nejpravděpodobnější . Snahou tak je chápat parametr jako náhodnou veličinu s příslušným rozdělením pravděpodobnosti. Tzv. propagace neučitostí modelem (viz [5]) je pak zpravidla3 prováděna metodou Monte Carlo, při níž se vygeneruje „dostatečně mnoho realizací náhodných veličin reprezentujících parametry modelu, které poslouží k různým vyhodnocením modelu. Získáme tak pravděpodobnostní rozdělení hodnot řešení. Použití metody Monte Carlo „přináší další neurčitosti, neboť se jedná o přibližnou metodu. Se vzrůstajícím počtem generovaných realizací náhodných veličin a vyhodnocení modelu konverguje získané rozdělení hodnot ke skutečnému rozdělení hodnot řešení modelu. Přesnost získané aproximace pravděpodobnostního rozdělení testujeme porovnáním základních statistických parametrů (střední hodnota, rozptyl) u dvou různých skupin vygenerovaných řešení nebo použitím odhadu chyby Monte Carlo metody. Uvažujme model tvaru Y = f(X1, ..., Xm), 3 U (velmi) jednoduchých modelů, případně u výpočetně náročných, tomu tak být nemusí. 4 kde Y je hledané řešení, X1, ..., Xm (m ∈ N) jsou parametry modelu reprezentované náhodnými veličinami s pravděpodobnostními rozděleními p1(x1), ..., pm(xm). Označme X = (X1, ..., Xm), x = (x1, ..., xm) a p(x) sdruženou hustotu pravděpodobnosti náhodného vektoru X. Nechť dále M je počet vygenerovaných realizací náhodného vektoru X, X(k) jeho k-tá realizace. Pak pro odhad střední hodnoty funkce f(X) máme: E (f(X)) = 1 M · M k=1 f X(k) . (1.4) Tento odhad konverguje (podle silného zákona velkých čísel) ke skutečné střední hodnotě pro M → ∞, jestliže |f(x)| · p(x) dx < ∞. Rozptyl vM takto určených odhadů středních hodnot pro dané M zapíšeme jako: vM = V 1 M · M k=1 f X(k) = 1 M · V (f(X)) , kde jsme využili vlastnosti rozptylu, že pro nezávislé náhodné veličiny Z1, Z2 a konstanty a, b platí: V (a · Z1 + b · Z2) = a2 · V (Z1) + b2 · V (Z2), a současně předpokladu (f(x))2 · p(x) dx < ∞. Hodnotu √ vM nazýváme chybou Monte Carlo metody (viz [8]). Aplikací nevychýleného odhadu pro rozptyl V (f(X)) můžeme vM odhadnout vztahem vM = 1 M · (M − 1) · M k=1 f X(k) − E(f(X)) 2 , (1.5) případně vM = 1 M · (M − 1) · M k=1 f2 X(k) − M · E(f(X)) 2 . (1.6) 1.1.2 Bayesova metoda Bayesova metoda nám umožňuje na základě apriorní informace o hodnotě parametru (resp. řešení modelu) a experimentálních výsledků vyvodit posteriorní informaci o parametru (resp. řešení modelu). Mějme opět model tvaru Y = f(X1, ..., Xm) se stejným označením jako dříve. Dále označme Oj hodnotu j-tého (j = 1, ..., n, kden ∈ N) pozorování (experimentálního výsledku). Monte Carlo simulací jsme získali pravděpodobnostní rozdělení hodnot Y skládajícího se z celkem M vygenerovaných (a následně vyhodnocených) řešení, které označíme Y(k) pro k = 1, ..., M. Pravděpodobnost p Y(k) , že náhodně 5 zvolené řešení je rovno Y(k) pro konkrétní k, je 1 M . Dále definujme pravděpodobnost p O|Y(k) vyjadřující, že pozorování O (O = (O1, ..., On)) odpovídá predikci modelu Y(k), opět pro konkrétní k. Za předpokladu, že chyby pozorování jsou nezávislé a normálně rozdělené s nulovou střední hodnotou a rozptylem σ2 , zmíněnou pravděpodobnost spočítáme podle vztahu4 p O|Y(k) = n j=1 1 √ 2πσ2 · exp (Oj − Y(k))2 2 · σ . (1.7) Užitím Bayesova vzorce můžeme definovat tzv. posteriorní pravděpodobnost p Y(k)|O vygenerované hodnoty řešení Y(k): p Y(k)|O = p O|Y(k) · p Y(k) M i=1 p O|Y(i) · p Y(i) , (1.8) a získat tak nové rozložení hodnot řešení modelu Y více odpovídající uskutečněným pozorováním se střední hodnotou: E(Y ) = M i=1 p Y(i)|O · Y(i) a rozptylem V(Y ) = M i=1 p Y(i)|O · (Y(i) − E(Y ))2 . Zcela analogicky získáváme též nová rozdělení pravděpodobnosti pro parametry X1, ..., Xm. Jejich střední hodnoty a rozptyly počítáme dle stejných vztahů (pouhou záměnou veličiny Y za zvolený parametr). 1.2 Zobrazovací metody a lokální analýza citlivosti Analýzu neurčitosti obvykle doplňuje analýza citlivosti. Díky té se dozvíme, jak které parametry ovlivňují řešení modelu. Metody analýzy citlivosti můžeme rozdělit do tří skupin na metody zobrazovací, lokální a globální. Účelem zobrazovacích metod je identifikace nejvýznamnějších parametrů, které nejvíce ovlivňují řešení modelu. Zobrazovací metody bývají tzv. „one-at-a-time , což znamená, že v jednom kroku zkoumají vždy vliv jednoho parametru na řešení modelu, zatímco ostatní parametry nabývají svých středních hodnot (a jsou pro tento krok konstantní). Používají se především při analýze citlivosti modelů s velkým počtem paramerů (až stovky). Získaný výsledek je kvalitativní. Zjistíme, které parametry jsou citlivější (tj. ovlivňují více řešení modelu) než jiné, nedozvíme se však, o kolik. Výhodou je nižší výpočetní náročnost (oproti ostatním metodám analýzy citlivosti), nevýhodou pak zmíněná „nedostatečná informace o citlivosti parametrů. Více o zobrazovacích metodách analýzy citlivosti lze najít např. v [13]. Lokální analýza citlivosti byla již stručně představena v [5]. Deterministický matematický model bývá obvykle zadán systémem algebraických rovnic 0 = f(y(x), x), (1.9) 4 exp(z)=ez 6 kde y(x) = (y1(x), y2(x), ..., yn(x))T je n-rozměrný vektor hledaného řešení (závisle proměnných), x = (x1, x2, ..., xm)T je m-rozměrný vektor parametrů modelu a f = (f1, f2, ..., fk)T vektorová funkce; nebo systémem diferenciálních rovnic, jako je např. počáteční problém: dy(t, x) dt = f(y(t, x), x), y(0) = y0, (1.10) kde y(t, x) je hledané řešení, t čas, x je m-rozměrný vektor parametrů modelu, f je krozměrná vektorová funkce a y0 vektor počátečních podmínek. Jelikož je rovnice (1.9) speciálním případem rovnice (1.10) pro konstantní t, budeme nadále uvažovat pouze rovnici (1.10). Analogicky k [5] definujeme lokální citlivosti prvního řádu výrazem ∂yi(t, x) ∂xj pro i = 1, ..., n, j = 1, ..., m a lokální citlivosti druhého řádu výrazem ∂2 yi(t, x) ∂xj∂xl pro i = 1, ..., n, j, l = 1, ..., m. Nakonec definujme matici citlivostí S, která je tvořena lokálními citlivostmi prvního řádu, tedy S = ∂yi(t, x) ∂xj pro i = 1, ..., n, j = 1, ..., m. Pokud nejsme schopni určit popsané citlivosti analyticky pomocí uvedených vztahů, používáme některou z numerických metod. Uveďme pro ilustraci dvě nejjednodušší metody, další je možné nalézt v [13]. 1.2.1 Metoda konečných diferencí Nejjednodušší metoda analýzy lokální citlivosti spočívá v nahrazení derivace diferencí. ∂y(t, x) ∂xj ≈ y(t, x1, ..., xj + ∆xj, ..., xm) − y(t, x) ∆xj , j = 1, ..., m. (1.11) Přesnost této metody závisí na velikosti diference ∆xj. Čím větší hodnota, tím méně je výsledek přesný zejména u nelineárních modelů. Pokud je diference příliš malá, zvyšuje se vliv zaokrouhlovací chyby počítače. V praxi se často volí diference rovná 1 % velikosti příslušného parametru, nicméně nalezení optimální (resp. akceptovatelné) diference je zpravidla otázkou metody pokus-omyl [13]. 1.2.2 Přímá metoda Jestliže zderivujeme rovnici (1.10) podle parametru xj dostaneme (za použití pravidla o derivaci složené funkce) systém d dt ∂y ∂xj = ∂f ∂y · ∂y ∂xj + ∂f ∂xj , ∂ ∂xj y(0) = S0. (1.12) 7 Systém (1.12) můžeme dále přepsat do maticové formy d dt S = J · S + F, S(0) = S0, (1.13) kde S0 je nulová matice, J a F matice zadané vztahy J = ∂fl ∂yi , F = ∂fl ∂xj , pro l = 1, ..., k, i = 1, ..., n, j = 1, ..., m. Přímé metody jsou založeny na řešení rovnice (1.12). Pro numerické řešení je třeba znát matice J a F, a proto se rovnice (1.12) musí řešit společně s rovnicí (1.10). Pro asymptoticky stabilní stacionární bod nebo v případě řešení systému (1.9), můžeme rovnici (1.13) transformovat na tvar (viz [13]): S = −J−1 · F. (1.14) 1.2.3 Lokální citlivosti a analýza neurčitosti Závěrem části o lokálních metodách analýzy citlivosti si ještě ukažme jejich možné využití při analýze neurčitosti. Lokální citlivosti lze totiž použít k lineárnímu odhadu rozptylu hledaného řešení a vlivu jednotlivých parametrů na tento rozptyl. Rozptyl řešení yi(t, x) způsobený rozptylem parametru xj můžeme odhadnout vztahem σj 2 (yi) = ∂yi ∂xj · σ(xj) 2 . (1.15) Za předpokladu, že mezi parametry x1, ..., xm nedochází k interakcím, pak můžeme celkový rozptyl řešení yi vyjádřit jako součet právě vypočítaných rozptylů: σ2 (yi) = m j=1 σj 2 (yi). (1.16) Z rovnic (1.15) a (1.16) je následně možné určit příspěvek parametru xj do celkového rozptylu řešení yi způsobeného všemi parametry zlomkem Sij = σj 2 (yi) σ2(yi) , (1.17) případně v procentech S%ij = Sij · 100 %. (1.18) V následující části se budeme věnovat metodám globální analýzy citlivosti, které určí vliv parametrů na rozptyl v hledaném řešení (mnohem) přesněji. Přesto může být výhodné v některých případech použít představené lokální metody k zisku alespoň řádové přesnosti (a tedy prvotní představy). 8 1.3 Globální analýza citlivosti Analýza lokální citlivosti zkoumá lokální efekt změny parametrů (kolem jednoho bodu v m-rozměrném prostoru parametrů) na řešení modelu. Pokud však parametry nabývají širšího spektra hodnot, využívají se spíše metody globální analýzy citlivosti. Vraťme se pro jednoduchost k modelu popsaném rovnicí Y = f(X1, ..., Xm), (1.19) kde Y je hledané řešení a X1, ..., Xm(m ∈ N) parametry modelu reprezentované náhodnými veličinami s pravděpodobnostními rozděleními p1(x1), ..., pm(xm). Metodou Monte Carlo můžeme získat odhad pravděpodobnostního rozdělení pY (y) pro řešení Y a s ním i odhady střední hodnoty a rozptylu. V této části nás bude zajímat, jak se na zmíněném rozptylu podílejí parametry X = (X1, ..., Xm). Předpokládejme, že jsme provedli Monte Carlo simulaci, a máme tak k dispozici M realizací pro každou z náhodných veličin X1, ..., Xm a Y označené x11, ..., x1M , ..., xm1, ..., xmM a y1, ..., yM . Nejprve zmíníme dvě hojně používané metody v globální analýze citlivosti, které jsou však závislé na typu modelu. 1.3.1 Standardizované regresní koeficienty Jestliže je uvažovaný model (skoro) lineární, můžeme provést regresní analýzu a popsat jej rovnicí Y = b0 + m i=1 bi · Xi + ε, (1.20) kde bi jsou regresní koeficienty a ε chyba regresního modelu (rozdíl skutečné hodnoty původního modelu a regresního odhadu). Hodnoty regresních koeficientů získáme např. metodou nejmenších čtverců. Označme: y =    y1 ... yM    , x =    1 x11 · · · xm1 ... ... ... 1 x1M · · · xmM    , b =    b0 ... bm    , ε =    ε1 ... εM    . Minimalizací výrazu S(b) = (y − x · b)T · (y − x · b) dostaneme (viz ODKAZ) b = (xT · x)−1 · xT · y. (1.21) Označme dále SStot = M k=1 (yk − y)2 , SSreg = M k=1 ( ˆyk − y)2 , SSres = M k=1 ( ˆyk − yk)2 , kde ˆyk značí odhad hodnoty yk získaný z regresního modelu a y průměr hodnot yk (pro všechna k = 1, ..., M). Adekvátnost regresního modelu pak můžeme posoudit z hodnoty koeficientu determinace definovaného vzathem 9 R2 = SSreg SStot . Koeficient determinace nabývá hodnoty z intervalu [0, 1] a udává, jaký podíl neurčitosti (rozptylu) se podařilo regresí vysvětlit. Regresní model můžeme algebraicky přeformulovat na tvar Y − y ˆs = m i=1 bi · ˆsi ˆs · Xi − xi ˆsi , (1.22) kde Y = b0 + m i=1 bi · Xi, y = M k=1 yk M , ˆs = M k=1 (yk − y)2 M − 1 , xi = M k=1 xik M , ˆsi = M k=1 (xik − xi)2 M − 1 . Nyní koeficienty bi· ˆsi ˆs za předpokladu, že Xi (pro i = 1, ..., m) jsou nezávislé, vyjadřují vliv proměnné Xi na hledané řešení Y . Nazýváme je standardizovanými regresními koeficienty a značíme βi. Platí, že β0 = 0 a (viz např. [14]): m i=1 (βi)2 = R2 . Čím větší je hodnota R2 , tím lépe popisuje regresní model závislost řešení původního modelu na jeho parametrech. K posouzení významnosti regresních koeficientů slouží několik statistických testů (viz [13]), které však vychází z předpokladů, jež v tomto případě (globální analýzy citlivosti deterministických modelů) neplatí. Pro dané nastavení parametrů x1, ..., xm totiž dostaneme vždy stejné řešení y. 1.3.2 Spearmanův koeficient pořadové korelace V případě, že je model (1.19) (výrazně) nelineární a koeficient determinace lineárního regresního modelu malý, zavádí se často transformace pořadí. Po vygenerování M realizací každé z náhodných veličin X1, ..., Xm, Y se přiřadí pořadí realizacím veličiny Y (od 1 do M podle pořadí velikosti původní hodnoty ze všech realizací této veličiny) a provede se lineární regrese. Podobně je možné transformovat i náhodné veličiny X1, ..., Xm a provádět lineární regresi až poté. Další možný postup je přímý výpočet Spearmanových koeficientů pořadové korelace. Jednotlivým vygenerovaným realizacím (u příslušné náhodné veličiny) přiřadíme pořadí podle velikosti (od 1 do M). Koeficient pořadové korelace ri vyjadřující závislost mezi Y a Xi následně vypočítáme podle vztahu ri = 1 − 6 · M k=1 dk 2 M3 − M , (1.23) kde dk je rozdíl v přiřazaném pořadí k-té realizace Xi a k-té realizace Y . Koeficient korelace je číslo z intervalu [−1, 1]. Hodnota blízká jedné značí rostoucí závislost mezi Xi a Y , hodnota blízká minus jedné klesající (tj. s rostoucí hodnotou Xi klesá hodnota Y ). 10 Za předpokladu, že je model monotónní a že mezi vstupními proměnnými modelu nedochází k žádným interakcím, lze získat příspěvek Si parametru Xi do celkového rozptylu výstupu modelu Y jako (viz [9]): Si = ri 2 m j=1 rj 2 . (1.24) V případě použití pořadové transformace můžeme k posouzení adekvátnosti modelu a významnosti získáných regresních koeficientů použít opět koeficient determinace a statistické testy významnosti. Při užití přímého výpočtu Spearmanova koeficientu pořadové korelace se též nabízí některé statistické testy významnosti (viz např. [9], [19]). Nevýhodou transformace podle pořadí je zvýšení relativních vlivů parametrů prvního řádu, tj. dojde k potlačení vlivů způsobených interakcemi, což může vést k přehlédnutí vlivů parametrů, které působí především v interakcích. Tato nevýhoda je tím významnější, čím víc má model parametrů [13]. 1.3.3 Metody založené na rozptylu Předchozí dvě metody globální analýzy citlivosti jsou závislé na typu modelu. Následující dvě metody patří do skupiny metod založených na rozkladu rozptylu a na typu modelu nezávisí. Princip spočívá v rozložení rozptylu řešení V (Y ) na rozptyly tohoto řešení způsobené jak jednotlivými parametry, tak jejimi interakcemi. Vlivem prvního řádu parametru Xi na řešení Y rozumíme výraz Vi = VXi (EX−i (Y |Xi)), (1.25) kde zápis X−i vyjadřuje všechny parametry Xj, j = 1, ..., m kromě Xi. Upozorněme, že nemůžeme vliv prvního řádu definovat výrazem VX−i (Y |Xi), neboť jednak bychom určovali rozptyl řešení při zafixovaném parametru Xi na jedné hodnotě (ten však nabývá intervalu hodnot) a navíc by se nám mohlo stát, že VX−i (Y |Xi) > V (Y ). Naproti tomu pro dříve definovaný vliv parametru Xi výrazem (1.25) platí požadovaná nerovnost Vi ≤ V (Y ). Relativní hodnotu vlivu prvního řádu (v provnání s celkovým rozptylem řešení) pak nazveme indexem citlivosti prvního řádu a označíme Si: Si = VXi (EX−i (Y |Xi)) V (Y ) . (1.26) Za povšimnutí stojí, že v případě lineárního modelu dostáváme: Si = βi 2 , kde βi je standardizovaný regresní koeficient (viz část 1.3.1). Vliv druhého řádu parametrů Xi a Xj (pro i = j) na řešení Y definujeme výrazem Vij = VXi,j (EX−(i,j) (Y |Xi, Xj)) − Vi − Vj, (1.27) kde podobně X−(i,j) vyjadřuje všechny parametry Xk, k = 1, ..., m kromě Xi a Xj. Od celkového rozptylu způsobeného parametry Xi a Xj odečítáme na pravé straně (1.27) vlivy prvního řádu, abychom dostali pouze vliv způsobený interakcí mezi Xi a Xj. Následně můžeme definovat index citlivosti druhého řádu Sij: Sij = Vij V (Y ) . (1.28) 11 Tímto způsobem bychom mohli pokračovat až po vliv (index citlivosti) m-tého řádu. Celkově jsme tak popsali vlivy parametrů X1, ..., Xm na rozptyl řešení Y , přičemž platí (viz např. [15]): V (Y ) = m i=1 Vi + m−1 i=1 m j=2 i 1/2 + (1/ Pi )∗ a r c s i n ( s i n ( omega [ 1 ] ∗ s ) ) : G[ 2 ] : = s −> 1/2 + (1/ Pi )∗ a r c s i n ( s i n ( omega [ 2 ] ∗ s ) ) : with ( S t a t i s t i c s ) : S := RandomVariable ( Uniform(−Pi , Pi ) ) : Ss := Sample (S , M) : for i from 1 to M do A[ i ] := e v a l f (G[ 1 ] ( Ss [ i ] ) ) ; B[ i ] := e v a l f (G[ 2 ] ( Ss [ i ] ) ) ; end do : body := seq ( [A[ i ] , B[ i ] ] , i = 1 . . M) : p l o t ( [ body ] , s t y l e = point , axes = boxed , l a b e l s = [ t y p e s e t ( x [ 1 ] ) , t y p e s e t ( x [ 2 ] ) ] ) ; Z obrázku 1.1 vidíme, že v případě (a) je vyplněnost prostoru parametrů velmi „špatná , což je způsobeno nevhodnou volbou frekvencí ω1, ω2. Čím „více jsou ω1, ω2 nepoměrné, tím lepší dostáváme výsledek. (a) ω1 = 3, ω2 = 7 (b) ω1 = 11, ω2 = 21 (c) ω1 = π, ω2 = 21 Obrázek 1.1: Vykreslení grafů z příkladu 1.4. 7 U(0, 1) značí uniformní (tj. rovnoměrné) rozdělení na intervalu (0, 1). 16 Doplňme ještě zdrojový kód pro vykreslení případu (a) v Matlabu: M=5000; Ss = u n i f r n d (−pi , pi , [M 1 ] ) ; omega1=3; omega2=7; x1 = . 5 + asin ( sin ( omega1∗Ss ) ) / pi ; x2 = . 5 + asin ( sin ( omega2∗Ss ) ) / pi ; plot ( x1 , x2 , ’ . ’ ) xlabel ( ’ x [ 1 ] ’ ) ylabel ( ’ x [ 2 ] ’ ) Vyplněnost prostoru parametrů (resp. zachování sdruženého pravděpodobnostního rozdělení p(x) pro parametry X1, ..., Xm) je též ovlivněna volbou transformačních funkcí Gi. V literatuře jsou uváděné zejména transformace pro případ, že parametry X1, ..., Xm jsou (rovnoměrně) rozdělené na intervalu (0, 1): Gi(sin(ωi · s)) = x∗ i · (1 + ν∗ i · sin(ωi · s)), (1.48) Gi(sin(ωi · s)) = 1 2 + 1 π · arcsin(sin(ω1 · s)), (1.49) kde x∗ i je nominální (střední) hodnota parametru Xi a ν∗ i definuje koncové body intervalu neurčitosti. Pokud potřebujeme vytvořit transformaci Gi pro rozdělení pravděpodobnosti pi(x) parametru Xi, můžeme využít navrhnované rovnice, kterou by taková transformace měla splňovat (viz [13]): π · 1 − x2 i · pi(Gi) · dGi(xi) dxi = 1. (1.50) Všimněme si, že transformační funkce (1.49) splňuje rovnici (1.50). Příklad 1.5: Uvažujme opět model tvaru Y = f(X1, X2), pro nějž X1 ∼ U(0, 1), X2 ∼ U(0, 1). Nechť ω1 = 11, ω2 = 21. Vykreslete graf získaných bodů [x1, x2] po vygenerovaní M = 5000 realizací veličiny S ∼ U(−π, π) a aplikaci funkcí G1, G2 pro (a) G1(sin(ω1 · s)) = 1 2 · (1 + sin(ω1 · s)), G2(sin(ω2 · s)) = 1 2 · (1 + sin(ω2 · s)), (b) G1 = 1 2 + 1 π · arcsin(sin(ω1 · s)), G2 = 1 2 + 1 π · arcsin(sin(ω2 · s)). Řešení: Zdrojový kód k vykreslení grafů je analogický příkladu 1.4. (a) ω1 = 11, ω2 = 21 (b) ω1 = π, ω2 = 21 Obrázek 1.2: Vykreslení grafů z příkladu 1.5. 17 Počet vygenerovaných hodnot parametru S by měl splňovat Shanonův8 vzorkovací teorém (viz [18]): „Přesná rekonstrukce spojitého, frekvenčně omezeného, signálu z jeho vzorků je možná tehdy, pokud byl vzorkován frekvencí alespoň dvakrát vyšší, než je maximální frekvence rekonstruovaného signálu. V našem případě proto nastavujeme hodnotu M minimálně na (viz [13]): M = 2 · N · max i=1,...,m {ωi} + 1. Ukázali jsme, jak s pomocí rozkladu modelové funkce f do fourierovy řady odhadnout celkový rozptyl funkce f v závislosti na parametrech X1, ..., Xm a rozptyly způsobené pouze jednotlivými parametry (tj. vlivy prvního řádu). Ještě bychom potřebovali odhadnout celkový vliv parametru Xi včetně jeho interakcí s ostatními parametry. K tomu využijeme stejné myšlenky jako dříve. Parametru Xi přiřadíme frekvenci ωi „velmi odlišnou od všech ostatních frekvencí příslušných zbylým parametrům, které budou navíc „téměř shodné . Například bychom tedy volili: ωi = 20 a ωj = 1 pro j = i. Jak jsme však viděli dříve, volit stejná (resp. „málo nepoměrná ) ωi pro i = 1, ..., m není vhodné kvůli zachování pravděpodobnostního rozdělení parametrů X1, ..., Xm. Zpravidla proto volíme vyšší hodnotu pro ωi a ostatní ωj (j = i) rozdělíme rovnoměrně v intervalu [1, ωi 2·N ], viz [13]. Vliv V−i všech ostatních parametrů pak odhadneme jako: V−i = 2 · N p=1 m j=1 j=i a2 p·ωj + b2 p·ωj . (1.51) Výhodou FAST metody je, že pro výpočet koeficientů citlivosti prvního řádu stačí jediný soubor vyhodnocení modelu (díky nastavení ω1, ..., ωm). Tato výhoda se smazává ve chvíli, kdy chceme počítat i celkové indexy citlivosti. Pro ně je totiž třeba vždy přenastavit frekvence ω1, ..., ωm (pro každý parametr Xi). Pokud však vytváříme frekvence ω1, ..., ωm podle výše popsaného postupu, lze pro jedno jejich nastavení počítat současně Si a STi . 1.4 Ilustrační příklady V této části si ukážeme čtyři jednoduché příklady, na nichž demonstrujeme dříve uváděné metody analýzy neurčitosti a citlivosti. 1.4.1 Příklad 1 Uvažujme model tvaru: Y = X1 + X2 + X3, (1.52) kde X1 ∼ U(1 2 , 3 2 ), X2 ∼ U(3 2 , 9 2 ), X3 ∼ U(9 2 , 27 2 ). Analýza neurčitosti V tomto jednoduchém příkladu bychom mohli střední hodnotu, rozptyl i pravděpodobnostní rozdělení veličiny Y vypočítat přímo na základě znalostí ze základního kurzu teorie 8 Též známý jako Nyquistův či Kotělnikovův. 18 pravděpodobnosti. Pokud využijeme nějaký systém počítačové algebry, tj. např. Maple, získáme požadované následovně: with ( S t a t i s t i c s ) : X[ 1 ] := RandomVariable ( Uniform (1/2 , 3 / 2 ) ) : X[ 2 ] := RandomVariable ( Uniform (3/2 , 9 / 2 ) ) : X[ 3 ] := RandomVariable ( Uniform (9/2 , 2 7 / 2 ) ) : Y := X[ 1 ] + X[ 2 ] + X [ 3 ] : Mean(Y) ; # s t r e d n i hodnota Y Variance (Y) ; # r o z p t y l Y PDF(Y, x ) ; # pravdepodobnostni d i s t r i b u c n i funkce Y Dostaneme: E(Y ) = 13, V (Y ) = 91 12 . Lokální analýza citlivosti Přímo z definice získáváme: S1 = ∂Y ∂X1 = 1, S2 = ∂Y ∂X2 = 1, S3 = ∂Y ∂X3 = 1. Při použití nepřímé metody (viz 1.2.1) je nutné zvolit bod [x1, x2, x3], v němž lokální citlivosti numericky aproximujeme. Nechť [x1, x2, x3] = [E(X1), E(X2), E(X3)] = [1, 3, 9]. Velikost diference zvolíme na obvyklé 1 % velikosti parametru (tj. ∆x1 = 0.01, ∆x2 = 0.03, ∆x3 = 0.09) a pro výpočet citlivosti prvního parametru použijeme vztah S1 = E(|x1 ± ∆x1 + x2 + x3 − x1 − x2 − x3|) x1 + x2 + x3 · 100, (1.53) který jednak vztahuje citlivost na jednotku velikosti parametru (tj. počítá relativní citlivost) a současně uvažuje jak vliv přičtení konečné diference, tak jejího odečtení9 . Hodnoty S2 a S3 se získají analogicky. Celkem dostaneme: S1 = E(0, 01; 0, 01) 13 · 100 = 1 13 , S2 = E(0, 03; 0, 03) 13 · 100 = 3 13 , S3 = ... = 9 13 . Vidíme, že získáváme jiné výsledky. Metodou konečných diferencí jsme de facto určili ne citlivost, ale elasticitu. Ta bývá definována vztahem: ∂Y ∂Xi · Xi Y resp. ∂ ln(Y ) ∂ ln(Xi) . (1.54) Stejných hodnot jak z definice citlivosti bychom dosáhli, kdybychom uvažovali diferenci rovnou jedné pro všechny tři parametry (tj. ∆x1 = ∆x2 = ∆x3 = 1) a nevztahovali výsledek k hodnotě řešení. Je ke zvážení, který přístup (citlivost/elasticita) je vhodnější k posouzení (lokálního) vlivu parametru na řešení. 9 a tyto dva vlivy zprůměruje (aplikací střední hodnoty) 19 Globální analýza citlivosti Model (1.52) je lineární, takže R2 = 1. Standardizované regresní koeficienty můžeme určit z rovnice Y − E(Y ) V (Y ) = β1 · X1 − E(X1) V (X1) + β2 · X2 − E(X2) V (X2) + β3 · X3 − E(X3) V (X3) . Dostaneme tak β1 = 1 √ 91 , β2 = 3 √ 91 , β3 = 9 √ 91 . Použitím rozkladu rozptylu dostáváme: V1 = VX1 (E(Y |X1)) = V (X1) = 1 12 , V2 = VX2 (E(Y |X2)) = V (X2) = 9 12 , V3 = VX3 (E(Y |X3)) = V (X3) = 81 12 . a následně S1 = 1 91 , S2 = 9 91 , S3 = 81 91 . 1.4.2 Příklad 2 Uvažujme model tvaru: Y = X1 + X4 2 , (1.55) kde X1 ∼ U(0, 5), X2 ∼ U(0, 5). Analýza neurčitosti Analýzu neurčitosti provedeme stejným způsobem jako v předešlém příkladu. Tentokrát získáme: E(Y ) = 255 2 , V (Y ) = 1000075 36 . Lokální analýza citlivosti Přímo z definice získáváme: S1 = ∂Y ∂X1 = 1, S2 = ∂Y ∂X2 = 4 · X3 2 . Tentokrát vychází index citlivosti S2 závislý na parametru X2, je tedy podstatné, v jakém bodě lokální citlivost vyhodnotíme. 20 Globální analýza citlivosti Tentokrát již náš model lineární není (v proměnný X1, X2), v případě počítání regresních koeficientů jej musíme vytvořit. Vygenerujme proto M = 105 realizací veličiny X1 a stejný počet realizací veličiny X2. Následně z nich vytvořme M řešení modelu (1.55). V systému Maple získáme lineární model příkazem LinearFit, v Matlabu můžeme použít příkaz regress. Maple M := 10ˆ5: with ( S t a t i s t i c s ) : X[ 1 ] := RandomVariable ( Uniform ( 0 , 5 ) ) : X[ 2 ] := RandomVariable ( Uniform ( 0 , 5 ) ) : X1s := Sample (X[ 1 ] , M) : # generovani M r e a l i z a c i v e l i c i n y X[ 1 ] X2s := Sample (X[ 2 ] , M) : # generovani M r e a l i z a c i v e l i c i n y X[ 2 ] Ys := Vector (M) : for i from 1 to M do Ys [ i ] := X1s [ i ]+X2s [ i ] ˆ 4 end do : # uprava formatu hodnot X1s a X2s pro potreby prikazu L i n e a r F i t s l o u p e c 1 := LinearAlgebra [ Transpose ] ( X1s ) : s l o u p e c 2 := LinearAlgebra [ Transpose ] ( X2s ) : Matice := Matrix (M, 2 , [ sloupec1 , s l o u p e c 2 ] ) : # l i n e a r n i r e g r e s e f i t := L i n e a r F i t ( [ 1 , x [ 1 ] , x [ 2 ] ] , Matice , Ys , [ x [ 1 ] , x [ 2 ] ] , output = solutionmodule ) : funkce := f i t :− R e s u l t s ( [ l e a s t s q u a r e s f u n c t i o n ] ) ; # l i n e a r n i model Matlab M=100000; X1s = u n i f r n d ( 0 , 5 , [M 1 ] ) ; % g e n e r o v a n i M r e a l i z a c i v e l i c i n y X[ 1 ] X2s = u n i f r n d ( 0 , 5 , [M 1 ] ) ; % g e n e r o v a n i M r e a l i z a c i v e l i c i n y X[ 2 ] Ys = X1s + X2s . ˆ 4 ; % model b=r e g r e s s (Ys , [ ones (M, 1 ) , X1s , X2s ] ) % r e g r e s n i k o e f i c i e n t y Získáme lineární model tvaru: Y . = 0, 59 · X1 + 99, 96 · X2 − 123, 97. (1.56) Dalšími výpočty dostáváme: R2 . = 0, 75, β1 . = 0, 01, β2 . = 0, 87. Jestliže přiřadíme všem realizacím (u příslušných veličin) pořadí10 (podle velikosti) a provedeme lineární regresi znovu, dostaneme koeficient determinace R2 . = 0, 98 a model tvaru: Y . = 0, 08 · X1 + 0, 99 · X2 − 3390, 25. (1.57) Další možností je výpočet Spearmanových koeficientů podle vztahu (1.23), po němž ob- držíme: r1 . = 0, 07, r2 . = 0, 99, S1 . = 0, 01, S2 . = 0, 99. Pro posouzení kvalitativního vlivu proměnných X1 a X2 na řešení Y můžeme zobrazit grafy závisloti Y na X1 a Y na X2, které nám v tomto případě velmi pomohou. Současně je možné zobrazit generované realizace veličin Y , X1 a X2 s lineární funkcí (1.56). K výpočtu Spearmanových koeficientů a vykreslení grafů z obrázku 1.3 můžeme použít následující kód: 10 V systému Maple k tomu lze využít příkaz Rank, v Matlabu příkaz sort. 21 (a) Generované body a lineární model (1.56) (b) Závislost Y na X1 (c) Závislost Y na X2 Obrázek 1.3: Vykreslení grafů z ilustračního příkladu 2. Maple R1 := Rank( X1s ) : # poradi prvku X1s podle v e l i k o s t i R2 := Rank( X2s ) : # poradi prvku X2s podle v e l i k o s t i RY := Rank( Ys ) : # poradi prvku Ys podle v e l i k o s t i # Spearmanovy k o e f i c i e n t y r [ 1 ] := e v a l f (1−6∗add ( (RY[ i ]−R1 [ i ] ) ˆ 2 , i = 1 . . M) / (Mˆ3−M) ) ; r [ 2 ] := e v a l f (1−6∗add ( (RY[ i ]−R2 [ i ] ) ˆ 2 , i = 1 . . M) / (Mˆ3−M) ) ; # Vygenerovane body a l i n e a r n i model body := seq ( [ X1s [ i ] , X2s [ i ] , Ys [ i ] ] , i = 1 . . M) : p l o t 1 := plot3d ( funkce , x [ 1 ] = 0 . . 5 , x [ 2 ] = 0 . . 5 , axes = frame , l a b e l s = [ t y p e s e t ( x [ 1 ] ) , t y p e s e t ( x [ 2 ] ) , t y p e s e t ( y ) ] ) : p l o t 2 := p l o t s [ p o i n t p l o t 3 d ] ( [ body ] , axes = frame ) : p l o t s [ d i s p l a y ] ( { p l o t 1 , p l o t 2 } ) ; # Graf z a v i s l o s t i Y na X[ 1 ] body2 := seq ( [ X1s [ i ] , Ys [ i ] ] , i = 1 . . M) : p l o t s [ p o i n t p l o t ] ( [ body2 ] , axes = frame , l a b e l s = [ t y p e s e t ( x [ 1 ] ) , t y p e s e t ( y ) ] ) ; # Graf z a v i s l o s t i Y na X[ 2 ] body3 := seq ( [ X2s [ i ] , Ys [ i ] ] , i = 1 . . 1 0 0 0 0 ) : p l o t s [ p o i n t p l o t ] ( [ body3 ] , axes = frame , l a b e l s = [ t y p e s e t ( x [ 2 ] ) , t y p e s e t ( y ) ] ) ; Matlab % Spearmanovy k o e f i c i e n t y r 1 = c o r r ( X1s , Ys , ’ type ’ , ’ Spearman ’ ) r 2 = c o r r ( X2s , Ys , ’ type ’ , ’ Spearman ’ ) % Vygenerovane body a l i n e a r n i model [ x , y ] = meshgrid ( 0 : 0 . 1 : 5 , 0 : 0 . 1 : 5 ) ; surf ( x , y , b ( 1 ) + b (2)∗ x + b (3)∗ y ) % l i n . model t i t l e ( ’ Graf vygenerovanych bodu a l i n e a r n i model ’ ) xlabel ( ’X[ 1 ] ’ ) ylabel ( ’X[ 2 ] ’ ) zlabel ( ’Y ’ ) set ( get ( gca , ’ ZLabel ’ ) , ’ Rotation ’ , 0 . 0 ) hold on plot3 ( X1s , X2s , Ys , ’ . ’ ) % vygenerovane body figure % novy o b r a z e k subplot ( 1 , 2 , 1 ) ; plot ( X1s , Ys , ’ . ’ ) ; % Graf z a v i s l o s t i Y na X[ 1 ] t i t l e ( ’ Graf z a v i s l o s t i Y na X[ 1 ] ’ ) xlabel ( ’X[ 1 ] ’ ) ylabel ( ’Y ’ ) set ( get ( gca , ’ YLabel ’ ) , ’ Rotation ’ , 0 . 0 ) subplot ( 1 , 2 , 2 ) ; plot ( X2s , Ys , ’ . ’ ) % Graf z a v i s l o s t i Y na X[ 2 ] t i t l e ( ’ Graf z a v i s l o s t i Y na X[ 2 ] ’ ) xlabel ( ’X[ 2 ] ’ ) ylabel ( ’Y ’ ) set ( get ( gca , ’ YLabel ’ ) , ’ Rotation ’ , 0 . 0 ) Použitím rozkladu rozptylu dostáváme: 22 V1 = VX1 (E(Y |X1)) = V (X1) = 25 12 , V2 = VX2 (E(Y |X2)) = V (X4 2 ) = 250000 9 . a následně S1 = V1 V (Y ) = 3 40003 . = 0, 0001; S2 = V2 V (Y ) = 40000 40003 . = 0, 9999. 1.4.3 Příklad 3 Uvažujme model tvaru: Y = sin(X1) + 7 · sin2 (X2) + 1 10 · X3 · sin(X1), (1.58) kde X1 ∼ U(−π, π), X2 ∼ U(−π, π), X3 ∼ U(−π, π). Analýza neurčitosti E(Y ) = 7 2 , V (Y ) = 1 1800 · π8 + 1 50 · π4 + 53 8 . Pravděpodobnostní rozdělení veličiny Y můžeme odhadnout vygenerováním realizací pro X1, X2, X3, následným vyhodnocením a zobrazením histogramu (v Maple příkaz Histogram, v Matlabu příkaz hist). Pro 105 vyhodnocení získáme v Maple histogram z obrázku 1.4 Obrázek 1.4: Histogram z ilustračního příkladu 3. Lokální analýza citlivosti Přímo z definice získáváme: S1 = ∂Y ∂X1 = cos(X1) + 1 10 · X4 3 · cos(X1), S2 = ∂Y ∂X2 = 14 · sin(X2) · cos(X2), S3 = ∂Y ∂X3 = 2 5 · X3 3 · sin(X1). 23 Globální analýza citlivosti Pokusíme-li se znovu vytvořit lineární model, dostaneme R2 . = 0, 19. Přiřazení pořadí tentokrát nepomůže, R2 „zůstane na stejné hodnotě. Použitím rozkladu rozptylu dostáváme: V1 = VX1 (E(Y |X1)) = V sin(X1) + π4 50 · sin(X1) = 1 5000 · π8 + 1 50 · π4 + 1 2 , V2 = VX2 (E(Y |X2)) = V −7 · cos2 (X2) = 49 8 , V3 = VX3 (E(Y |X3)) = V 7 2 = 0. a následně S1 = V1 V (Y ) . = 0, 3139; S2 = V2 V (Y ) . = 0, 4424; S3 = V3 V (Y ) = 0. Přestože jsme právě vypočítali globální indexy citlivosti 1. řádu, vypočítejme je nyní pomocí Sobol’ovy metody. Postupně dostáváme: E(Y ) . = 3, 4958; V (Y ) . = 13, 9476; (V (Y ) . = 13, 8446) ; V1 . = 4, 3680; V2 . = 6, 1619; V3 . = −0, 0561; S1 . = 0, 3131; S2 . = 0, 4418; S3 . = −0, 0040. Jestliže aplikujeme korekční člen (1.38), jehož hodnota je přibližně −0, 0571, získáme: V1 . = 4, 4250; V2 . = 6, 2190; V3 . = 0, 0010; S1 . = 0, 3173; S2 . = 0, 4459; S3 . = 0, 0001; Předchozí výsledky lze získat následujícími zdrojovými kódy: 24 Maple M := 10ˆ5: with ( S t a t i s t i c s ) : X[ 1 ] := RandomVariable ( Uniform(−Pi , Pi ) ) : X[ 2 ] := RandomVariable ( Uniform(−Pi , Pi ) ) : X[ 3 ] := RandomVariable ( Uniform(−Pi , Pi ) ) : X1s := Sample (X[ 1 ] , M) : # M r e a l i z a c i v e l i c i n y X[ 1 ] X2s := Sample (X[ 2 ] , M) : # M r e a l i z a c i v e l i c i n y X[ 2 ] X3s := Sample (X[ 3 ] , M) : # M r e a l i z a c i v e l i c i n y X[ 3 ] Ys := Vector (M) : for i to M do Ys [ i ] := s i n ( X1s [ i ])+7∗ s i n ( X2s [ i ])ˆ2+(1/10)∗ X3s [ i ]ˆ4∗ s i n ( X1s [ i ] ) end do : f [ 0 ] := Mean( Ys ) : # odhad s t r e d n i hodnoty V[ 0 ] := Variance ( Ys ) : # odhad r o z p t y l u X1s2 := Sample (X[ 1 ] , M) : # ” druha skupina ” M r e a l i z a c i v e l i c i n y X[ 1 ] X2s2 := Sample (X[ 2 ] , M) : # ” druha skupina ” M r e a l i z a c i v e l i c i n y X[ 2 ] X3s2 := Sample (X[ 3 ] , M) : # ” druha skupina ” M r e a l i z a c i v e l i c i n y X[ 3 ] Ys2 := Vector (M) : # c i t l i v o s t parametru X[ 1 ] for i to M do Ys2 [ i ] := s i n ( X1s [ i ])+7∗ s i n ( X2s2 [ i ])ˆ2+(1/10)∗ X3s2 [ i ]ˆ4∗ s i n ( X1s [ i ] ) end do : V[ 1 ] := add ( Ys [ i ] ∗ Ys2 [ i ] , i = 1 . . M)/M−f [ 0 ] ˆ 2 ; S [ 1 ] = V[ 1 ] /V [ 0 ] ; # c i t l i v o s t parametru X[ 2 ] for i to M do Ys3 [ i ] := s i n ( X1s2 [ i ])+7∗ s i n ( X2s [ i ])ˆ2+(1/10)∗ X3s2 [ i ]ˆ4∗ s i n ( X1s2 [ i ] ) end do : V[ 2 ] := add ( Ys [ i ] ∗ Ys3 [ i ] , i = 1 . . M)/M−f [ 0 ] ˆ 2 ; S [ 2 ] = V[ 2 ] /V [ 0 ] ; # c i t l i v o s t parametru X[ 3 ] for i to M do Ys4 [ i ] := s i n ( X1s2 [ i ])+7∗ s i n ( X2s2 [ i ])ˆ2+(1/10)∗ X3s [ i ]ˆ4∗ s i n ( X1s2 [ i ] ) end do : V[ 3 ] := add ( Ys [ i ] ∗ Ys4 [ i ] , i = 1 . . M)/M−f [ 0 ] ˆ 2 ; S [ 3 ] = V[ 3 ] /V [ 0 ] ; # k o r e k c n i c l e n for i to M do Ysk [ i ] := s i n ( X1s2 [ i ])+7∗ s i n ( X2s2 [ i ])ˆ2+(1/10)∗ X3s2 [ i ]ˆ4∗ s i n ( X1s2 [ i ] ) end do : V[ k ] := add ( Ys [ i ] ∗ Ysk [ i ] , i = 1 . . M)/M−f [ 0 ] ˆ 2 ; Matlab M=100000; % M r e a l i z a c i j e d n o t l i v y c h v e l i c i n X1s = u n i f r n d (−pi , pi , [M 1 ] ) ; X2s = u n i f r n d (−pi , pi , [M 1 ] ) ; X3s = u n i f r n d (−pi , pi , [M 1 ] ) ; Ys = sin ( X1s ) + 7∗ sin ( X2s ) . ˆ 2 + (1/10)∗ X3s . ˆ 4 . ∗ sin ( X1s ) ; mean Ys = mean( Ys ) ; % odhad s t r e d n i hodnoty var Ys = var ( Ys ) ; % odhad r o z p t y l u % ” druha s k u p i n a ” M r e a l i z a c i j e d n o t l i v y c h v e l i c i n X1s2 = u n i f r n d (−pi , pi , [M 1 ] ) ; X2s2 = u n i f r n d (−pi , pi , [M 1 ] ) ; X3s2 = u n i f r n d (−pi , pi , [M 1 ] ) ; % ” druha r e s e n i ” pro j e d n o t l i v e v l i v y Ys2 = sin ( X1s ) + 7∗ sin ( X2s2 ) . ˆ 2 + (1/10)∗ X3s2 . ˆ 4 . ∗ sin ( X1s ) ; Ys3 = sin ( X1s2 ) + 7∗ sin ( X2s ) . ˆ 2 + (1/10)∗ X3s2 . ˆ 4 . ∗ sin ( X1s2 ) ; Ys4 = sin ( X1s2 ) + 7∗ sin ( X2s2 ) . ˆ 2 + (1/10)∗ X3s . ˆ 4 . ∗ sin ( X1s2 ) ; Ysk = sin ( X1s2 ) + 7∗ sin ( X2s2 ) . ˆ 2 + (1/10)∗ X3s2 . ˆ 4 . ∗ sin ( X1s2 ) ; % c i t l i v o s t parametru X[ 1 ] V1 = sum( Ys . ∗ Ys2 )/M − mean Ys ˆ 2 ; S1 = V1/ var Ys % c i t l i v o s t parametru X[ 2 ] V2 = sum( Ys . ∗ Ys3 )/M − mean Ys ˆ 2 ; S2 = V2/ var Ys % c i t l i v o s t parametru X[ 3 ] V3 = sum( Ys . ∗ Ys4 )/M − mean Ys ˆ 2 ; S3 = V3/ var Ys % k o r e k c n i c l e n Vk = sum( Ys . ∗ Ysk )/M − mean Ys ˆ 2 ; S3 = Vk/ var Ys 1.4.4 Příklad 4 Jako závěrečný ilustrační příklad uvažujme chemickou reakci a A + b B r −→ c C, při níž z látek A a B vzniká látka C. Rychlost reakce r je definována jako (viz [7]): 25 r = − 1 a · d[A] dt = − 1 b · d[B] dt = 1 c · d[C] dt , kde symboly v hranatých závorkách značí koncentraci příslušné látky, a platí pro ni GuldbergůvWaagův zákon (viz [1]): r = k · [A]α · [B]β . (1.59) Konstanta k z rovnice (1.59) se nazývá rychlostní konstanta. Koeficienty α, β jsou obecně různé od stochiometrických koeficientů a, b v rovnici chemické reakce a musí být určeny experimentálně. Pro jednoduchost zvolme reakci A + A r −→ C. a předpokládejme, že α = 1, β = 1. Rychlostní konstanta je tepelně závislá a platí pro ni Arrheniova rovnice (viz [17]) k = X1 · e −X2 T . (1.60) Nechť v našem případě X1 ∼ U(8, 97·106 ; 3, 59·107 ) dm3 mol s−1 , X2 ∼ U(0, 1000) K−1 a T = 300 K. Pro počáteční koncentraci Y0 látky A v čase t = 0 s platí: Y0 ∼ U(1, 0; 1, 2) mol dm−3 . Pro hodnotu koncentrace Y (t) látky A v čase t dostáváme modelovou rovnici dY (t) dt = −2 · k · Y 2 (t), Y (0) = Y0. (1.61) Díky jednoduchosti modelové rovnice 1.61 jsme schopni určit její řešení, a získat tak model Y (t) = 2 · X1 · e −X2 T · t + 1 Y0 −1 . (1.62) Analýza neurčitosti Řešení modelu je závislé na čase. K prvotnímu nahlédnutí, jak se řešení chová, můžeme vykreslit grafy hraničních situací, tj. maximální Y (t) a minimální Y (t) pro daná t. Maximální Y (t) dostaneme maximalizací veličin X2 a Y0 a minimalizací veličiny X1 (viz předpis řešení (1.62)). Analogicky získáme minimální Y (t). Řešení se tak pro libovolné realizace veličin X1, X2, Y0 pohybuje mezi křivkami na obrázku 1.5.(a). Na tomtéž grafu jsou také odhady střední hodnoty (v některých časových bodech). Od počátku až do času t = 10−4 je zřejmý rozptyl hodnot způsobený neurčitostmi v parametrech X1, X2 a Y0. Dále od tohoto časového bodu převažuje vliv času t na hodnotu Y (t) a neurčitosti v parametrech tak mohou být zanedbány. Zaměříme se proto pouze na vyšetření intervalu t ∈ [10−10 , 10−4 ], jak je tomu i na obázku 1.5. Na grafu 1.5.(b) jsou zobrazeny odhady střední hodnoty a směrodatné odchylky řešení Y (t). Můžeme pozorovat, že největší směrodatné odchylky dosahuje Y (t) kolem času t = 10−7 . K tomuto grafu lze dospět následujícími zdrojovými kódy. 26 (a) Maximální a minimální řešení (1.62) (b) Odhad střední hodnoty a směrodatné odchylky řešení (1.62) Obrázek 1.5: Grafy analýzy neurčitosti ilustračního příkladu 4. Maple M := 10ˆ5: with ( S t a t i s t i c s ) : X[ 1 ] := RandomVariable ( Uniform ( 8 . 9 7 ∗ 1 0 ˆ 6 , 3 . 5 9 ∗ 1 0 ˆ 7 ) ) : X[ 2 ] := RandomVariable ( Uniform ( 0 , 1 0 0 0 ) ) : Y0 := RandomVariable ( Uniform ( 1 , 1 . 2 ) ) : X1s := Sample (X[ 1 ] , M) : # M r e a l i z a c i v e l i c i n y X[ 1 ] X2s := Sample (X[ 2 ] , M) : # M r e a l i z a c i v e l i c i n y X[ 2 ] Y0s := Sample (Y0 , M) : # M r e a l i z a c i v e l i c i n y X[ 3 ] Ys := Vector (M) : for i to M do Ys [ i ] := Y0s [ i ] / ( 2 ∗ X1s [ i ] ∗ exp ( −(1/300)∗ X2s [ i ] ) ∗ t ∗Y0s [ i ]+1) end do : # vypocet s t r . hodnot a smerodat . odchylek for i from 1 to 21 do t := 10ˆ( −10)∗2ˆ( i −1); mean Ys [ i ] := Mean( Ys ) ; dev Ys [ i ] := StandardDeviation ( Ys ) ; end do : t := ’ t ’ : # ” priprava ” bodu k v y k r e s l e n i str hodnoty := seq ([10ˆ( −10)∗2ˆ( i −1) , mean Ys [ i ] ] , i = 1 . . 2 1 ) ; odchylky := seq ([10ˆ( −10)∗2ˆ( i −1) , dev Ys [ i ] ] , i = 1 . . 2 1 ) ; # v y k r e s l e n i p l o t 1 := p l o t s [ p o i n t p l o t ] ( [ str hodnoty ] ) ; p l o t 2 := p l o t s [ p o i n t p l o t ] ( [ odchylky ] , symbol = c r o s s ) ; p l o t s [ d i s p l a y ] ( { p l o t 1 , p l o t 2 , p l o t s [ s e m i l o g p l o t ] ( 0 , view = [10ˆ( −10) . . 3∗10ˆ( −4) , 0 . . 1 . 2 ] , l a b e l s = [ t y p e s e t ( t ∗ [ s ] ) , t y p e s e t (Y( t ) ∗ [ mol/dm ˆ 3 ] ) ] ) } ) ; Matlab M=100000; % M r e a l i z a c i j e d n o t l i v y c h v e l i c i n X1s = u n i f r n d ( 8 . 9 7 ∗ 1 0 ˆ 6 , 3.59∗10ˆ7 , [M 1 ] ) ; X2s = u n i f r n d ( 0 , 1000 , [M 1 ] ) ; Y0s = u n i f r n d ( 1 , 1 . 2 , [M 1 ] ) ; mean Ys=zeros ( 1 , 2 1 , ’ double ’ ) ’ ; dev Ys=zeros ( 1 , 2 1 , ’ double ’ ) ’ ; time=zeros ( 1 , 2 1 , ’ double ’ ) ’ ; % v y p o c e t s t r e d n i c h hodnot a smerodatnych o d c h y l e k for i =1:21 t ( i ) = 10ˆ( −10)∗2ˆ( i −1); Ys = Y0s . / ( 2 ∗ X1s . ∗ exp( −(1/300)∗ X2s )∗ t ( i ) . ∗ Y0s +1); mean Ys ( i ) = mean( Ys ) ; dev Ys ( i ) = std ( Ys ) ; end semilogx ( t , mean Ys , ’ . ’ ) ; % v y k r e s l e n i s t r . hodnot hold on semilogx ( t , dev Ys , ’+ ’ ) ; % v y k r e s l e n i s t d . o d c h y l e k xlabel ( ’ t [ s ] ’ ) ; ylabel ( ’Y( t ) [ mol/dmˆ 3 ] ’ ) ; 27 Lokální analýza citlivosti Lokální citlivosti bychom získali přímo z definice díky předpisu řešení (1.62). Tentokrát jsou citlivosti závislé nejen na všech parametrech X1, X2, Y0, ale též na čase t. Globální analýza citlivosti Analýzu globální citlivosti je třeba vyšetřit opět v závislosti na čase. Vzhledem k tomu, že se nám nepodaří určit globální citlivosti analyticky, musíme je počítat numericky pro konkrétní časové body. Aplikujme tentokrát metodu FAST konkrétně na bod t = 10−7 , kde dosahuje směrodatná odchylka (a potažmo i rozptyl) řešení nejvyšších hodnot. Nejprve je třeba zvolit transformační funkce Gi a úhlové frekvence ωi (pro i = 1, 2, 3). Využijme dříve uvedené funkce Gi(s) = 1 2 + 1 π · arcsin(sin(ωi · s)). Jelikož tato funkce trasformuje interval (−π, π) na interval (0, 1), musíme ji pro naše potřeby upravit. Všechny náhodné veličiny v tomto příkladu mají rovnoměrné rozdělení, takže stačí zapsanou funkci Gi vynásobit délkou příslušného intervalu neurčitosti a přičíst dolní mez tohoto intervalu. Dostaneme tedy G1(s) = 2, 693 · 107 · 1 2 + 1 π · arcsin(sin(ω1 · s)) + 8, 97 · 106 , G2(s) = 1000 · 1 2 + 1 π · arcsin(sin(ω2 · s)) , G3(s) = 0, 2 · 1 2 + 1 π · arcsin(sin(ω3 · s)) + 1. Úhlové frekvence ωi nechť splňují ω1 = 11, ω2 = 21, ω3 = 31. Transformovali jsme tak modelovou rovnici (1.62) na rovnici Y (t) = 2 · G1(S) · e −G2(S) T · t + 1 G3(S) −1 , (1.63) kde S ∼ U(−π, π). Počet generovaných bodů ponechme stejný, tj. M = 105 . Nechť dále počet harmonických složek N = 6. Získáme tak: E Y 10−7 . = 0, 5801; V Y 10−7 . = 0, 0599; V1 . = 0, 0070; V2 . = 0, 0513; V3 . = 0, 0008; S1 . = 0, 1177; S2 . = 0, 8573; S3 . = 0.0128. Pro výpočet lze použít zobrazené zdrojové kódy. 28 Maple M := 10ˆ5: N:=6: with ( S t a t i s t i c s ) : S := RandomVariable ( Uniform(−Pi , Pi ) ) : Ss := Sample (S , M) : # t r a n s f o r m a c n i funkce G[ 1 ] ( s ):=2.693∗10ˆ(7)∗(1/(2)+1/( Pi )∗ a r c s i n ( s i n ( omega [ 1 ] ∗ s ) ) ) + 8 . 9 7 ∗ 1 0 ˆ ( 6 ) : G[ 2 ] ( s ):=1000∗(1/(2)+1/( Pi )∗ a r c s i n ( s i n ( omega [ 2 ] ∗ s ) ) ) : G[ 3 ] ( s ):=0.2∗(1/(2)+1/( Pi )∗ a r c s i n ( s i n ( omega [ 3 ] ∗ s ) ) ) + 1 : # uhlove f r e k v e n c e omega [ 1 ] := 1 1 ; omega [ 2 ] := 2 1 ; omega [ 3 ] := 3 1 : Ys := Vector (M) : for i from 1 to M do Ys [ i ] := e v a l f (G[ 3 ] ( Ss [ i ] ) / ( 2 ∗G[ 1 ] ( Ss [ i ] ) ∗ exp ( −(1/300)∗G[ 2 ] ( Ss [ i ]))∗10ˆ( −7)∗G[ 3 ] ( Ss [ i ] ) + 1 ) ) end do : f [ 0 ] := Mean( Ys ) ; # odhad s t r e d n i hodnoty V[ 0 ] := Variance ( Ys ) ; # odhad r o z p t y l u # vypocet Fourierovych k o e f i c i e n t u for j from 1 to N do fk [ j ] := add ( Ys [ i ] ∗ cos ( j ∗omega [ 1 ] ∗ Ss [ i ] ) , i = 1 . . M)/M: fk [ j+N] := add ( Ys [ i ] ∗ s i n ( j ∗omega [ 1 ] ∗ Ss [ i ] ) , i = 1 . . M)/M: fk2 [ j ] := add ( Ys [ i ] ∗ cos ( j ∗omega [ 2 ] ∗ Ss [ i ] ) , i = 1 . . M)/M: fk2 [ j+N] := add ( Ys [ i ] ∗ s i n ( j ∗omega [ 2 ] ∗ Ss [ i ] ) , i = 1 . . M)/M: fk3 [ j ] := add ( Ys [ i ] ∗ cos ( j ∗omega [ 3 ] ∗ Ss [ i ] ) , i = 1 . . M)/M: fk3 [ j+N] := add ( Ys [ i ] ∗ s i n ( j ∗omega [ 3 ] ∗ Ss [ i ] ) , i = 1 . . M)/M: end do : # c i t l i v o s t parametru X[ 1 ] V[ 1 ] := 2∗add ( fk [ i ] ˆ 2 , i = 1 . . 1 2 ) ; S [ 1 ] = V[ 1 ] /V [ 0 ] ; # c i t l i v o s t parametru X[ 2 ] V[ 2 ] := 2∗add ( fk2 [ i ] ˆ 2 , i = 1 . . 1 2 ) ; S [ 2 ] = V[ 2 ] /V [ 0 ] ; # c i t l i v o s t parametru Y[ 0 ] V[ 3 ] := 2∗add ( fk3 [ i ] ˆ 2 , i = 1 . . 1 2 ) ; S [ 3 ] = V[ 3 ] /V [ 0 ] ; Matlab M=100000; N=6; omega1=11; omega2=21; omega3=31; % u h l o v e f r e k v e n c e Ss = u n i f r n d (−pi , pi , [M 1 ] ) ; x1 = 2 . 6 9 3 ∗ 1 0 ˆ 7 ∗ ( . 5 + asin ( sin ( omega1∗Ss ) ) / pi )+8.97∗10ˆ6; x2 = 1000∗(.5 + asin ( sin ( omega2∗Ss ) ) / pi ) ; y0 = . 2 ∗ ( . 5 + asin ( sin ( omega3∗Ss ) ) / pi )+1; % modelova r o v n i c e Y = y0 . / ( 2 ∗ x1 . ∗ exp(−x2 /300)∗10ˆ( −7).∗ y0 + 1 ) ; mean Y = mean(Y) % odhad s t r e d n i hodnoty var Y = var (Y) % odhad r o z p t y l u fk=zeros (1 ,2∗N, ’ double ’ ) ’ ; fk2=zeros (1 ,2∗N, ’ double ’ ) ’ ; fk3=zeros (1 ,2∗N, ’ double ’ ) ’ ; % v y p o c e t F o u r i e r o v y c h k o e f i c i e n t u for j =1:N fk ( j ) = sum(Y. ∗ cos ( j ∗omega1∗Ss ) ) /M; fk ( j+N) = sum(Y. ∗ sin ( j ∗omega1∗Ss ) ) /M; fk2 ( j ) = sum(Y. ∗ cos ( j ∗omega2∗Ss ) ) /M; fk2 ( j+N) = sum(Y. ∗ sin ( j ∗omega2∗Ss ) ) /M; fk3 ( j ) = sum(Y. ∗ cos ( j ∗omega3∗Ss ) ) /M; fk3 ( j+N) = sum(Y. ∗ sin ( j ∗omega3∗Ss ) ) /M; end % c i t l i v o s t parametru X[ 1 ] V1 = 2∗sum( fk . ˆ 2 ) S1 = V1/ var Y % c i t l i v o s t parametru X[ 2 ] V2 = 2∗sum( fk2 . ˆ 2 ) S2 = V2/ var Y % c i t l i v o s t parametru Y[ 0 ] V3 = 2∗sum( fk3 . ˆ 2 ) S3 = V3/ var Y 29 1.5 Příklady k procvičení Příklad 1.6: Uvažujme model Y = X1 + X2, kde X1 ∼ N(1, 4), X2 ∼ N(2, 9). Určete E(Y ), V (Y ) a indexy globální citlivosti prvního řádu.11 Příklad 1.7: Nechť Y = 2X1 · X2, kde X1 ∼ U(0, 3), X2 ∼ U(0, 3). Proveďte analýzu neurčitosti i analýzu citlivosti řešení uvedeného modelu, tj. určete E(Y ), V (Y ), p(y) a indexy globální i lokální citlivosti prvního řádu. Příklad 1.8: Uvažujme model Y = X4 1 X2 2 , kde X1 ∼ U(1 2 , 3 2 ), X2 ∼ U(1 2 , 3 2 ). Určete E(Y ), V (Y ) a indexy globální i lokální citlivosti prvního řádu. Příklad 1.9: Nechť Y = 8 i=1 Xi, kde X1 ∼ N(0, 1), X2 ∼ N(2, 4), X3 ∼ N(1, 9), X4 ∼ N(1, 25), X5 ∼ N(3, 1), X6 ∼ N(4, 1), X7 ∼ N(1, 4), X8 ∼ N(5, 25). Určete, které tři parametry nejméně ovlivňují rozptyl řešení Y . Příklad 1.10: Uvažujme model Y = sin(X1) + 7 · sin2 (X2) + 1 10 · X3 · sin(X1) z třetího ilustračního příkladu se stejně rozdělenými parametry X1, X2, X3. Vypočítejte indexy globální citlivosti prvního řádu pomocí metody FAST. Příklad 1.11: Vraťte se k modelu chemické rovnice ze čtvrtého ilustračního příkladu. Vypočítejte indexy globální citlivosti prvního řádu pomocí Sobol’ovy metody. Dále vytvořte lineární modely (a to jak „klasické , tak pro přiřazené pořadí) pro různé časové body v intervalu [10−10 , 10−4 ] a zkoumejte závislost koeficientu determinace na čase. 11 N(µ, σ2 ) značí normální rozdělení se střední hodnotou µ a rozptylem σ2 . 30 Literatura [1] Burda, J.: Chemická kinetika [online]. [cit 2011-09-25]. Dostupný z WWW: . [2] Drosg, M.: Dealing with uncertainties. Springer (2009) [3] Homma, T., Saltelli, A.: Importance measures in global sensitivity analysis of nonlinear models. Reliability Engineering and System Sa]ety (1996) [4] Horová, I., Zelinka, J.: Numerické metody. Masarykova univerzita, Brno (2004) [5] Hřebíček, J., Pospíšil, Z., Urbánek, J.: Úvod do matematického modelování s využitím Maple. CERM, Brno (2010) [6] Hušková, M.: Bayesovské metody. Karlova univerzita, Praha (1985) [7] IUPAC: Gold Book [online]. [cit 2011-09-25]. Dostupný z WWW: . [8] Komárek, A.: Bayesovské metody (druhá část) [online]. [cit 2011-09-01]. Dostupný z WWW: . [9] Luo, Y., Yang, X.: A multimedia environmental model of chemical distribution: Fate, transport, and uncertainty analysis Chemosphere, 66 (8), p. 1396 – 1407 (2007) [10] Muehlenstaedt, T., Roustant, 0., Carraro, L., Kuhnt, S.: Data-driven Kriging models based on FANOVA-decomposition [online]. [cit 2011-09-09]. Dostupný z WWW: . [11] Qian, S. S., Stow, C. A., Borsuk, M. E.: On Monte Carlo methods for Bayesian inference [online]. [cit 2011-09-01]. Dostupný z WWW: . [12] Roustant, 0.: Analytical computation of Sobol indices with a Gaussian process metamodel [online]. [cit 2011-09-09]. Dostupný z WWW: . [13] Saltelli, A., Chan, K., Scott, E. M.: Sensitivity analysis. Wiley (2000) [14] Saltelli, A., Ratto, M., Andres, T., Campolongo, F., Cariboni, J., Gatelli, D., Saisana, M., Tarantola, S.: Global sensitivity analysis. Wiley (2008) 31 [15] Saltelli, A., Tarantola, S., Campolongo, F.: Sensitivity Analysis as an Ingredient of Modeling. Statistical Science (2000) [16] Shonkwiler, R. W., Herod, J.: Mathematical Biology. An Introduction with Maple and Matlab. Springer (2009) [17] Wikipedia: Arrhenius equation [online]. [cit 2011-09-25]. Dostupný z WWW: . [18] Wikipedia: Shannonův teorém [online]. [cit 2011-09-18]. Dostupný z WWW: . [19] Wikipedia: Spearman’s rank correlation coefficient [online]. [cit 2011-09-08]. Dostupný z WWW: . 32