Jak na nejistoty metodou Monte Carlo jednoduše a bez drahých programů

Rok / Year: 2014

Svazek / Volume: 16

Číslo / Number: 2

Jak na nejistoty metodou Monte Carlo jednoduše a bez drahých programů How to calculate uncertainties by means of Monte Carlo method in an easy manner and without any expensive software Martin Šíra [email protected] Český Metrologický Insitut, Brno

Abstrakt: Dokument se zabývá výpočtem nejistoty měřené veličiny metodou Monte Carlo z praktického pohledu. Jsou řešeny hlavně problémy výběru vhodného programu, generování náhodných čísel požadovaných rozdělení a zpracování výsledků. Součástí článku jsou dva řešené příklady. Článek předpokládá alespoň základní znalosti v oblasti metrologie a výpočtu nejistot.

Abstract: A calculation of uncertainties by means of Monte Carlo method is explained in the presented paper. Selection of the computer software, generation of random numbers of required probability distributions and evaluation of results is discussed. A basic knowledge of metrology and uncertainty calculation is required.

VOL.16, NO.2, APRIL 2014

Jak na nejistoty metodou Monte Carlo jednoduše a bez drahých programů Martin Šíra Český Metrologický Insitut, Brno Email: [email protected]

2

Abstrakt – Dokument se zabývá výpočtem nejistoty měřené veličiny metodou Monte Carlo z praktického pohledu. Jsou řešeny hlavně problémy výběru vhodného programu, generování náhodných čísel požadovaných rozdělení a zpracování výsledků. Součástí článku jsou dva řešené příklady. Článek předpokládá alespoň základní znalosti v oblasti metrologie a výpočtu nejistot.

1

Výběr programu

Nejprve je třeba zvolit si vhodný program. Většina by jako první zkusila Microsoft Excel. Bohužel do verze 2007 měl generátor pseudonáhodných čísel použitý v Excelu vady, které jej diskreditují [5]. Poté již byly chyby z větší části odstraněny, ale stále generátor nedosahuje kvalit specializovaných programů. Také limit počtu řádků starších verzí Excelu na 65 536 je pro MMC omezující, a je vhodné použít vyšší verze, kde je limit přibližně 106 řádků. Přesto bych ani nejnovější verze Excelu nedoporučil, neboť práce s ním je těžkopádná. Například vytvořit dynamicky počítaný histogram není jednoduchý úkol. Vypočítat interval nejmenšího pokrytí pravděpodobnostní funkce výstupní veličiny je již komplikované. Excel lze tedy použít pro vyzkoušení MMC, obtížnější ale bude aplikovat ho pro smysluplný výpočet z praxe. Přesto z důvodu rozšířenosti programu v článku popíšu, jak použít Excel 2010. Obdoba výpočtu bude použitelná i v programech LibreOffice nebo OpenOffice.

Úvod

Ve známém dokumentu Guide to The Expression of Uncertainty In Measurement (GUM) [1] je popsána metoda na výpočet nejistot, často nazývaná Gum Uncertainty Framework (GUF). Bohužel tato metoda má mnoho omezení, a proto byl vydán první dodatek [2], ve kterém je uveden výpočet nejistot pomocí metody Monte Carlo (MMC). Postup byl dále rozšířen na více výstupních veličin v druhém dodatku [3]. MMC je velmi jednoduchý numerický algoritmus, jehož nevýhodou je, že už nestačí tužka a papír. Je potřeba vhodný program s kvalitním generátorem čísel, který dokáže zpracovat desítky tisíc až milióny čísel. Oproti tomu pomocí MMC lze numericky vyřešit i nejistoty, které nejsou řešitelné analyticky nebo metodou GUF. Není třeba počítat parciální derivace a řešit stupně volnosti. Metoda GUF je založena na mnoha aproximacích, jejichž vhodnost se obvykle kontroluje právě pomocí MMC. Například již v tak základním případu, jako je výpočet nejistoty poměru dvou veličin [4], může při použití GUF dojít k příliš velkým aproximacím, a je třeba použít MMC (nebo jinou metodu).

Velmi vhodný program pro MMC je například Matlab [6]. Pro svou jednoduchou syntaxi, velký výpočetní výkon a komplexní grafické uživatelské prostředí je takřka ideálním nástrojem. Jeho nevýhoda je ovšem cena. Základní komerční verze stojí více než 50 000 Kč. Ale naštěstí existuje celý svět programů s otevřeným zdrojovým kódem. Například GNU Octave [7] je program zdarma, se syntaxí a použitím téměř identickým jako Matlab, a výpočetní výkon je v mnoha oblastech s Matlabem srovnatelný. Neobsahuje ovšem natolik komfortní uživatelské prostředí jako Matlab. Proto v článku podrobně popíšu i jak nainstalovat a použít GNU Octave. Veškeré uvedené příklady ale budou plně funkční i v Matlabu.

Požadavky na programové vybavení odradí od MMC spoustu metrologů. Proto jsem se snažil sepsat velmi jednoduchý návod, jak na MMC pomocí snadno dostupných programů. V článku budu předpokládat základní znalosti z oblasti nejistot a práce s počítačem.

Dále existuje spousta dalších programů více či méně vhodných pro výpočty MMC. Ty více vhodné budou například statistický program R (zdarma), méně vhodný bude celkem oblíbený program na zpracování dat Origin (placený). Obdobné programy již nechám na zkušenosti čtenáře.

Snažil jsem se napsat článek přístupnou formou, proto není metoda MMC popsána exaktně a úplně a dopouštím se mnoha zjednodušení. Pro přesný popis a definice je nejlepší prostudovat výše zmíněné dodatky.

Další kapitoly se budou věnovat MMC a výpočtům pomocí této metody. Pokud si vyberete jen jeden program, ve kterém chcete MMC vyzkoušet, můžete přeskočit kapitoly, které se ho netýkají.

V článku jsou různým písmem rozlišeny veličiny matematických vzorců a proměnné programů. Tedy veličina je zapsána takto: X, ale proměnná takto: X.

86


Obrázek 1: Instalace GNU Octave. Doporučuji zvolit rozhraní pro vykreslování grafů Gnuplot.

3

Obrázek 2: Instalace GNU Octave. Doporučuji zvolit instalaci do adresáře neobsahující mezery v názvu.

Příprava pro použití MMC: Instalace výstupní veličiny. Podobně jako u metody GUF, je třeba před samotným výpočtem nejistot určit následující: GNU Octave 1. vstupní veličiny Xi ,

Program GNU Octave je dostupný na stránkách SourceForge [8]. Jelikož je možné stáhnout mnoho verzí a program se neustále vyvíjí, vybral jsem jednu konkrétní verzi a umístil ji na své stránky [9]. Veškerý popis práce se bude vztahovat k této jediné verzi. Není třeba bát se licencí, program je zdarma i pro komerční použití. Zadejte adresu internetových stránek do internetového prohlížeče a zvolte odkaz instalačního souboru. Proveďte instalaci spuštěním staženého souboru. Instalace se provede obvyklou opakovanou volbou tlačítka Next. Nejprve je třeba zvolit vlastnosti CPU, obvykle není třeba nic měnit a pokračujte volbou tlačítka Next. Pak se volí součásti, které budou nainstalovány. Opět není třeba nic měnit. Další volba se týká grafického rozhraní pro vykreslování grafů, doporučuji zvolit položku Gnuplot, viz obr. 1. V dalším kroku se volí, kde budou soubory programu nainstalovány. Doporučuji nainstalovat GNU Octave do adresáře bez mezer (mohli byste narazit na omezení dané zpětnou kompatibilitou platformy Windows), viz obr. 2. Po volbě Next a Install se provede instalace programu. Hned po instalaci se GNU Octave spustí. Pokud ne, spusťte jej z nabídky Start operačního systému. Uvidíte příkazovou řádku (obdobně jako v Matlabu). Zkusit vykreslit graf můžete zadáním následujícího příkazu:

2. výstupní veličinu Y , 3. model měření Y = f (X1 , X2 , . . . , Xi ). U metody GUF dále je třeba znát střední hodnoty vstupních veličin xi , jejich nejistoty u(xi ) a stupně volnosti ν(xi ). Výsledkem GUF je střední hodnota y, nejistota uc (y) a stupeň volnosti ν(y). Naproti tomu MMC vyžaduje znalost hustot pravděpodobnosti vstupních veličin, a výsledkem je pravděpodobnostní funkce výstupní veličiny (viz obr. 4). Algoritmus MMC (viz obr. 5) je následující: 1. Zvolí se dostatečně velký počet opakování N . 2. Pro každou vstupní veličinu Xi je vygenerováno pseudonáhodné číslo podle její hustoty pravděpodobnosti. 3. Vygenerovaná pseudonáhodná čísla jsou dosazena do modelu měření. 4. Body 2 a 3 jsou opakovány N krát. 5. Výsledkem MMC je N vypočtených hodnot výstupní veličiny Y tvořící pravděpodobnostní funkci výstupní veličiny.

plot ([1 2 3]) 6. Pokud je rozdělení výstupní veličiny Y normální, nejistota u(y) je vypočtena jako výběrová směrodatná odchylka: s PN 2 i=1 (yi − y) (1) u(Y ) = N −1

a potvrďte klávesou Enter. Vykreslení prvního grafu po instalaci může trvat delší dobu. Pokud vidíte graf s jednou čarou (obr. 3), instalace proběhla v pořádku.

4

Krátký úvod do metody Monte Carlo

Princip MMC je generování pseudonáhodných čísel podle hustoty pravděpodobnosti vstupních veličin, jejich zadání do modelu měření a vypočtení pravděpodobnostní funkce

Jinak je nejistota určena hledáním nejkratšího intervalu pokrytí dle zvolené pravděpodobnosti pokrytí.

87


Obrázek 3: Příkazová řádka se spuštěným programem Octave a vykresleným grafem.

Obrázek 4: Porovnání vstupních a výstupních údajů metod GUF a MMC. Obrázek 5: Algoritmus výpočtu nejistot metodou Monte Carlo.

88


4.1

zvolit N ≥ 104 . Dovolím si říct, že i pro nejkomplikovanější výpočty nejistot je N = 106 naprosto dostatečné. Pro komplikované výpočty může být účelné snížit N na minimum. V tom případě lze použít adaptivní metodu určení N , kdy se sleduje konvergence nejistoty výstupní veličiny Y . Adaptivní metoda je popsaná v prvním dodatku GUMu [2].

Určení hustot pravděpodobnosti vstupních veličin

Jak zjistíme hustoty pravděpodobnosti vstupních veličin? Postupujeme podobně jako v metodě GUF. Někdy se může stát, že jsou přímo známy. V opačném případě jsou nejobvyklejší případy: Normální rozdělení Pokud je vstupní veličina a její nejistota vyčtena například z kalibračního listu a nejsou přítomny doplňující informace, můžeme předpokládat, že hodnoty byly zpracovány podle správných metrologických postupů, a tedy vstupní veličina má normální rozdělení se střední hodnotou rovnou hodnotě uvedené v kalibračním listu a směrodatnou odchylkou danou standardní nejistotou dle kalibračního listu (pozor, v kalibračních listech se obvykle uvádí nejistota s pravděpodobností pokrytí 95,45 % (koeficient rozšíření k = 2, tedy standardní nejistota je rovna polovině uvedené nejistoty). Většina běžně měřených veličin má rovnoměrné rozdělení.

4.3

Následuje přehled rovnic buněk pro vygenerování pseudonáhodných čísel o očekávané hodnotě pX a standardní nejistotě uX, s případným stupněm volnosti nX, nebo mezemi a, b (přičemž a
Studentovo rozdělení Pokud vstupní veličina pochází z opakovaného měření, předpokládáme Studentovo rozdělení se střední hodnotou danou průměrem naměřených hodnot, se stupněm volnosti rovném počtu měření mínus jedna a škálované nejistotou vypočtenou jako nejistota stanovovaná způsobem A. Také některé kalibrační listy uvádí podrobnější informace, jako je efektivní stupeň volnosti a koeficient rozšíření, tedy opět předpokládáme Studentovo rozdělení škálované nejistotou uvedenou v kalibračním listu. Jako stupeň volnosti Studentova rozdělení se použije efektivní stupeň volnosti dle kalibračního listu.

Studentovo rozdělení =T.INV(NÁHČÍSLO();nX)*uX+pX Rovnoměrné rozdělení =(NÁHČÍSLO())*(b-a)+a Trojúhelníkové rozdělení Do buňky A1 je třeba zadat vzorec: =(pX-a)/(b-a) Do sloupce B zadat: =NÁHČÍSLO() a do sloupce C zadat ($A$1 je odkaz na buňku sloupce A, B1 je odkaz na buňky sloupce B): =KDYŽ(B1<$A$1;a+ODMOCNINA(B1*(b-a)*(pX-a)); b-ODMOCNINA((1-B1)*(b-a)*(b-pX)))

Rovnoměrné rozdělení Pokud u vstupní veličiny známe pouze maximální a minimální hodnotu, mezi kterými předpokládáme stejnou pravděpodobnost veličiny, předpokládáme rovnoměrné rozdělení s mezemi danými minimální a maximální hodnotou. Příkladem je třeba nejistota způsobená zaokrouhlením na poslední platnou číslici při zobrazení měřené hodnoty na číslicovém přístroji, aproximace nejistoty odhadem apod.

Ukázkový soubor pro Excel s vygenerovanými náhodnými čísly podle různých rozdělení si můžete stáhnout na mých stránkách [9] 4.4

Generování pseudonáhodných čísel dle požadovaného rozdělení v GNU Octave a Matlabu

Následuje přehled příkazů pro vygenerování N pseudonáhodných čísel X o očekávané hodnotě pX a standardní nejistotě uX, s případným stupněm volnosti nX, nebo mezemi a, b (přičemž a
Trojúhelníkové rozdělení Pokud u vstupní veličiny známe nejpravděpodobnější hodnotu a meze, mimo které je pravděpodobnost vstupní veličiny nulová, můžeme předpokládat trojúhelníkové rozdělení. Příkladem z praxe může být regulátor udržující veličinu v daných mezích a nejpravděpodobněji na zadané hodnotě, stabilita přístrojů mezi kalibracemi apod.

Normální rozdělení X=normrnd(pX,uX,1,N); nebo X=randn(1,N).*uX+pX;

Existuje mnoho dalších rozdělení pravděpodobnosti používaných pro vstupní veličiny, které jsou již mimo rozsah tohoto článku. 4.2

Generování pseudonáhodných čísel dle požadovaného rozdělení v Excelu

Studentovo rozdělení X=trnd(nX,1,N).*uX+pX; Rovnoměrné rozdělení X=unifrnd(a,b,1,N) nebo X=rand(1,N).*(b-a)+a;

Určení počtu opakování N

Jak určit počet opakování N ? Jelikož se nejistoty uvádí maximálně na dvě platná místa, ve většině případů stačí

89


Trojúhelníkové rozdělení u=rand(1,N); lim=(pX-a)/(b-a); X=[a+sqrt(u(u=lim))*(b-a)*(b-pX))];

snadno pro každou hodnotu yk určit všechny intervaly pokrytí odpovídající požadované pravděpodobnosti pokrytí a z nich zvolit nejkratší interval pokrytí, jehož meze odpovídají mezím nejistoty výsledné veličiny Y . Graficky je postup znázorněn na obr. 6 a 7. Zvolí se yk , které nyní slouží jako „leváÿ mez intervalu pokrytí yL k . Dle KDF se určí odpovídající „levostrannáÿ pravděpodobnost pL k , k ní se přičte p a tak se získá „pravostrannáÿ pravděpodobnost pP k . Opět pomocí křivky KDF se z pP k určí „praváÿ hodnota intervalu pokrytí yP k . Tento postup provedeme pro všechny yk . Tím získáme sadu intervalů pokrytí hyL k , yP k i pro požadovanou pravděpodobnost pokrytí p. Závislost délky intervalů pokrytí

Ověřit výsledné rozdělení lze například vykreslením histogramu pomocí příkazu: hist(X,50) 4.5

Určení nejistoty z pravděpodobnostní funkce výstupní veličiny

Jak se určí nejistota výstupní veličiny ze sady N čísel výstupní veličiny vypočtených pomocí MMC? Ve většině případů výsledná pravděpodobnostní funkce odpovídá normálnímu rozdělení. V takovém případě je standardní nejistota pro pravděpodobnost pokrytí 68,27 % rovna výběrové směrodatné odchylce, neboli: s PN 2 i=1 (yi − y) . (2) u(y) = N −1

lk = yP k − yL k

na „levostrannéÿ pravděpodobnosti pL k je pro normální rozdělení zobrazena na obr. 6. Jako výsledný interval pokrytí se zvolí zvolí nejkratší interval pokrytí. Pokud existuje více nejkratších intervalů, zvolí se ten, který je symetrický vůči očekávané hodnotě. Tímto způsobem získáme „levéÿ a „pravéÿ meze nejkratšího intervalu pokrytí pro zvolenou pravděpodobnost pokrytí, které odpovídají „levéÿ a „pravéÿ mezi nejistoty. Pokud jsou meze symetrické, můžeme, jak je běžné, uvádět jen jednu hodnotu. Často ale můžeme narazit na nesymetrická rozdělení, pro které intervaly pokrytí a tedy ani nejistoty výstupní veličiny nejsou symetrické. To je třeba případ, kdy výstupní veličina má rozdělení gama, jako na obr. 7.

Pro normální rozdělení také platí, že rozšířená nejistota pro pravděpodobnost pokrytí 95,45 % (neboli k = 2 v metodě GUF) je dvojnásobkem standardní nejistoty (k = 1 v metodě GUF). V případě složitějších modelů měření může být rozdělení pravděpodobnosti výstupní veličiny rozdílné od normálního rozdělení. V tom případě je třeba vypočítat nejkratší interval pokrytí s požadovanou pravděpodobností pokrytí p (např. p = 68,27 %, p = 95,45 % atd.). Nejprve se z výsledku MMC, tedy vektoru hodnot (y1 , . . . , yi , . . . , yN ), určí kumulativní distribuční funkce (KDF) výstupní náhodné veličiny, což je funkce udávající pro každou hodnotu yk pravděpodobnost, že náhodná veličina Y je menší nebo rovna yk , neboli: KDF(yk ) = P(Y ≤ yk ) .

5

N (yi < yk ) , N

První příklad

Jako první vyzkoušíme nejjednodušší příklad výpočtu nejistot. Mějme aditivní model měření, například dva odpory zapojené do série (nebo dvě koncové měrky za sebou) a zajímá nás celkový odpor (nebo celková délka), tedy součet dvou veličin, a nejistota součtu veličin. U každého odporu známe jeho hodnotu a nejistotu z kalibračního listu. Můžeme tedy předpokládat, že nejistoty mají normální rozdělení (jelikož kalibrační list neuvedl bližší informace). První odpor, neboli veličina X1 , má hodnotu 200 Ω a nejistotu udanou v kalibračním listě 50 Ω pro pravděpodobnost pokrytí 95,45 % (koeficient rozšíření k = 2), tedy standardní nejistota je 25 Ω (koeficient rozšíření k = 1). Druhý odpor, neboli veličina X2 , má hodnotu 300 Ω a nejistotu 60 Ω pro pravděpodobnost pokrytí 95,45 %, tedy standardní nejistota je 30 Ω. Zadání příkladu je shrnuto v tabulce 1. Model měření je tedy:

(3)

Nejprve se zvolí hodnoty yk jako vektor vzestupně seřazených rovnoměrně rozložených čísel v rozsahu výsledku MMC. Pro tyto hodnoty se vypočte KDF. Tedy pokud hodnoty výsledku MMC leží např. v intervalu od 4.5 do 6.92, je vhodné za hodnoty yk zvolit např. (4.50, 4.51, 4.52, . . . , 6.91, 6.92). Čím menší bude krok hodnot yk , tím menší bude výsledná numerická chyba, ale výpočet bude trvat déle. Poté se hodnoty yi se seřadí vzestupně a pro každou yk se vypočte poměr počtu hodnot yi menších než yk ku počtu všech hodnot yi , tedy: KDF(yk ) =

(5)

Y = X1 + X2 .

(6)

Spočítat výsledek a nejistotu Y metodou GUF je snadné. Očekávaná hodnota y je vypočtena následovně:

(4)

kde N (yi < yk ) je počet hodnot z výsledku MMC, které jsou menší než yk . Pro případ normálního rozdělení je výsledná kumulativní distribuční funkce zobrazena na obr. 6. Nyní lze

y = 200 + 300 = 500 .

(7)

Citlivostní koeficienty (parciální derivace modelu měření podle vstupních veličin) jsou rovny jedné, tedy výsledná

90


vstupní veličina X1 X2

Tabulka 1: Hodnoty vstupních veličin pro příklad č. 1 průměrná nejistota dle standardní rozdělení hodnota kalibračního listu nejistota pravděpodobnosti (k = 2) (k = 1) 200 50 25 normální 300 60 30 normální

nejistota s pravděpodobností pokrytí 68,27 % (koeficient rozšíření k = 1) je: p p u(y) = u(x1 )2 + u(x2 )2 = 252 + 302 = 39,1 (8) 5.1

vypočítat jejich standardní nejistotu (za předpokladu normálního rozdělení) pomocí funkce STDEVA. Tím se ujistíme, že pseudonáhodná čísla vygenerovaná podle tabulky 1 mají požadované parametry. Funkci STDEVA můžeme použít pro výpočet nejistoty výstupní veličiny Y , protože má normální rozdělení. Pro jiná rozdělení je situace komplikovanější. Normální rozdělení můžeme ověřit zobrazením histogramu vygenerovaných náhodných čísel, což odpovídá pravděpodobnostní funkci veličiny (správné by bylo použít statistický test, ale to je mimo rozsah tohoto článku). Histogram vytvoříme následovně. Zvolíme si šířku třídy, třeba 10, a zapíšeme tuto hodnotu do buňky I3. Poté do buňky J3 zapíšeme dolní mez první třídy: 370. Do buňky J4 vypočteme horní mez první třídy histogramu zadáním vzorce:

Výpočet příkladu 1 metodou MC v Excelu

Výpočet příkladu 1 v Excelu je pro snažší orientaci v buňkách zobrazen na obr. 8. Soubor se zpracovaným příkladem si můžete stáhnout z mých internetových stránek [9]. Nejprve je třeba vygenerovat pseudonáhodná čísla. Zvolme si počet opakování metody MC N = 10 000. Do buněk A2 až A10001 zadejte vzorec: =NORMINV(NÁHČÍSLO();200;25) kterým se vygenerují pseudonáhodná čísla vstupní veličiny X1 . Obdobně do buněk B2 až B10001 zadejte vzorec:

=J3+$I$3 Tento vzorec aplikujeme pro následující buňky J5 až J29. Tím jsme sestavili řadu čísel od 370 do 630 s krokem 10. Nyní vypočteme četnosti pro jednotlivé třídy. Je třeba označit buňky K3 až K29 a poté do řádku vzorců zadejte vzorec:

=NORMINV(NÁHČÍSLO();300;30) kterým se vygenerují pseudonáhodná čísla vstupní veličiny X2 . Funkce NÁHČÍSLO generuje pseudonáhodná čísla v intervalu h0, 1) s rovnoměrným rozdělením. Abychom získali pseudonáhodná čísla s normálním rozdělením, použijeme tato pseudonáhodná čísla jako vstupní parametr do inverzní funkce k distribuční funkci normálního rozdělení NORMINV. Druhý a třetí parametr této funkce je požadovaná střední hodnota normálního rozdělení a směrodatná odchylka. Hodnoty výstupní veličiny yi vypočteme dle modelu, tedy buňka C2 obsahuje vzorec:

=ČETNOSTI(C2:C10001;J3;J29) a zmáčkněte Ctrl+Shift+Enter (zadání maticové operace). Vzorec nesmíte zadat psaním přímo do buňky, a také před psaním vzorce do řádku vzorců musí být označeny všechny buňky K3 až K29, jinak výsledek nebude vypočten správně. Jako poslední krok vytvořte graf ze sloupců J a K. Měli byste získat křivku blízkou Gaussově křivce (což odpovídá pravděpodobnostní funkci normálního rozdělení), viz obr 8.

=A2+B2 Vzorec je třeba zkopírovat pro buňky C3 až C10001. Tím jsme získali 10 000 výsledných hodnot výstupní veličiny yi , které použijeme pro výpočet očekávané hodnoty a nejistoty. Do buňky G3 zadejte vzorec pro výpočet průměru:

5.2

Výpočet příkladu 1 metodou MC v GNU Octave a Matlabu

Celý výpočet příkladu 1 v GNU Octave nebo Matlabu je na 6 řádků příkazů i s vykreslením histogramu výstupní veličiny, neboli pravděpodobnostní funkci výstupní veličiny. Řádky s příkazy můžete zadat do GNU Octave nebo Matlabu postupně:

=PRŮMĚR(C2:C10001) a do buňky G8 zadejte vzorec pro výpočet výběrové směrodatné odchylky:

1

=STDEVA(C2:C10001) Buňka G3 tedy obsahuje očekávanou hodnotu y a buňka G8 obsahuje standardní nejistotu u(y) (za předpokladu normálního rozdělení) s pravděpodobností pokrytí 68,27 %. Pro kontrolu je dobré do buněk E3 a F3 zadat výpočet průměrných hodnot veličin X1 a X2 a do buněk E8 a F8

2 3 4 5 6

91

N =1 e6 ; pX1 =200; uX1 =25; pX2 =300; uX2 =30; X1 = normrnd ( pX1 , uX1 ,1 , N ) ; X2 = normrnd ( pX2 , uX2 ,1 , N ) ; Y = X1 + X2 ; hist (Y ,50) uY = std ( Y )


100

140

y yL

120

yP

yP

100

60

ˇcetnost

ˇcetnost

80

y

yL

40

80 60 40

20

20 0

-4

-3

-2

-1

0 yi

1

2

3

0

4

2

4

6

8

10

12

14

yi

1

1 pP

0.8

0.8 p (68,27%)

P(Y ≤ yk )

P(Y ≤ yk )

0

0.6 0.4 0.2 pL

0.6 p (68,27%)

0.4 0.2

0 -4

pP

-3

-2

yL -1

0 yk

yP 1 2

0 3

4

pL yL 0 2

yP 4

6

8

10

12

14

yk

4.5

25

4 20

3

15 lk

lk

3.5

2.5

10

2 1.5 1 -0.05

5 0

0.05 0.1 0.15 0.2 0.25 0.3 0.35 pL k

-0.05

Obrázek 6: Určení nejkratšího intervalu pokrytí veličiny Y reprezentované 10 000 čísly s normálním rozdělením. Nahoře: Histogram. Uprostřed: Kumulativní distribuční funkce. Dole: Závislost délky intervalu pokrytí na „levostrannéÿ pravděpodobnosti. Čárkovaně je vyznačen průměr dat, čerchovaně jsou vyznačeny meze nejkratšího intervalu pokrytí. Hvězdičkou je označen nejkratší interval pokrytí.

0

0.05 0.1 0.15 0.2 0.25 0.3 0.35 pL k

Obrázek 7: Určení nejkratšího intervalu pokrytí veličiny Y reprezentované 10 000 čísly s gamma rozdělením. Nahoře: Histogram. Uprostřed: Kumulativní distribuční funkce. Dole: Závislost délky intervalu pokrytí na „levostrannéÿ pravděpodobnosti. Čárkovaně je vyznačen průměr dat, čerchovaně jsou vyznačeny meze nejkratšího intervalu pokrytí. Hvězdičkou je označen nejkratší interval pokrytí.

92


Obrázek 8: Výpočet příkladu 1 v Excelu. Poslední dva řádky nejsou ukončené středníkem, jinak by GNU Octave nebo Matlab nezobrazily výsledek na obrazovku. Také si můžete vytvořit třeba v poznámkovém bloku soubor nazvaný priklad1.m a zapište do něj příkazy. Soubor uložte do vámi zvoleného adresáře (například c:\). V GNU Octave nebo Matlabu pak pomocí příkazu cd přejděte do zvoleného adresáře (cd c:\) a spusťte soubor příkazem priklad1 (tj. název souboru bez přípony .m). Co znamenají jednotlivé příkazy? Příkaz na řádku č. 1 nastaví hodnoty proměnných. N udává počet vygenerovaných pseudonáhodných čísel, neboli ve svém důsledku počet opakování metody Monte Carlo. pX1 je průměrná hodnota veličiny X1 , uX1 je nejistota veličiny X1 , obdobně pro pX2 a uX2. Řádek č. 2 vygeneruje do proměnné X1 matici o velikosti 1 × N obsahující pseudonáhodná čísla podle normálního rozdělení se střední hodnotou 200 a směrodatnou odchylkou 25. Podobně řádek č. 3 pro proměnnou X2. V řádku č. 4 aplikujeme model měření. Jelikož GNU Octave a Matlab chápou všechny proměnné jako matice (vektory), proměnná Y je matice o velikosti 1×N obsahující prvky yi = x1i + x2i , kde i = 1 . . . N . Histogram, neboli pravděpodobnostní funkce, je vykreslen v řádku č. 5, standardní nejistota je vypočtena v řádku č. 6. Jelikož rozdělení pravděpodobnosti výstupní veličiny je dle pohledu na histogram normální, stačí použít funkci std počítající směrodatnou odchylku.

6

Tedy rozdělení odporů v jednotlivých třídách přesnosti neodpovídají normálnímu rozdělení (viz obr. 9). Mějme dva odpory druhé třídy o nominálních hodnotách 200 Ω a 300 Ω. Celkový odpor bude opět jako v prvním příkladu roven: y = 200 + 300 = 500 .

(9)

Nejistota celkového odporu je ale komplikovanější. Kdybychom pracovali pouze v rámci GUF, mohli bychom zvolit přiblížení a uvažovat, že standardní nejistota odporů druhé třídy je stejná jako směrodatná odchylka všech vyrobených odporů, tedy 25 Ω. Výsledná standardní nejistota výstupní veličiny pak vychází 42,4 Ω. Výsledek je ovšem podhodnocen. Tento příklad není řešitelný metodou GUF, jelikož rozdělení pravděpodobnosti vstupních i výstupních veličin jsou příliš komplikované. Výpočet tohoto příkladu v Excelu pomocí MMC je možný, ale je již poměrně složité vše zadat do buněk. Proto je výpočet příkladu 2 pomocí MMC v Excelu pouze naznačen a celý je vyřešen pouze pro GNU Octave a Matlab. 6.1

Výpočet příkladu 2 v Excelu

Částečně zpracovaný příklad v Excelu lze opět stáhnout na internetových stránkách [9]. Speciální rozdělení pravděpodobnosti vstupních veličin X1 a X2 je provedeno výběrem z normálního rozdělení. Všechna vygenerovaná pseudonáhodná čísla, která jsou v intervalu

Druhý příklad

hXi − 30, Xi + 30i ,

Druhý příklad je podobný prvnímu, tj. mějme dva odpory zapojené do série a zajímá nás celkový odpor. Neznáme ale nejistotu použitých odporů, pouze třídu odporů specifikovanou výrobcem. Výrobce odporů má sériovou výrobní linku. Vyrobené odpory mají normální rozdělení hodnot se směrodatnou odchylkou 25 Ω. Ale výrobce prodává přesnější odpory za vyšší cenu. To dělá tak, že každý vyrobený odpor změří. Pokud se jeho hodnota odchyluje od nominální hodnoty o méně než 30 Ω, zařadí ho do první třídy. Zbývající vyrobené odpory zařadí do druhé třídy.

(10)

jsou zahozena. Tím je získáno rozdělení pravděpodobnosti odporů druhé třídy. To ale znamená, že počet pseudonáhodných čísel výsledného rozdělení pravděpodobnosti je také náhodný. Jelikož práce v Excelu je komplikovaná, příklad je zpracovaný jen částečně. Není zajištěn počet opakování MMC, jelikož počet pseudonáhodných čísel vstupních veličin je náhodný. Také není zpracován algoritmus nalezení nejmenšího intervalu pokrytí pravděpodobnostní funkce

93


0.03

6.2

p

výsledné veličiny. Proto není vypočtena výsledná nejistota výstupní veličiny Y . Příklad je kompletně zpracován jen pro GNU Octave a Matlab, viz následující kapitola.

0 0.03

Výpočet příkladu 2 v GNU Octave a Matlabu

3 4 5 6 7 8 9

p

p

0.010 0.005

p

0.000 0.015

p

0 0.03 Y = X1 + X2 0.02 0.01 0 100 200 300

400 Y (Ω)

500

600

Kdybychom v řádcích č. 2 a 5 vygenerovali pouze N čísel, po aplikaci výběru na řádcích č. 3 a 6 by nám zbylo méně než N čísel. Navíc, jelikož používáme pseudonáhodná čísla, proměnné X1 a X2 by nemusely vždy obsahovat stejný počet čísel. Proto je proveden výběr prvních N čísel na řádcích č. 4 a 7. Tento postup si můžeme dovolit, protože používáme kvalitní generátor pseudonáhodných čísel. V řádku č. 8 se aplikuje model měření. Řádek č. 9 vykreslí histogram, neboli pravděpodobnostní funkci výsledné veličiny Y , viz obr. 10. Z histogramu výstupní veličiny Y je zřejmé, že rozdělení pravděpodobnosti neodpovídá normálnímu rozdělení, a proto nemůžeme vypočítat nejistotu pomocí směrodatné odchylky. Musíme najít nejkratší interval pokrytí, který obsahuje 68,27 % vygenerovaných pseudonáhodných čísel. Tento interval pokrytí pak odpovídá standardní nejistotě veličiny Y pro pravděpodobnost pokrytí 68,27 %. K tomu se hodí funkce scovint, která je v souboru scovint.m. Tuto funkci pro GNU Octave lze stáhnout z internetových stránek [9]. Nejistotu u(y) vypočteme následovně (druhý příkaz není ukončen středníkem):

0.015

0.000 0.015 0.010 0.005 0.000

X2

Obrázek 10: Pravděpodobnostní funkce vstupních veličin X1 a X2 a výsledné veličiny Y z příkladu 2.

Příkazy v řádku č. 1 určují hodnoty proměnných. N udává počet opakování MMC, tedy počet požadovaných pseudonáhodných čísel. Sice neexistuje funkce, která by přímo vygenerovala pseudonáhodná čísla dle požadovaného rozdělení pravděpodobnosti, ale lze postupovat následovně. Nejprve je na řádku č. 2 vygenerováno dostatečné množství (mnohem více než N ) pseudonáhodných čísel s normálním rozdělením, průměrnou hodnotou 200 a směrodatnou odchylkou 25. Poté jsou na řádku č. 3 vybrána pouze čísla, která jsou větší než 30 nebo menší než 30, tedy větší než průměrná hodnota pX1 plus výrobní limit limX a menší než pX1 mínus limX. Jelikož potřebujeme jen N pseudonáhodných čísel, je na řádku č. 4 vybráno prvních N čísel z proměnné X1. Obdobně jsou vygenerovány hodnoty pro proměnnou X2 na řádcích č. 5–7.

0.010 0.005

0.02 0.01

Nrez =1 e6 ; N =1 e4 ; pX1 =200; pX2 =300; uX =25; limX =30; X1 = normrnd ( pX1 , uX ,1 , Nrez ) ; X1 = X1 ( X1 < pX1 - limX | X1 > pX1 + limX ) ; X1 = X1 (1: N ) ; X2 = normrnd ( pX2 , uX ,1 , Nrez ) ; X2 = X2 ( X2 < pX2 - limX | X2 > pX2 + limX ) ; X2 = X2 (1: N ) ; Y = X1 + X2 ; hist (Y ,50) ;

p

2

X1

0.01

Příklad lze vypočítat pomocí následujících příkazů: 1

0.02

[uY,ulY,urY]=scovint(Y,0.6827) uY=uY./2 První parametr funkce scovint jsou hodnoty výstupní veličiny získane z MMC, druhý parametr je pravděpodobnost pokrytí. Funkce vrací tři parametry. První je délka nejkratšího intervalu pokrytí, druhý parametr je levá mez intervalu pokrytí, třetí parametr je pravá mez intervalu pokrytí. Existují výstupní veličiny, jejichž pravděpodobnostní funkce nejsou symetrické a je třeba uvádět nejistotu ve formě levé i pravé meze. V tomto příkladu je výsledkem symetrická pravděpodobnostní funkce, tedy stačí výsledný nejkratší interval podělit dvěma a získáme standardní nejistotu. Výsledná nejistota je 61,3 Ω pro pravděpodobnost

140 160 180 200 220 240 260 Y (Ω) Obrázek 9: Nahoře: Pravděpodobnostní funkce všech vyrobených odporů s nominální hodnotou 200 Ω. Uprostřed: Pravděpodobnostní funkce těchto odporů zařazených do první třídy přesnosti. Dole: Pravděpodobnostní funkce odporů zařazených do druhé třídy přesnosti.

94

VOL.16, NO.2, APRIL 2014 pokrytí 68,27 %, nebo 106 Ω pro pravděpodobnost pokrytí 95,45 %. Je vidět, že na rozdíl od normálního rozdělení druhé číslo není dvojnásobkem prvního. Nejistota vypočtená metodou GUF je tedy podhodnocena o 31 %. Pokud by se nepočítal nejkratší interval pokrytí, ale jen směrodatná odchylka, došlo by pro případ pravděpodobnosti pokrytí 95,45 % k nadhodnocení o 14 %.

7

Závěr

Z příkladů je vidět, že MMC je poměrně jednoduchý ale velmi účinný nástroj k výpočtu i komplikovaných nejistot. Jednoduchost použití je ale závislé na použitém programovém vybavení. Jak bylo ukázáno na druhém příkladě, komplikovanější zadání je obtížně řešitelné v běžných kancelářských programech, ale při použití specializovaných matematických programů lze snadno získat správný výsledek, a to i bez počítání parciálních derivací a stupňů volnosti. Doufám, že se mi podařilo sestavit co nejjednodušší návod jak začít s metodou Monte Carlo pro výpočet nejistot. Součástí článku nejsou další aspekty MMC, jako jsou např. korelované vstupní veličiny a určování nejistot pro více výstupních veličin. Také nebyly řešeny všechny předpoklady použití MMC.

Literatura [1]

JCGM. Evaluation of measurement data - Guide to the expression of uncertainty in measurement. 1995. ISBN 92-67-10188-9.

[2]

JCGM. Evaluation of measurement data - Supplement 1 to the “Guide to the expression of uncertainty in measurement” - Propagation of distributions using a Monte Carlo method. 2008.

[3]

JCGM. Evaluation of measurement data – Supplement 2 to the “Guide to the expression of uncertainty in measurement” – Extension to any number of output quantities. 2011.

[4]

COX, M. G. Propagation of distributions by a Monte Carlo method, with an application to ratio models. The European Physical Journal Special Topics. 2009, roč. 172, č. 1, s. 153–162. Dostupné také z: hhttp://www. springerlink . com /index / 10 . 1140 / epjst/ e2009 01048-0i. ISSN 1951-6355.

[5]

ROTZ, Wendy; FALK, Eric; JOSHEE, Archana. A Comparison of Random Number Generators Used in Business–2004 Update. In. Proceedings of the Survey Research Methods Section of the American Statistical Association. Č. 2. 2004, s. 1316–1319.

[6]

Matlab. [online]. [cit. 2014-01-20]. Dostupné z: hhttp: //www.mathworks.com/products/matlab/i.

[7]

GNU Octave. [online]. [cit. 2014-01-20]. Dostupné z: hhttp://www.gnu.org/software/octave/i.

[8] [9]

95

Octave-Forge. [online]. [cit. 2014-01-20]. Dostupné z: hhttp://octave.sourceforge.neti. Instalační soubor GNU Octave a ukázkové příklady a skripty. [online]. [cit. 2014-01-20]. Dostupné z: hhttp: //kaero.wz.cz/mmcjednoduse.htmli.

Jak na nejistoty metodou Monte Carlo jednoduše a bez drahých programů

Recommend Documents