📊 Metody Monte Carlo: státnice

1. Co je vlastně metoda Monte Carlo

Monte Carlo (MC) je výpočetní metoda, která modeluje náhodné děje pomocí počítače a výsledky následně statisticky zpracovává, podobně jako fyzik vyhodnocuje měření v laboratoři. Hodí se ve dvou typech úloh:

• Skutečně náhodné procesy: kde stochasticita je v zadání: rozpad radioaktivních jader, čekací doby v bance, krach trhu, transport částic skrz stínění reaktoru. Tam si MC „povídá“ s realitou v její přirozené řeči.
• Deterministické úlohy: kde si náhodu přivedeme jako nástroj. Klasický učebnicový příklad je výpočet integrálu nebo součtu řady. Hodnota je předem daná, ale spočítat ji analyticky je obtížné, takže ji „naměříme“ simulací.

Princip vzdáleně připomíná zjišťování poměru obvodu kruhu k jeho průměru tak, že do čtverce s vepsaným kruhem hážeme náhodné body a sledujeme, kolik procent z nich padne dovnitř. To je nejjednodušší ukázka MC: kvadratura plochy. Stejnou logikou se ovšem počítají věci nesrovnatelně složitější (vícerozměrné integrály, optimalizace stochastických procesů, finanční modely, simulace neutronového toku).

2. Proč právě rovnoměrné rozdělení $R(0,1)$

$R(0,1)$ je v MC svatý grál. Veškerá další rozdělení (diskrétní, exponenciální, normální, Weibullovo, libovolně exotická) se z něj odvozují transformací. Stačí tedy mít kvalitní zdroj rovnoměrných čísel a o zbytek se postará matematika.

Vlastnosti $Y \sim R(0,1)$, které využíváme:

• Střední hodnota $EY = \int_0^1 x\,dx = \tfrac{1}{2}$.
• Rozptyl $DY = \int_0^1 (x - \tfrac{1}{2})^2\,dx = \tfrac{1}{12}$.
• Distribuční funkce $F(x) = x$ pro $x \in (0,1)$. Z toho plyne klíčový fakt: pokud máme libovolnou distribuční funkci $G$, pak $G^{-1}(Y)$ má rozdělení popsané $G$. To je princip inverzní transformace, kterou rozebíráme v Otázce 3.

3. Kde rovnoměrná čísla brát

Tabulky náhodných čísel

Historicky se sestavovaly tabulky tisíců náhodných číslic (například RAND Corporation vydala v roce 1955 knihu „A Million Random Digits with 100,000 Normal Deviates“). Dnes mají význam spíš muzejní, ale ukazují, že problém zdroje náhodnosti je hodně starý.

Hardwarové (fyzikální) generátory

Vychází z fyzikálního zdroje šumu: tepelný šum rezistoru, šum lavinové diody, časy detekce kvantové fluktuace, lavinová ionizace u Geigerova čítače. Vzorkováním šumu dostáváme „pravé“ náhodné bity.

Výhoda: teoreticky perfektní náhoda, nelze ji předpovědět.
Nevýhoda: výstup nelze reprodukovat. To je pro vědecké výpočty problém. Pokud najdeme v simulaci anomálii, chceme ji umět zopakovat a debugovat. S hardwarovým generátorem to nejde, protože tatáž posloupnost se už nevrátí.

Softwarové (pseudonáhodné) generátory

Algoritmus, který z počátečního stavu (seedu) deterministicky generuje posloupnost čísel, jež se chová „jako náhodná“. Stejný seed dá vždy stejnou posloupnost, takže simulace je reprodukovatelná. Toto je dnes naprostý standard, i ten random.random() v Pythonu je softwarový generátor.

Pojem „pseudo“ neříká, že je to horší, jen že posloupnost je předem daná. Důležité je, aby žádný rozumný statistický test nedokázal posloupnost rozlišit od skutečně náhodné.

4. Co musí dobrý PRNG splňovat (formální požadavky)

Označme $\theta = \{x_i\}_{i=0}^{\infty}$ pseudonáhodnou posloupnost a $X^{(n)} = (x_1, \dots, x_n)$ libovolný výběr o $n$ prvcích. Aby $\theta$ obstála jako model výběru z $Y \sim R(0,1)$, musí pro každý výběr platit:

1. Konvergence střední hodnoty:
   $\displaystyle \lim_{n\to\infty}\frac{1}{n}\sum_{i=1}^{n} x_i = \frac{1}{2}$
   Logika: $EY = \tfrac{1}{2}$, aritmetický průměr musí dlouhodobě konvergovat k téhle hodnotě. Pokud by průměr šel třeba k 0,55, generátor systematicky preferuje vyšší čísla.
2. Konvergence rozptylu:
   $\displaystyle \lim_{n\to\infty}\frac{1}{n}\sum_{i=1}^{n}(x_i - \bar{x})^2 = \frac{1}{12}$
   Pokud by rozptyl byl menší, čísla by se kupila kolem 0,5; pokud by byl větší, hustila by se u krajů.
3. Rovnoměrnost ve smyslu objemu: Pro libovolnou měřitelnou podmnožinu $I_v \subset (0,1)$ s objemem $V(I_v) = v$ musí počet bodů, které do ní padnou, splňovat
   $\displaystyle \lim_{n\to\infty}\frac{N_v(X^{(n)})}{n} = v$
   Lidsky: do 10% výseku interválu musí padnout 10 % bodů. To je test, který frekvenční histogramy.
4. $k$-rovnoměrnost: Vezmeme $k$-tice po sobě jdoucích čísel a interpretujeme je jako body v $k$-rozměrné jednotkové krychli $C_k = \langle 0,1 \rangle^k$. Pro libovolnou měřitelnou podmnožinu $G \subset C_k$ musí
   $\displaystyle \lim_{n\to\infty}\frac{N_G}{n} = V(G)$
   Tohle je kritický bod, který odhalí slabší generátory. Posloupnost může vypadat v 1D rovnoměrně, ale jakmile $k$-tice tvoří v $k$-D viditelné mřížky či shluky, generátor padá. Klasický průšvih RANDU (viz níže) se odhalil právě tímto.
5. Nulová autokorelace se zpožděním $t$: Pro libovolné $t > 0$ musí platit
   $\displaystyle \lim_{n\to\infty}\frac{1}{n}\sum_{i=1}^{n}\left(x_{i-t} - \tfrac{1}{2}\right)\left(x_i - \tfrac{1}{2}\right) = 0$
   Lidsky: hodnota $x_i$ nesmí být statisticky predikovatelná z $x_{i-t}$. U deterministických PRNG je $x_i$ samozřejmě plně určeno celou historií, ale statisticky to nesmí být vidět.

5. Lineární kongruenční generátor (LCG)

Vůbec nejjednodušší a historicky první rozšířený PRNG. Pochází z přelomu 40. a 50. let (Derrick Lehmer). Stačí mu jedno násobení, jedno sčítání a jedno modulo:

Význam parametrů:

• $x_0$ je seed (semínko). Počáteční hodnota, která určuje celou budoucí posloupnost. Stejný seed = stejná posloupnost.
• $K$ je multiplikátor. Klíčový parametr, špatně zvolený rozhodne o krátké periodě nebo silné korelaci.
• $B$ je inkrement. Pokud $B = 0$, mluvíme o čistě multiplikativním kongruenčním generátoru.
• $M$ je modul. Často $M = 2^n$ kvůli rychlosti (modulo přes mocninu dvojky je jen ořezání bitů).

Perioda je délka cyklu, po kterém se posloupnost začne opakovat. U LCG je nejvýš $M$. Existují Hull–Dobellovy podmínky, které zaručují, že perioda dosahuje maxima (vztahy mezi $K$, $B$, $M$ a jejich společnými děliteli).

6. Mersenne Twister (MT19937)

Dnes nejpoužívanější PRNG. Najdeš ho jako výchozí v Pythonu (random), Ru, MATLABu, PHP. Navrhli ho Matsumoto a Nishimura v roce 1997. Stojí na principu lineárního posuvného registru se zpětnou vazbou, do něhož je vložen geniální trik s maskováním bitů.

Parametry varianty MT19937 (číslo 19937 souvisí s periodou, viz dále):

• $w = 32$ bitů na slovo.
• $n = 624$ slov ve vnitřním stavu (registru).
• $m = 397$ posun pro zpětnou vazbu.
• $r = 31$ hranice maskování.
• $A$ konstantní matice $w \times w$ s prvky z $\{0,1\}$, ale ve skutečnosti realizovaná jako rychlá bitová operace (shift + XOR).

Rekurentní vzorec, který generuje nové slovo:

• $x_k^{\{u(w-r)\}}$ horních $w - r$ bitů slova $x_k$.
• $x_{k+1}^{\{l(r)\}}$ dolních $r$ bitů slova $x_{k+1}$.
• $|$ označuje zřetězení bitů (spojení).
• $\oplus$ je bitový XOR (sčítání modulo 2).

Proč je MT tak dobrý:

• Perioda je rovna $2^{19937} - 1$, což je Mersennovo prvočíslo (odsud název). Číslo má skoro 6000 desetinných cifer. Žádný reálný výpočet ho neumí vyčerpat ani za miliardy let.
• $k$-rovnoměrnost až do $k = 623$ rozměrů. Nadroviny, na kterých $k$-tice leží, jsou tak husté, že se v praxi neprojeví.
• Rychlost je výborná. Algoritmus používá jen bitové operace, žádné násobení velkých čísel.

7. Transformace na obecný rovnoměrný interval $R(a,b)$

Často potřebujeme náhodu z jiného intervalu než $(0,1)$. Vzorec je triviální:

Roztáhneme interval $(0,1)$ na šířku $b - a$ a posuneme začátek na $a$.

8. Předjímka inverzní metody

Pro úplnost obecný princip, který využijeme v Otázce 3. Pokud máme hustotu $f(x)$ a distribuční funkci $F(x) = \int_{-\infty}^{x} f(t)\,dt$, a chceme generovat realizaci náhodné veličiny s rozdělením $F$, položíme $F(\xi) = \Upsilon$ a vyřešíme $\xi = F^{-1}(\Upsilon)$.

9. Golombovy požadavky na binární posloupnosti

Když se generátor implementuje jako binární posloupnost (bit po bitu), používáme se zpřísňující kritéria od Solomona Golomba. Jsou to nutné, ale ne dostačující podmínky kvality:

1. Vyváženost. V jedné periodě se počet jedniček a nul liší nejvýš o jedna. Polovina bitů musí být 1, polovina 0.
2. Distribuce běhů (runs). Běh je nepřerušená řada stejných bitů. V náhodné posloupnosti by polovina běhů měla mít délku 1, čtvrtina délku 2, osmina délku 3 atd. Geometrické rozdělení.
3. Ideální autokorelační funkce. Když posloupnost porovnáme samu se sebou posunutou o $t$, počet shod se má lišit od počtu neshod o konstantu nezávislou na $t$ (kromě případu $t = 0$, kde se vše shoduje).

Co testy vlastně dělají a co neumějí

Pokud generátor selže byť v jediném testu, statisticky se nechová jako výběr z $R(0,1)$ a nesmíme ho použít. Pokud projde všemi testy, nelze tvrdit, že je správný. Statistické testy jsou nutné, ne dostačující podmínky kvality. Nový test může odhalit nedostatek, který předchozí přehlédly. Klasický příklad: RANDU prošel jednorozměrnými testy bez problému a teprve trojrozměrný test odhalil katastrofu (15 nadrovin).

Testy provádíme jako statistické testy hypotéz. Pracovní hypotéza $H_0$ zní: „posloupnost se chová jako výběr z $R(0,1)$“. Spočítáme testovou statistiku a porovnáme ji s kritickou hodnotou, která odpovídá zvolené hladině významnosti, typicky $\alpha = 0{,}05$ (95% spolehlivost). Pokud statistika kritickou hodnotu překročí, hypotézu zamítáme: generátor je podezřelý.

1. Frekvenční test (test dvojic 0/1)

Nejjednodušší test. Ověřuje, že v binárním zápisu posloupnosti se 0 a 1 vyskytují přibližně stejně často. U skutečně rovnoměrného generátoru je pravděpodobnost 0 i 1 na libovolném bitu rovna $1/2$, takže v posloupnosti délky $n$ bitů by se mělo objevit zhruba $n/2$ jedniček.

Označme $n_0$, $n_1$ počty nul a jedniček ($n = n_0 + n_1$). Testovou statistiku definujeme jako:

• Rozdělení statistiky: pro dostatečně velký vzorek ($n > 10$) se $X_1$ řídí $\chi^2$ rozdělením s jedním stupněm volnosti. Stupeň volnosti je 1, protože ze dvou kategorií (0 a 1) je jen jedna „volná“: jakmile známe $n_0$, je $n_1 = n - n_0$ určeno.
• Kritická hodnota: $X_1 < 3{,}84$ s 95% spolehlivostí. Hodnota 3,84 je 95% kvantil $\chi^2_1$, tedy plocha pod hustotou $\chi^2_1$ od 0 do 3,84 činí přesně 0,95.

Pohled přes binomické rozdělení

Pokud máme $m$ slov délky $n$ bitů, je celkový počet jedniček náhodná veličina s binomickým rozdělením s parametry $mn$ a $p = 1/2$. Pro velká $m$ se binomické rozdělení aproximuje normálním a počet jedniček s pravděpodobností 0,95 leží v intervalu:

Tento intervalový pohled je ekvivalentní statistice $X_1$, jen vyjádřený jinou formou.

2. Test sérií (rovnoměrnost dvojic)

Frekvenční test snadno propustí patologie jako 0101010101...: tato posloupnost má sice 50 % jedniček, ale je striktně periodická. Test sérií jde o úroveň výš a kontroluje, jestli i dvojice po sobě jdoucích bitů (00, 01, 10, 11) jsou rovnoměrně rozloženy.

Označme $n_0$, $n_1$ jako dříve a $n_{00}, n_{01}, n_{10}, n_{11}$ počty výskytů jednotlivých dvojic. Statistika:

• Rozdělení: $\chi^2$ se dvěma stupni volnosti. Máme čtyři typy dvojic, ale jsou vázány na celkový počet a na výsledky frekvenčního testu, takže zbývají 2 volné parametry.
• Kritická hodnota: $X_2 < 5{,}99$. Opět 95% kvantil $\chi^2_2$.

Druhá varianta (intervalová)

Místo bitů můžeme test sérií provádět nad reálnými čísly: rozdělíme $(0,1)$ na $r$ stejně velkých nepřekrývajících se podintervalů $I_1, \dots, I_r$ a sledujeme, kolikrát padly vygenerované hodnoty do každého z nich. Pokud bychom navíc sledovali dvojice po sobě jdoucích čísel, dostaneme se k testování dvourozměrné rovnoměrnosti (které je úvod ke $k$-rovnoměrnosti).

3. Pokerový test

Generalizace testu sérií: místo dvojic sleduje $m$-tice po sobě jdoucích bitů. Jméno pochází z analogie s kartami v pokeru, kde také sledujeme kombinace „karet v ruce“ určité délky.

Postup:

1. Posloupnost $X$ převedeme do binárního tvaru a získáme řetězec $s$ o délce $N$ bitů.
2. Zvolíme délku úseku $m$ a rozdělíme $s$ na $k = \lfloor N/m \rfloor$ nepřekrývajících se úseků.
3. Existuje $2^m$ možných kombinací úseků délky $m$. Označme $n_i$ počet skutečných výskytů $i$-tého typu úseku.
4. Pro statistickou průkaznost musí platit $k \ge 5 \cdot 2^m$, aby na každý možný úsek připadalo v průměru aspoň 5 výskytů.

Rozdělení: $\chi^2$ s $2^m - 1$ stupni volnosti. Jeden stupeň „ujídá“ vazba na celkový počet úseků $k$.

4. Test běhů (Runs test)

Sleduje strukturu nepřerušených sekvencí stejných bitů, tzv. běhů. Test ověřuje, zda jejich distribuce odpovídá Golombovým požadavkům (které jsme zavedli v Otázce 1).

Pojmy:

• Blok nepřerušená řada jedniček. Pokud má délku $i$, počítáme ho do $B_i$.
• Mezera nepřerušená řada nul délky $i$, počítáme do $M_i$.
• Běh souhrnný název pro blok nebo mezeru.

Teoretický očekávaný počet běhů délky $i$ v posloupnosti délky $n$ bitů:

Pro statistickou průkaznost nás zajímají jen běhy, jejichž očekávaný počet je $\ge 5$. Najdeme tedy největší $k$ tak, aby $e_k \ge 5$. Delší běhy do testu nezahrnujeme.

Testová statistika:

Rozdělení: $\chi^2$ s $2k - 2$ stupni volnosti.

5. Autokorelační test

Hledá skryté závislosti mezi bity, které jsou od sebe vzdálené o $d$ pozic. U opravdové náhody musí být znalost bitu $s_i$ úplně nesouvisející s bitem $s_{i+d}$ pro každé $d > 0$.

Postup:

1. Zvolíme zpoždění $d$ tak, aby $1 \le d \le \lfloor n/2 \rfloor$. $d$ nesmí být rovno periodě generátoru (to by trivializovalo test).
2. Spočítáme počet párů, kde se bity liší:
   $A(d) = \displaystyle\sum_{i=1}^{n - d - 1} s_i \oplus s_{i+d}$
   $\oplus$ je XOR, který vrátí 1 právě tehdy, jsou-li bity různé.
3. U skutečně náhodné posloupnosti se polovina dvojic shoduje a polovina liší, takže očekávaná hodnota $A(d)$ je $(n-d)/2$. Odchylku měří statistika:
   $X_5 = 2\,\dfrac{A(d) - \frac{n - d}{2}}{\sqrt{n - d}}$

Rozdělení: $X_5 \sim N(0, 1)$ pro $n - d > 10$. Generátor projde, pokud $|X_5| \le 1{,}96$ (kritická hodnota pro 95% oboustranný interval normálního rozdělení).

6. Maurerův univerzální statistický test

Nejsilnější ze standardních testů. Místo aby hledal konkrétní vzor (rovnoměrnost nul/jedniček, dlouhé běhy), zkoumá celkovou entropii posloupnosti přes informační teorii.

Základní myšlenka:

• Skutečně náhodná posloupnost má maximální entropii a je nezkomprimovatelná bez ztráty informace.
• Pokud by v datech existoval jakýkoli skrytý řád (i velmi jemný), dobrý kompresní algoritmus by ho odhalil a soubor zmenšil.
• Maurer formálně používá vzdálenosti mezi opakováními bitových vzorů. Pokud se vzory opakují příliš pravidelně, posloupnost je nedostatečně náhodná.

Proč je univerzální a silný:

• Postihuje široké spektrum nedostatků, které jednodušší testy minou.
• Velmi efektivně odhaluje skryté dlouhé periodicity.
• Pokud generátor projde Maurerovým testem, je to silná indikace, že vypadá jako $R(0,1)$.

Vlastnosti generátorů a jejich porovnání

Lineární kongruenční generátor (LCG)

• Snadná implementace. Jedno násobení, jedno sčítání, jedno modulo. Pokud $M = 2^n$, je modulo zdarma (ořezání bitů). Výpočetně nenáročný.
• Krátká perioda. Nejvýše $M$, prakticky často méně. Při $M = 2^{31}$ je perioda pod 2 miliardy, což moderní simulace lehce vyčerpá.
• Korelace hodnot. Po sobě jdoucí čísla nejsou dostatečně nezávislá, generátor selhává v testu $k$-rovnoměrnosti.
• Nadroviny. Geometrický důsledek: pokud $k$-tice po sobě jdoucích čísel chápeme jako body v $k$-rozměrném prostoru, body neleží rovnoměrně, ale shlukují se na $(k{-}1)$-rozměrných nadrovinách. U RANDU jen 15 nadrovin v 3D, což je destruktivní.

Mersenne Twister (MT19937)

• Parametry standardní verze: $w = 32$, $n = 624$, $m = 397$, $r = 31$.
• Obrovská perioda $2^{19937} - 1$ (Mersennovo prvočíslo, odsud název).
• $k$-rovnoměrnost až do $k = 623$ rozměrů.
• Rychlý. Bitové posuny a XOR, žádné drahé operace.
• Není kryptograficky bezpečný. Z 624 po sobě jdoucích výstupů lze obnovit vnitřní stav a od té chvíle generátor přesně predikovat. To pro vědecké simulace nevadí, pro generování klíčů ano.

Tabulka srovnání

Princip: $R(0,1)$ jako univerzální stavební kámen

Mám-li kvalitní generátor čísel z $R(0,1)$, dokážu z něj sestrojit realizace prakticky libovolného rozdělení. To je krásná vlastnost a důvod, proč MC vystačí s jedním základním PRNG.

Rozdělíme úvahu na dvě části, protože diskrétní a spojité veličiny vyžadují jiné techniky:

• Diskrétní veličiny: hodnoty $x_k$ s pravděpodobnostmi $p_k$. Interval $(0,1)$ se rozkrájí na „škatulky“ délky $p_k$ a sleduje se, do které $\Upsilon$ padne.
• Spojité veličiny: hustota $f(x)$, distribuční funkce $F(x) = \int_{-\infty}^{x} f(t)\,dt$. Používáme inverzní transformaci, zamítací metodu, nebo speciální triky pro vybraná rozdělení (Box–Muller pro normální).

A) Diskrétní náhodné veličiny

Chceme vyrobit realizaci náhodné veličiny $\xi$ popsané rozdělením pravděpodobnosti:

kde $x_k$ jsou možné hodnoty a $I$ je množina indexů (často $I = \{0, 1, \dots, n\}$ nebo $I = \mathbb{N}$).

Základní algoritmus

Geometrická idea: úsečku $(0,1)$ rozdělíme na podintervaly o délkách $p_0, p_1, p_2, \dots$ Vygenerujeme $\Upsilon \sim R(0,1)$ a zjistíme, do kterého podintervalu padl. Implementace nevyžaduje předpočítané hranice, ale postupně odečítá $p_k$ od pomocné proměnné $S$:

Složitost algoritmu

Spočítejme střední dobu výpočtu $\tau$. První krok proběhne vždy a stojí $\tau_R$ (čas na vygenerování $\Upsilon$). Kroky 2 a 3 se opakují podle toho, jak daleko je hledaný index $k$. Pravděpodobnost, že výpočet skončí přesně po $k+1$ iteracích, je $p_k$.

Označme $\tau_p$ čas na získání hodnoty $p_k$ a $\tau_+$ čas jedné odečítací operace. Pak:

Protože $\sum_{k} k\,p_k = E\xi$ (definice střední hodnoty diskrétní veličiny), můžeme psát:

Závěr je důležitý: střední doba výpočtu roste lineárně se střední hodnotou indexu $E\xi$. Pokud máme rozdělení s velkou $E\xi$ (např. Poissonovo s $\lambda = 1000$), algoritmus v průměru projde tisícem iterací než najde výsledek. To je důvod pro optimalizace, které následují.

Modifikovaný algoritmus 1: setřídění sestupně

Pokud místo přirozeného pořadí indexů použijeme pořadí, kde největší pravděpodobnosti jsou na začátku, sníží se $E\xi$ a tím i průměrný počet iterací. Stačí jednou setřídit posloupnost $p_k$ sestupně:

Algoritmus je pak identický, jen pracuje s přeindexovanou posloupností.

Modifikovaný algoritmus 2: start od vrcholu (pro „kopcovitá“ rozdělení)

U některých rozdělení (typicky Poissonovo s velkým $\lambda$) jsou pravděpodobnosti $p_k$ nejprve malé, kolem vrcholu $h$ velké a pak zase klesají. Setřídění by zničilo přirozený index. Místo toho startujeme rovnou u vrcholu $h$ a podle $\Upsilon$ se rozhodneme, zda jdeme doleva nebo doprava.

Předpočítáme součet pravděpodobností od počátku po vrchol:

Tři konkrétní diskrétní rozdělení

1. Bernoulliho rozdělení s parametrem $p$

Pokus, který skončí buď úspěchem (hodnota 1) s pravděpodobností $p$, nebo neúspěchem (0) s pravděpodobností $1 - p$. Nejjednodušší možný diskrétní model: pojistka, hod mincí s nestejnoměrnou váhou, klinický test.

2. Rovnoměrné rozdělení na $\{1, 2, \dots, n\}$

Spravedlivá kostka, losování. Každá hodnota má stejnou pravděpodobnost $1/n$.

3. Poissonovo rozdělení s parametrem $\lambda$

Modeluje počet událostí v jednotce času (počet zákazníků za hodinu, počet alfa-rozpadů za sekundu). Má hustotu $P(\xi = k) = \dfrac{\lambda^k}{k!}e^{-\lambda}$. Naivně bychom mohli použít základní algoritmus, ale existuje elegantní trik založený na vztahu mezi Poissonovým a exponenciálním rozdělením: čekací doby mezi událostmi jsou exponenciální. Místo sčítání času se ekvivalentně násobí náhodná čísla.

Proč to funguje: čekací doby mezi Poissonovými událostmi mají exponenciální rozdělení s parametrem $\lambda$, generujeme je inverzní transformací jako $\eta_i = -\tfrac{1}{\lambda}\ln \Upsilon_i$. Suma $\sum \eta_i$ překročí jednotku $t = 1$ právě tehdy, když součin $\prod \Upsilon_i$ klesne pod $e^{-\lambda}$. Z toho plyne podmínka „násob, dokud $S > q$“.

B) Spojité náhodné veličiny

1. Inverzní transformace (základní metoda)

Předpokládejme, že chceme generovat realizaci spojité náhodné veličiny s distribuční funkcí $F(x)$. Klíčový fakt: pokud $\Upsilon \sim R(0,1)$, pak $F^{-1}(\Upsilon)$ má rozdělení popsané $F$. To je matematická hodnota inverzní metody.

Postup:

1. Generuj $\Upsilon \sim R(0,1)$.
2. Najdi řešení $\xi$ rovnice $F(\xi) = \Upsilon$, tedy $\xi = F^{-1}(\Upsilon)$.
3. $\xi$ je realizace náhodné veličiny s distribuční funkcí $F$.

Inverzní metoda funguje vždy, když umíme $F^{-1}$ vypočítat. Pro některá rozdělení existuje hezký uzavřený tvar, pro jiná nikoli.

2. Zamítací metoda (Rejection method)

Pro mnoho rozdělení je $F$ buď neuzavíratelná v elementárních funkcích (např. normální rozdělení), nebo příliš komplikovaná pro snadné invertování. V takovém případě používáme zamítací metodu, která pracuje přímo s hustotou $f$.

Předpoklady:

• Hustota $f(x)$ má nosič v omezeném intervalu $\langle a, b \rangle$ (nebo se na něj přibližně omezí).
• $f(x)$ je shora omezená konstantou $C$ (tedy $f(x) \le C$ pro všechna $x \in \langle a, b \rangle$).

Efektivita zamítací metody se měří jako poměr plochy pod křivkou k ploše obdélníka, tedy $\dfrac{1}{C(b-a)}$ (protože pod hustotou je vždy plocha 1). Pokud máme špatně zvolený $C$, generujeme zbytečně mnoho odmítnutých vzorků.

3. Normální rozdělení $N(0,1)$

Normální (Gaussovo) rozdělení je nejdůležitější rozdělení v MC i ve statistice obecně. Distribuční funkce není v uzavřeném tvaru (vyjadřuje se přes nestandardní funkci $\Phi$), takže inverzní metoda v přímé formě selhává. Existují ale dva oblíbené postupy.

Přechod na $N(\mu, \sigma^2)$

Máme-li $\xi \sim N(0,1)$, můžeme snadno škálovat na libovolné normální rozdělení transformací $\eta = \mu + \sigma \xi$. To má rozdělení $N(\mu, \sigma^2)$.

Kontext: kde se simulace transportu používá

Monte Carlo simulace transportu částic je standardní nástroj jaderné fyziky a radiační ochrany. Modeluje průchod toku elementárních částic (fotony, elektrony, neutrony, protony, ionty) materiálním prostředím. Typické úlohy:

• Výpočet stínění: kolik vrstev olova, betonu nebo polyethylenu potřebujeme, aby záření z reaktoru nebo izotopového zdroje bezpečně neproniklo ven? Jak silnou ochranu musí mít lékař při radiologickém zákroku?
• Optimalizace experimentů a jaderných zařízení: kde nejlépe umístit detektor, jakým směrem orientovat svazek, jak vybrat aktivní materiál palivového článku.
• Interakce záření s hmotou v dozimetrii: kolik dávky absorbuje tkáň, jak proběhne radioterapie nádoru, jaké jsou efekty kosmického záření na elektroniku družic.
• Simulace detektorů: jakou účinnost má scintilátor, jaké signály očekávat při různých částicových toků.

V této otázce uvažujeme homogenní tok, tedy částice stejného druhu. Trajektorie jedné částice je posloupností náhodných událostí: emise → přímočarý let → interakce s atomem → buď zánik, nebo pokračování (případně další větve).

Zjednodušující předpoklady modelu

• Mezi migrujícími částicemi nedochází k interakcím. Tento předpoklad platí, pokud hustota toku není extrémní (např. v reaktoru nebo ve stínění). U vysokých toků by se vzájemné interakce nezanedbatelně projevily.
• Pohyb je plně popsán interakcemi s atomy prostředí. Potřebujeme znát pravděpodobnosti jevů (pohlcení, rozptyl, štěpení) a hustotu atomů.
• Mezi interakcemi částice letí rovnoměrně přímočaře. Zakřivení trajektorií gravitací nebo magnetickými poli zanedbáváme (kromě speciálních úloh).
• Sekundární částice vznikají přesně v bodě interakce. Žádné prodlení, žádné rozprostření zdroje.
• Interakce jsou dobře definované jevy s pravděpodobnostmi známými experimentálně nebo z teorie.

Účinné průřezy: měřítko pravděpodobnosti interakce

Účinný průřez (cross section) $\sigma$ je centrální pojem. Geometricky si představ atomy jako kuličky a procházející částici jako kámen. Plocha „terče“, který kámen zasáhne s nějakou pravděpodobností, je účinný průřez. Jednotkou je obvykle barn ($1\,\text{b} = 10^{-28}\,\text{m}^2$).

Diferenciální účinný průřez

Popisuje pravděpodobnost, že částice v daném bodě prostoru se po interakci odchýlí do konkrétního prostorového úhlu.

Význam parametrů:

• $\mathbf{r}$: bod v prostoru, kde k interakci dochází.
• $\boldsymbol{\omega}'$: jednotkový směrový vektor částice před interakcí.
• $\boldsymbol{\omega}$: jednotkový směrový vektor po interakci.
• $E$ (nebo $W$): energie částice.
• $d\Omega = \sin\vartheta\,d\vartheta\,d\varphi$: element prostorového úhlu na povrchu jednotkové koule.
• $\sigma_a$: mikroskopický účinný průřez. Je to podmíněná pravděpodobnost interakce vztažená na jediný atomový terč při jednotkové hustotě dopadajícího toku.

Makroskopický účinný průřez

Mikroskopický průřez popisuje jeden atom. V reálu má prostředí mnoho atomů na jednotku objemu, takže makroskopický průřez je:

kde $N$ je počet terčů na jednotku objemu. Makroskopický průřez vyjadřuje celkovou pravděpodobnost interakce v daném materiálu, jeho rozměr je $\text{m}^{-1}$ (převrácená délka). Intuitivně: $\sigma^{(m)}$ určuje, jak často částice naráží do něčeho.

Volná dráha a optická vzdálenost

Vzdálenost mezi dvěma po sobě jdoucími interakcemi nazýváme volná dráha $\ell$. Je to náhodná veličina s exponenciálním rozdělením. V homogenním prostředí (kde $\sigma^{(m)}$ nezávisí na $\mathbf{r}$) je hustota pravděpodobnosti:

V nehomogenním prostředí, kde se $\sigma^{(m)}(\mathbf{r})$ mění podél dráhy, dostáváme obecnější vztah:

Distribuční funkce:

kde $\tau(t)$ je optická vzdálenost (nebo koeficient pohlcení):

Optická vzdálenost je tedy „integrál tuhosti prostředí“ po dráze. Pravděpodobnost, že částice přežije bez interakce dráhu délky $t^*$, je:

Rozdělení podle typu interakce

Když částice dorazí do bodu interakce, musíme určit, co se s ní stane. Klíčový parametr je střední hodnota počtu vylétávajících (sekundárních) částic, kterou označujeme $\nu$.

Pokud při interakcích vznikají identické sekundární částice (např. neutrony při štěpení), můžeme všechny dílčí druhy interakcí $i = 1, \dots, m$ sloučit do jedné „sumární“ interakce, jejíž parametry jsou váženými průměry:

Algoritmus pro $\nu \le 1$ (pohlcení a rozptyl, bez štěpení)

Sledujeme jednu částici od emise do zániku (pohlcení nebo opuštění oblasti).

Algoritmus pro $\nu > 1$ (štěpení / multiplikace)

Pokud v kroku 4 nastane interakce s $\nu > 1$, situace je složitější: počítač neumí sledovat víc trajektorií najednou, ale částice se větví na strom potomků. Musíme algoritmus upravit a používat zásobník (kontejner), kam si odložíme nezpracované větve.

Počet sekundárních částic $n$ se losuje stochasticky tak, aby jeho střední hodnota byla přesně $\nu$:

• $n = \lfloor \nu \rfloor$ s pravděpodobností $1 - \{\nu\}$.
• $n = \lfloor \nu \rfloor + 1$ s pravděpodobností $\{\nu\}$.

kde $\{\nu\}$ značí desetinnou část a $\lfloor \nu \rfloor$ celou část. Pro $\nu = 2{,}5$ tedy s pravděpodobností 0,5 vzniknou 2 částice a s 0,5 vzniknou 3 částice. Střední hodnota je $0{,}5 \cdot 2 + 0{,}5 \cdot 3 = 2{,}5 = \nu$.

Postup po štěpení:

• Spočítáme směrové vektory pro všech $n$ nově vzniklých částic (každá má svůj náhodný směr odvozený z diferenciálního průřezu).
• První částici simulujeme dál běžným způsobem (jako kdyby šlo o rozptyl).
• Ostatních $n - 1$ částic uložíme do kontejneru (zásobníku).
• Až aktuálně sledovaná větev zanikne (pohlcení nebo opuštění oblasti), počítač sáhne do zásobníku a vyzvedne další odloženou částici k simulaci.
• Výpočet pro danou primární částici končí ve chvíli, kdy je zásobník úplně prázdný.

Modelování v nehomogenním prostředí

V reálu se částice málokdy pohybuje pouze jedním materiálem. Typický scénář: foton vyletí z palivového proutku (uran), prochází vrstvou olova, pak betonem reaktorové budovy, pak vzduchem ven. V každém materiálu má jiný makroskopický průřez.

Tento problém se řeší modelem po částech homogenního prostředí. Prostor rozdělíme na oblasti $V_1, V_2, \dots$, kde v každé je $\sigma^{(m)}$ konstantní. Pokud částice na své dráze prochází postupně oblastmi $V_{i_1}, V_{i_2}, \dots, V_{i_s}$, optická vzdálenost se počítá jako suma příspěvků z jednotlivých bloků:

kde $\lambda_{ij} = \|\mathbf{r}_j - \mathbf{r}_{j-1}\|$ je geometrická délka úseku v $j$-té oblasti a $\sigma^{(m, ij)}$ je tamní makroskopický průřez. Celková optická vzdálenost do cílového bodu:

Distribuční funkci volné dráhy přes rozhraní materiálů zapisujeme multiplikativně:

Algoritmicky se to implementuje tak, že místo „uleť konstantní vzdálenost“ kontrolujeme, kdy částice překračuje rozhraní mezi materiály, a v každém úseku používáme příslušný $\sigma^{(m)}$. Geometrický modul simulace musí umět rychle detekovat protnutí trajektorie s povrchy oblastí.

Co se v simulaci vyhodnocuje (skóre)

Cílem MC transportu obvykle není jednotlivá trajektorie (ta je sama o sobě k ničemu), ale statistický odhad nějaké veličiny:

• Tok částic přes danou plochu (kolik částic za sekundu prošlo detektorem).
• Absorbovaná dávka v dané oblasti (kolik energie zde uvízlo).
• Distribuce energie částic v detektoru.
• Účinnost stínění (poměr prošlého toku k toku zdroje).
• Multiplikační faktor u štěpitelných materiálů (poměr počtu neutronů v generaci $n+1$ k počtu v generaci $n$).

Co je teorie hromadné obsluhy a proč ji potřebujeme

Teorie hromadné obsluhy (anglicky queueing theory, teorie front) studuje systémy, do kterých přicházejí požadavky a které je musí obsloužit. Pojem požadavku je obecný: zákazník v bance, paket v síťovém routeru, výrobek na lince, telefonát v call centru, čtenář v knihovně. Pojem obsluhy také: prodavač, server, robot na lince, operátor.

Cílem MC simulace je najít optimální nastavení systému. Konkrétně:

• Kolik obslužných linek potřebujeme, aby se fronta nezahltila?
• Jaký režim řazení požadavků (FIFO, LIFO, prioritní) je nejvhodnější?
• Jak dlouhý zásobník (frontu) mít? Co se stane, když je plný?
• Jak se chovat v špičkách, kdy intenzita příchodů kolísá?

Aplikace jsou všudypřítomné: call centra (kolik operátorů potřebujeme v dané hodině), banky (kolik přepážek), nemocnice (kolik lůžek na pohotovosti), webové servery (kolik instancí), výrobní linky (kolik strojů paralelně), letecké terminály (počet odbavovacích přepážek). Často je v sázce hodně peněz: každá nadbytečná linka stojí, každý ztracený zákazník je škoda.

Klíčové náhodné veličiny

Při simulaci pracujeme se dvěma fundamentálními náhodnými veličinami, jejichž rozdělení pravděpodobnosti musíme znát buď z teorie (analyticky) nebo z empirických dat (měření v reálném systému):

1. Čas příchodu události (požadavku) do systému, respektive interval mezi příchody $\eta_i$.
2. Doba potřebná k obsluze jednoho požadavku, označovaná $\tau_z$.

Sekundární náhodné veličiny:

• Doba čekání $\tau_p$ - doba, po kterou je požadavek ochotný čekat ve frontě, než z ní odejde neuspokojený.
• Atributy požadavku jako priorita, kategorie, velikost (pro pakety v síti), atd.

Popis systému: tři komponenty

Každý systém hromadné obsluhy se popisuje třemi částmi:

• Vstupní tok požadavků: popsán časy příchodů $t_i$, prioritou požadavků a tolerovanými dobami čekání $\tau_p$.
• Struktura systému: počet obslužných linek (uzlů) a kapacita zásobníku (maximální délka fronty).
• Režim: pravidla pro řazení požadavků do front a způsob jejich výběru. Klasické varianty:
  • FIFO (First In, First Out): kdo přijde dřív, dřív se obslouží.
  • LIFO (Last In, First Out): zásobníková fronta, poslední se obsluhuje první.
  • Prioritní fronta: požadavky s vyšší prioritou se předbíhají před nižší.
  • SJF (Shortest Job First): kratší úlohy mají přednost před dlouhými.

Každý jednotlivý požadavek je v simulaci popsán vektorem parametrů $(t_i, a_1, \dots, a_m)$, kde:

• $i$ je pořadové číslo požadavku.
• $t_i$ je čas jeho příchodu.
• $a_1, \dots, a_m$ jsou další atributy (priorita, typ, velikost úlohy, ...).

Kategorie systémů podle doby čekání $\tau_p$

Vlastnosti toku požadavků

Označme $\eta_i = t_i - t_{i-1}$ pro $i > 1$ a $\eta_1 = t_1$. To jsou intervaly mezi po sobě jdoucími příchody. Sledujeme zejména pravděpodobnost $P(k, t_0, t)$, že v časovém okně $\langle t_0, t_0 + t \rangle$ dorazí právě $k$ požadavků.

Typy toků

• Homogenní tok: jediným sledovaným parametrem požadavku je čas příchodu $t_i$. Žádné priority, žádné kategorie. V opačném případě je tok nehomogenní.
• Tok s omezenou závislostí: intervaly $\eta_i$ jsou nezávislé náhodné veličiny. Jejich sdružená hustota je součinem podmíněných hustot:
  $g_s(z_1, \dots, z_s) = \displaystyle\prod_{k=1}^{s} f_k(z_k)$
  kde $f_k$ je hustota $k$-tého intervalu za podmínky, že na začátku intervalu dorazil předchozí požadavek.
• Stacionární tok: počet požadavků $P(k, t)$ závisí pouze na délce intervalu $t$, ne na startovním čase $t_0$. Statistické vlastnosti se v čase nemění. Podmíněné hustoty jsou si rovny: $f_2 = f_3 = \dots = f$. Normovací konstanta $\lambda$, která vyjadřuje celkovou hustotu (intenzitu) toku, je převrácenou hodnotou střední délky intervalu:
  $\lambda = \left(\displaystyle\int_0^\infty x\,f(x)\,dx\right)^{-1}$
• Ordinární tok: pravděpodobnost, že by v extrémně krátkém čase ($t \to 0$) přišel více než jeden požadavek najednou, je zanedbatelná. Formálně $P(\ge 2 \text{ příchody v } \Delta t) = o(\Delta t)$. To znamená, že požadavky přicházejí „po jednom“, ne v dávkách.

Poissonův tok: nejdůležitější aproximace

V praxi se nejčastěji setkáváme s tokem, kde intervaly mezi příchody mají exponenciální rozdělení s hustotou $f(x) = \lambda e^{-\lambda x}$. Z toho automaticky plyne, že počet příchodů v daném intervalu má Poissonovo rozdělení.

Důvod, proč je Poissonův tok tak rozšířený: je to jediný tok, který současně splňuje stacionárnost, ordinaritu a nezávislost přírůstků. Centrální limitní věta v podobném smyslu „přitahuje“ reálné toky k Poissonovi.

Stacionární Poissonův tok

Pokud je intenzita $\lambda$ konstantní v čase:

Nestacionární Poissonův tok

V reálu intenzita kolísá: v bance přijde poledne víc lidí, na webové aplikaci je špička večer, server má jiný provoz v pracovní dny a o víkendu. Pak se intenzita stane funkcí času $\lambda(u)$ a vzorec se zobecní:

$\Lambda(t_0, t)$ je střední počet požadavků v daném intervalu, $\lambda(u)$ je okamžitá hustota toku.

Dva přístupy k simulaci v MC

Při implementaci máme dvě možnosti, jak zacházet s časem:

Metoda rozštěpení (splitting) pro vzácné jevy

Některé úlohy se zaměřují na kritické stavy, které nastávají velmi zřídka. Příklady:

• Pravděpodobnost, že v telefonní centrále současně čeká více než 1000 hovorů.
• Pravděpodobnost, že webový server přečerpá kapacitu paměti.
• Pravděpodobnost, že čekací doba na pohotovosti přesáhne kritickou mez.

U takových jevů selhává naivní MC: museli bychom simulovat miliardy běhů, než by jev vůbec jednou nastal. Náhrada: metoda rozštěpení (splitting), která uměle „posiluje“ pravděpodobnost zájmového jevu váženými klony trajektorií.

Co se v simulaci vyhodnocuje

Po dostatečném počtu nezávislých běhů MC simulace odhadneme řadu praktických veličin:

• Průměrná délka fronty a její maximum (peak).
• Průměrná doba čekání a doba obsluhy.
• Procento odmítnutých / odcházejících požadavků (důležité u systémů s odmítáním nebo smíšených).
• Vytížení linek (utilizace) v procentech času, kdy byly obsazené. Cílem je obvykle utilizace 70–80 %: nižší znamená plýtvání kapacitou, vyšší vede k frontám.
• Pravděpodobnost vzácných stavů (přes splitting).
• Distribuce doby pobytu v systému (čekání + obsluha).
• Náklady na provoz (počet linek × hodinová sazba) versus kvalita služby (procento spokojených zákazníků).

Souvislost s ostatními otázkami

Otázka 5 se silně překrývá s ostatními:

• Generování intervalů příchodu používá inverzní transformaci pro exponenciální rozdělení (Otázka 3): $\eta = -\tfrac{1}{\lambda}\ln \Upsilon$.
• Generování počtu příchodů v intervalu používá Poissonův algoritmus (Otázka 3).
• Metoda splitting je strukturálně podobná zpracování štěpení v transportu částic (Otázka 4): větvení trajektorií s vahami, zásobník na nezpracované větve.
• Vše stojí na kvalitním PRNG pro $R(0,1)$ (Otázky 1 a 2).