Metody Monte Carlo: Studijní materiál ke státnicím

Co je tu uvnitř. Pět otázek ke státnicím z Monte Carla. U každé je nejprve „kostra“, kterou musí zkoušející slyšet, pak technické detaily, vzorce a názorné příklady.

Jak to číst. Modré rámečky = lidská intuice. Žluté = příklad k zapamatování. Oranžové = časté chytáky. Zelené = co si odnést.

Seznam státnicových otázek

Popište metodu Monte Carlo. Vysvětlete obvyklé způsoby modelování a generování rovnoměrně rozdělených pseudonáhodných čísel.
Popište metodu Monte Carlo. Vysvětlete způsoby testování kvality generátorů pseudonáhodných čísel. Popište vlastnosti generátorů a porovnejte je.
Popište metodu Monte Carlo. Popište obecné principy generování obecně rozdělených náhodných čísel. Vyberte si dvě rozdělení a popište způsoby jejich generování.
Popište principy modelování transportu částic metodou Monte Carlo.
Popište principy modelování systémů hromadné obsluhy metodou Monte Carlo.

Otázka 1

Metoda Monte Carlo. Modelování a generování rovnoměrně rozdělených pseudonáhodných čísel.

1. Co je metoda Monte Carlo

Monte Carlo (MC) je výpočetní metoda, která modeluje náhodné děje a získává tak odhad veličiny, kterou by jinak bylo těžké nebo nemožné spočítat analyticky. Funguje i pro úlohy, které náhodné nejsou (např. výpočet integrálu), náhodu si tam záměrně přivedeme jako nástroj.

Tři kroky, ze kterých MC stojí:

Znát rozdělení pravděpodobnosti všech jevů v modelu (z teorie nebo empiricky).
Vygenerovat realizace náhodné veličiny, typicky přes $\Upsilon \sim R(0,1)$.
Transformovat $\Upsilon$ na požadované rozdělení a výsledky statisticky zpracovat.

Přesnost roste s počtem simulací (typicky chyba ~$1/\sqrt{N}$), takže často potřebujeme miliony opakování.

2. Zdroje náhodných čísel

Tabulky náhodných čísel, historické, dnes prakticky bezvýznamné.
Hardwarové generátory (bílý šum, šum diody, kvantové jevy). Skutečně náhodné, ale nereprodukovatelné, nelze opakovat tentýž výpočet.
Softwarové (pseudonáhodné) generátory: deterministické posloupnosti, které se chovají „jako náhodné“. Vlastnosti jdou matematicky studovat a výpočet jde zopakovat se stejným seedem. Dnes standard.

3. Požadavky na pseudonáhodnou posloupnost $\theta = \{x_i\}_{i=0}^{\infty}$

Posloupnost musí být statisticky nerozlišitelná od výběru z $Y \sim R(0,1)$. Z toho plyne:

Střední hodnota: $\displaystyle \lim_{n\to\infty}\frac{1}{n}\sum_{i=1}^{n} x_i = \frac{1}{2}$ (protože $EY = 1/2$).
Rozptyl: $\displaystyle \lim_{n\to\infty}\frac{1}{n}\sum_{i=1}^{n}(x_i-\bar{x})^2 = \frac{1}{12}$.
Rovnoměrnost v 1D: Pro libovolnou podmnožinu $I_v \subset (0,1)$ s objemem $v$ platí $\displaystyle \lim \frac{N_v}{n} = v$. (10 % intervalu → 10 % bodů.)
$k$-rovnoměrnost: $k$-tice po sobě jdoucích čísel musí rovnoměrně vyplnit jednotkovou krychli $C_k = \langle 0,1 \rangle^k$. Tohle je klíčové, generátor může být dobrý v 1D, ale fatálně selhat ve 2D nebo 3D.
Nulová autokorelace: $\displaystyle \lim \frac{1}{n}\sum (x_{i-t}-\tfrac{1}{2})(x_i-\tfrac{1}{2}) = 0$ pro každé $t > 0$. Znalost dřívějšího čísla nesmí pomoct předpovědět další.

Problém s implementací: Počítač zobrazí jen konečný počet reálných čísel, a navíc nerovnoměrně. Proto se počítá v celých číslech $y_i \in \{0,\dots,M-1\}$ (rychlejší aritmetika) a výsledek se získá jako $x_i = y_i / M$.

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

$x_{n+1} = (K \cdot x_n + B) \bmod M$

$x_0$ (seed): počáteční hodnota, určuje celou budoucí řadu.
$K, B, M$: konstanty (multiplikátor, inkrement, modul). Volí se pečlivě, často $M = 2^n$ kvůli rychlosti.
Výsledek: celé číslo $0 \le x_{n+1} < M$, převod na $(0,1)$ dělením $M$.

RANDU (klasický odstrašující případ): $M=2^{31}$, $K=65539$, $B=0$. V jednom rozměru vypadal v pohodě, ale když se vzaly trojice po sobě jdoucích čísel jako body v 3D, ležely jen na 15 nadrovinách. Naprosto kazil simulace, dlouhá léta se ale používal.

Whichmann–Hillův generátor: kombinace (součet modulo 1) čtyř LCG. Lepší než jeden LCG, ale stále má problémy s $k$-rovnoměrností a korelací.

5. Mersenne Twister (MT19937)

Dnes nejrozšířenější PRNG (Python, R, PHP, MATLAB). Založen na lineárních posuvných registrech.

$x_{k+n} = x_{k+m} \oplus \bigl(x_k^{\{u(w-r)\}} \,|\, x_{k+1}^{\{l(r)\}}\bigr)\,A$

$w$: délka slova (typicky 32 bitů).
$n = 624$: počet stavů v registru.
$m = 397$, $r = 31$: posuny a maskování.
$A$: konstantní matice $w \times w$ z $\{0,1\}$, realizovaná jako shift+XOR.
$\oplus$ je bitový XOR, $|$ je zřetězení.

Algoritmus „slepí“ horní bity slova $x_k$ s dolními bity $x_{k+1}$, projede je rychlou bitovou transformací (matice $A$) a sečte XORem s dávno použitým slovem $x_{k+m}$. To dramaticky prodlouží periodu a rozhází bity rovnoměrně.

Proč je tak dobrý:

Perioda $2^{19937}-1$, Mersennovo prvočíslo, odsud název.
$k$-rovnoměrnost až ve 623 rozměrech.
Rychlý: jen XOR a posuny.

Mersenne Twister není kryptograficky bezpečný. Když útočník vidí 624 po sobě jdoucích výstupů, zrekonstruuje vnitřní stav a předpoví všechno dál. Pro kryptografii se používají CSPRNG (např. ChaCha20).

6. Generování na obecný interval $R(a,b)$

$\Upsilon_{ab} = (b-a)\cdot \Upsilon + a, \qquad \Upsilon \sim R(0,1)$

Délka hovoru mezi 2 a 10 minutami. Padne $\Upsilon = 0.3$. Výsledek: $(10-2)\cdot 0.3 + 2 = 4{,}4$ min.

7. Golombovy požadavky (binární posloupnosti)

Nutná (ne dostačující) podmínka rovnoměrnosti pro binární posloupnosti:

Rovnováha: počet 0 a 1 se liší max. o jedna.
Běhy (runs): polovina běhů má délku 1, čtvrtina délku 2, osmina délku 3, …
Ideální autokorelace: při posunu vůči sobě je shoda minimální.

Princip MC = (a) máme rozdělení, (b) přes $\Upsilon \sim R(0,1)$ vyrábíme realizace, (c) opakujeme miliony krát. Pamatuj LCG (jednoduchý + nadroviny) a MT (perioda, ne pro kryptografii). K-rovnoměrnost je důvod, proč nestačí testovat průměr.

Otázka 2

Testování kvality generátorů, vlastnosti generátorů a jejich porovnání.

Pokud posloupnost selže v některém testu, není to model $R(0,1)$. Pokud projde všemi, ještě to neznamená, že je správná (testy jsou nutná, ne postačující podmínka). Hladina významnosti je obvykle $\alpha = 0{,}05$.

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

Ověřuje, že v binárním zápisu se 0 a 1 vyskytují přibližně stejně často.

$X_1 = \dfrac{(n_0 - n_1)^2}{n}$

Rozdělení: $\chi^2$ s 1 stupněm volnosti (známe-li $n_0$, $n_1 = n - n_0$).
Kritická hodnota: $X_1 < 3{,}84$ s 95% spolehlivostí.

Vygenerujeme 1000 bitů, najdeme 550 jedniček a 450 nul. $X_1 = \frac{(450-550)^2}{1000} = 10$. Protože $10 > 3{,}84$, generátor v tomto testu padá.

Číslo 3,84 je 95% kvantil $\chi^2_1$. Plocha pod křivkou rozdělení od 0 do 3,84 = 95 %. Vše nad 3,84 už spadá do „podezřelých 5 %“.

Pohled přes binomické rozdělení

Pro $m$ slov o $n$ bitech má počet jedniček binomické rozdělení s parametry $mn$ a $p = 1/2$. S 95% pravděpodobností by měl ležet v intervalu $\Bigl(\frac{nm}{2} - \sqrt{nm},\ \frac{nm}{2} + \sqrt{nm}\Bigr)$.

Test sérií (rovnoměrnost dvojic)

Frekvenční test neodhalí třeba posloupnost 010101..., ta má stejně 0 i 1, ale je periodická. Test sérií zkoumá, zda dvojice 00, 01, 10, 11 jsou přibližně stejně časté.

$X_2 = \dfrac{4}{n-1}(n_{00}^2 + n_{01}^2 + n_{10}^2 + n_{11}^2) - \dfrac{2}{n}(n_0^2 + n_1^2) + 1$

Rozdělení: $\chi^2$ se 2 stupni volnosti (4 typy dvojic, vazby na celkový počet).
Kritická hodnota: $X_2 < 5{,}99$.

Druhá varianta: rozdělí se $(0,1)$ na $r$ stejných podintervalů a sleduje četnost padání čísel do nich, připravuje testy ve více dimenzích.

Pokerový test

Sleduje četnost vzorů (úseků) délky $m$ bitů. Pokud by některá kombinace byla častější, generátor je předvídatelný.

Vyber $m$; rozděl posloupnost délky $N$ na $k = \lfloor N/m \rfloor$ úseků.
Existuje $2^m$ typů úseků; $n_i$ = počet výskytů $i$-tého typu.
Podmínka průkaznosti: $k \ge 5 \cdot 2^m$ (na každý typ aspoň 5 očekávaných výskytů).

$X_3 = \dfrac{2^m}{k} \sum_{i=1}^{2^m} n_i^2 \;-\; k$

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

Pro $m = 3$ máme $2^3 = 8$ typů trojic a 7 stupňů volnosti. Pokud generátor vyplivuje příliš mnoho „111“, $X_3$ rychle naroste a test selže.

Test běhů (Runs test)

Sleduje, zda délky nepřerušených sekvencí (běhů) bitů odpovídají Golombovým požadavkům.

Blok $B_i$: nepřerušená řada jedniček délky $i$.
Mezera $M_i$: nepřerušená řada nul délky $i$.

$e_i = \dfrac{n - i - 3}{2^{i+2}}$ (očekávaný počet běhů délky $i$)

Vezmeme největší $k$ takové, že $e_k \ge 5$.

$X_4 = \displaystyle\sum_{i=1}^{k} \dfrac{(B_i - e_i)^2}{e_i} + \sum_{i=1}^{k} \dfrac{(M_i - e_i)^2}{e_i}$

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

Hodíš mincí 10× a padne pořád orel? „Běh“ délky 10. Test očekává, že polovina běhů má délku 1, čtvrtina délku 2, atd. Když tyto poměry generátor nedrží, testem neprojde.

Autokorelační test

Hledá skryté závislosti mezi bity vzdálenými o $d$ pozic.

Zvolí se posun $1 \le d \le \lfloor n/2 \rfloor$ (různý od periody).
Spočítá se $A(d) = \sum_{i=1}^{n-d-1} s_i \oplus s_{i+d}$ (kolikrát se bity liší).

$X_5 = 2\,\dfrac{A(d) - \frac{n-d}{2}}{\sqrt{n-d}}$

Rozdělení: $N(0,1)$ pro $n - d > 10$. Generátor projde, pokud $|X_5| \le 1{,}96$ (95% spolehlivost).

Maurerův univerzální statistický test

Nejsilnější z uvedených. Vychází z informační teorie: opravdu náhodná data nesmí jít zkomprimovat. Pokud lze, generátor obsahuje skrytý řád.

Univerzální = postihne široké spektrum nedostatků, které jednodušší testy minou.
Velmi dobře odhalí dlouhé periodicity.
Průchod Maurerovým testem je silný indikátor, že posloupnost vypadá jako $R(0,1)$.

Shrnutí kvality: LCG vs. Mersenne Twister

Vlastnost	LCG	Mersenne Twister
Složitost implementace	Velmi jednoduchá (jedno násobení + modulo)	Složitější (posuvné registry, matice $A$)
Perioda	Krátká, max $M$	$2^{19937} - 1$
$k$-rovnoměrnost	Selhává, body leží na nadrovinách	Dobrá až v 623 rozměrech
Korelace	Vyšší, viditelné vzory	Velmi nízká
Kryptografická bezpečnost	Ne	Ne (predikce po 624 výstupech)
Test průchodnost	Často padá v komplexních testech	Prochází téměř vším

Zapamatuj si 5 typů testů (frekvenční, sérií, pokerový, běhů, autokorelační) + Maurer. Kritické hodnoty $\chi^2$: 1 stupeň → 3,84; 2 stupně → 5,99. Plus principiální tabulku LCG vs. MT.

Otázka 3

Generování obecně rozdělených náhodných čísel. Příklady dvou rozdělení.

Princip je vždy stejný: máme „spravedlivé“ $\Upsilon \sim R(0,1)$ a chceme z něj vyrobit realizaci náhodné veličiny s úplně jiným rozdělením. Cestou je distribuční funkce, slouží jako převodník.

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

Chceme $\xi$ s rozdělením $P(\xi = x_k) = p_k$. Geometricky: rozřežeme úsečku $(0,1)$ na kousky o délkách $p_k$ a sledujeme, do kterého padne $\Upsilon$.

Základní algoritmus

Generuj $\Upsilon \sim R(0,1)$; polož $S \leftarrow \Upsilon$, $k \leftarrow 0$.
$S \leftarrow S - p_k$.
Je-li $S > 0$, zvětši $k$ o 1 a jdi na krok 2.
Výsledek: $\xi = x_k$.

Složitost

$\tau = \tau_R + (\tau_p + \tau_+)(1 + E\xi)$

Časová náročnost lineárně závisí na střední hodnotě indexu $E\xi$. Čím dál se průměrně musí algoritmus dostat, tím déle počítá.

Představ si „rozřezaný metr“. Pravděpodobnosti jsou kousky, $\Upsilon$ je drobek, který někam dopadl. Algoritmus postupně odřezává kousky, dokud drobek „nepropadne pod nulu“, tam je odpověď.

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

Setřídíme $p_{k_0} \ge p_{k_1} \ge p_{k_2} \ge \dots$ a pak běží stejný algoritmus. Nejčastější hodnoty se trefíme hned napoprvé.

Mám rozdělení Pivo $p_0 = 0{,}1$, Voda $p_1 = 0{,}2$, Rum $p_2 = 0{,}7$.
Bez setřídění: $E\xi = 0\cdot 0{,}1 + 1\cdot 0{,}2 + 2\cdot 0{,}7 = 1{,}6$, průměr 2,6 dotazů.
Po setřídění (Rum, Voda, Pivo): $E\xi = 0\cdot 0{,}7 + 1\cdot 0{,}2 + 2\cdot 0{,}1 = 0{,}4$, průměr 1,4 dotazu, skoro 2× rychleji.

Modifikovaný algoritmus 2: start od vrcholu

Pro rozdělení, kde $p_k$ rostou až k indexu $h$ a pak klesají (typicky Poissonovo s velkým $\lambda$). Předpočítáme $Q = \sum_{k=0}^{h} p_k$ a startujeme od $h$:

Generuj $\Upsilon$. Polož $p \leftarrow p_h$, $m \leftarrow h$, $S \leftarrow \Upsilon - Q$.
Když $S > 0$, hledej napravo ($m+1, m+2, \dots$).
Když $S \le 0$, hledej nalevo ($m-1, m-2, \dots$).

Vrtulník tě vysadí přímo na vrcholu kopce. Ze vstupenky $\Upsilon$ poznáš, jestli máš jít kousek nahoru, nebo dolů. Pro $\lambda = 1000$ ušetříš ~1000 iterací.

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

1. Bernoulli s parametrem $p$. Generuj $\Upsilon$. Je-li $\Upsilon < p$, výsledek 1; jinak 0.

Šance na výhru 30 % ($p=0{,}3$). Vše od 0,00 do 0,29 → vyhrál jsi.

2. Rovnoměrné na $\{1, \dots, n\}$. Generuj $\Upsilon$. Spočítej $s = \lfloor \Upsilon \cdot n \rfloor + 1$.

Kostka ($n=6$). $\Upsilon = 0{,}75$. $0{,}75 \cdot 6 = 4{,}5$, celá část 4, +1 = 5.

3. Poissonovo s parametrem $\lambda$.

Nastav $n \leftarrow -1$, $S \leftarrow 1$, $q \leftarrow e^{-\lambda}$.
Generuj $\Upsilon$.
$n \leftarrow n+1$, $S \leftarrow S \cdot \Upsilon$.
Je-li $S > q$, vrať se na krok 2. Jinak výsledek je $\xi = n$.

Trik: namísto sčítání času se násobí náhodná čísla < 1. Hodnota $S$ klesá, dokud nepropadne pod limit $q$.

Kavárna a Poisson, $\lambda = 2$ (2 zákazníci/hod). Limit $q = e^{-2} \approx 0{,}1353$.

Kolo 1: $\Upsilon_1 = 0{,}60$, $S = 0{,}60$. Je $0{,}60 > 0{,}1353$ → ANO. $n=0$.
Kolo 2: $\Upsilon_2 = 0{,}45$, $S = 0{,}60 \cdot 0{,}45 = 0{,}27$. ANO. $n=1$.
Kolo 3: $\Upsilon_3 = 0{,}30$, $S = 0{,}27 \cdot 0{,}30 = 0{,}081$. NE → STOP.

Výsledek: 2 zákazníci v této hodině.

B) Spojité náhodné veličiny

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

Distribuční funkce $F(x)$ mapuje hodnoty na interval $(0,1)$. Položíme $F(\xi) = \Upsilon$ a vyřešíme:

$\xi = F^{-1}(\Upsilon)$

Distribuční funkce je kopec stoupající od 0 do 1. „Vystřelíš“ vodorovně z osy $y$ z hodnoty $\Upsilon$ a tam, kde se trefíš do kopce, projektuješ dolů na osu $x$, to je $\xi$.

Exponenciální rozdělení $f(t) = \lambda e^{-\lambda t}$, $F(t) = 1 - e^{-\lambda t}$. Řešení $F(\xi) = \Upsilon$ dá $\xi = -\frac{1}{\lambda} \ln(1 - \Upsilon)$. Protože $1 - \Upsilon$ má stejné rozdělení jako $\Upsilon$, píšeme:

$\xi = -\dfrac{1}{\lambda} \ln \Upsilon$

Tento vzorec je nesmírně užitečný, exponenciální čekací doby jsou srdcem teorie front i transportu částic.

2. Zamítací metoda (Rejection)

Když rovnici $F(\xi) = \Upsilon$ nejde rozumně vyřešit. Předpoklad: hustota $f(x)$ je omezená konstantou $C$ na intervalu $\langle a, b \rangle$.

$\Upsilon_1 \sim R(a,b)$: náhodná pozice na ose $x$.
$\Upsilon_2 \sim R(0,C)$: náhodná výška.
Je-li $\Upsilon_2 > f(\Upsilon_1)$ (bod nad křivkou), zamítni a zkus znovu.
Jinak $\xi = \Upsilon_1$.

Stříleli bys do obdélníka přes grafem hustoty. Trefíš pod křivku → zásah platí. Trefíš nad → zahodíš a střílíš znovu. Efektivita závisí na tom, jak „těsně“ obdélník obaluje křivku.

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

A) Aproximace přes centrální limitní větu.

$\eta = \displaystyle\sum_{i=1}^{12} \Upsilon_i - 6$

Volíme $n=12$, protože tehdy $\sqrt{12/n} = 1$ a celý vzorec se zjednoduší. Dostaneme dobrou (ne přesnou) aproximaci $N(0,1)$.

B) Box–Mullerova transformace (přesná). Ze dvou $\Upsilon_1, \Upsilon_2 \sim R(0,1)$:

$\xi_1 = \sqrt{-2\ln\Upsilon_1}\cdot \sin(2\pi\Upsilon_2)$
$\xi_2 = \sqrt{-2\ln\Upsilon_1}\cdot \cos(2\pi\Upsilon_2)$

$\xi_1$ a $\xi_2$ jsou dvě nezávislé realizace $N(0,1)$. Žádné sčítání 12 čísel, navíc přesné.

Tři páky: inverzní transformace (umíš $F^{-1}$), zamítací metoda (neumíš $F^{-1}$, ale znáš $f$), a pro $N(0,1)$ Box–Muller nebo CLV se součtem 12 čísel. Klíčové diskrétní příklady: Bernoulli, rovnoměrné na $\{1,\dots,n\}$, Poisson.

Otázka 4

Modelování transportu částic metodou Monte Carlo.

Simulace průchodu záření (fotony, elektrony, neutrony) hmotou. Typické úlohy: výpočet stínění, optimalizace jaderných zařízení, dozimetrie, simulace detektorů. Trajektorie každé částice je posloupností náhodných událostí: uletí kus → narazí na atom → buď je pohlcena, rozptýlena, nebo způsobí štěpení.

Zjednodušující předpoklady

Mezi částicemi není interakce: platí, pokud hustota toku není extrémní.
Pohyb je popsán pouze interakcemi s atomy prostředí; potřebujeme znát pravděpodobnosti (pohlcení, rozptyl, štěpení) a hustotu atomů.
Mezi interakcemi částice letí přímočaře (zakřivení gravitací atd. zanedbáno).
Sekundární částice vznikají v bodě interakce.

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

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

$g(\mathbf{r}, \boldsymbol{\omega}', \boldsymbol{\omega}, E) = \dfrac{d\sigma_a}{d\Omega}$

$\mathbf{r}$: bod v prostoru.
$\boldsymbol{\omega}', \boldsymbol{\omega}$: směr před a po interakci.
$E$: energie částice.
$d\Omega = \sin\vartheta\, d\vartheta\, d\varphi$: element prostorového úhlu.
$\sigma_a$: mikroskopický průřez (pravděpodobnost na jeden atomový terč, normovaná na jednotkovou hustotu toku).

Makroskopický účinný průřez

$\sigma_a^{(m)}(\mathbf{r}) = \sigma_a \cdot N$

$N$ je počet terčů (atomů) na jednotku objemu. Vyjadřuje skutečnou pravděpodobnost interakce v daném materiálu.

Volná dráha a optická vzdálenost

Vzdálenost mezi dvěma interakcemi (volná dráha $\ell$) je náhodná veličina s exponenciálním rozdělením:

$f_\ell(t) = \sigma^{(m)}(\mathbf{r}(t))\cdot e^{-\int_0^t \sigma^{(m)}(\mathbf{r}(u))\,du}$

Distribuční funkce:

$F_\ell(t) = 1 - e^{-\tau(t)}, \qquad \tau(t) = \displaystyle\int_0^t \sigma^{(m)}(\mathbf{r}(u))\,du$

$\tau(t)$ je optická vzdálenost, integrál makroskopického průřezu po dráze. Pravděpodobnost, že částice přežije dráhu délky $t^*$, je $e^{-\tau(t^*)}$.

V homogenním prostředí má volná dráha exponenciální rozdělení s parametrem $\sigma^{(m)}$, takže ji generuješ vzorcem $s = -\frac{1}{\sigma^{(m)}} \ln \Upsilon$, přesně otázka 3, inverzní transformace.

Typy interakce: střední hodnota $\nu$

Interakce	$\nu$	Co se stane
Pohlcení (absorpce)	$\nu = 0$	Částice zaniká, sledování končí.
Rozptyl	$\nu = 1$	Mění směr a/nebo energii, letí dál.
Štěpení / multiplikace	$\nu > 1$	Vznikne víc nových částic.

Vznikají-li identické sekundární částice, sloučíme všechny dílčí interakce do jediné s váženými průměry:

$\sigma^{(m)} = \sum_i \sigma_i^{(m)}, \quad g = \dfrac{1}{\sigma^{(m)}}\sum_i \sigma_i^{(m)} g_i, \quad \nu = \dfrac{1}{\sigma^{(m)}}\sum_i \sigma_i^{(m)} \nu_i$

Algoritmus pro $\nu \le 1$ (pohlcení + rozptyl)

Emise: vygeneruj $\mathbf{r}_0$ a směr $\boldsymbol{\omega}_0$ podle rozdělení zdroje.
Transport: nasimuluj volnou dráhu $s$ (exp. rozdělení) a posuň částici: $\mathbf{r} = \mathbf{r}_0 + s\,\boldsymbol{\omega}_0$.
Test geometrie: je $\mathbf{r}$ mimo zkoumanou oblast? Pak částice utekla, jdi na 8.
Druh interakce: diskrétní výběr podle průřezů. Pokud pohlcení, jdi na 8.
Sekundární parametry: z diferenciálních průřezů simuluj úhel rozptylu $\vartheta$ ($\mu = \cos\vartheta$) a azimut $\varphi$.
Nový směr: přepočti směrový vektor $\boldsymbol{\omega}$.
Cyklus: $\mathbf{r}_0 \leftarrow \mathbf{r}$, $\boldsymbol{\omega}_0 \leftarrow \boldsymbol{\omega}$, jdi na 2.
Vyhodnocení: připočti skóre (např. dávku, počet rozptylů), připrav novou primární částici.

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

Počítač neumí sledovat dvě částice najednou, potřebujeme zásobník (kontejner).

Počet sekundárních částic $n$ se losuje z $\lfloor \nu \rfloor$ s pravděpodobností $1 - \{\nu\}$ a $\lfloor \nu \rfloor + 1$ s pravděpodobností $\{\nu\}$.
První částici simulujeme dál běžným způsobem; ostatní (2. až $n$-tou) uložíme do zásobníku.
Když aktuálně sledovaná částice zanikne nebo opustí oblast, vyzvedneme další ze zásobníku.
Výpočet pro jednu primární částici končí, až je zásobník prázdný.

Štěpení = strom trajektorií. Kontejner je „TODO list“ větví, které musíme ještě projít. Hloubkové prohledávání, jedna větev do konce, pak další.

Nehomogenní prostředí

Reálně částice prochází více materiály (Fe → beton → vzduch). Modelujeme po částech homogenní prostředí.

Optická vzdálenost se počítá sčítáním přes úseky:

$\tau_{ij} = \lambda_{ij} \cdot \sigma^{(m,ij)}$

kde $\lambda_{ij} = \|\mathbf{r}_j - \mathbf{r}_{j-1}\|$ je délka úseku v $j$-té oblasti. Celkem:

$\tau(\mathbf{r}) = \displaystyle\sum_{j=1}^{s-1} \tau_{ij} + \lambda_{i_s} \cdot \sigma^{(m,i_s)}$

A distribuční funkce volné dráhy je multiplikativní přes rozhraní:

$F_\ell(t) = 1 - e^{-\lambda_{i_s}\sigma^{(i_s)}} \cdot \displaystyle\prod_{j=1}^{s-1} e^{-\tau_{ij}}$

Stínění reaktoru. Foton vyletí z paliva, prochází vrstvou olova (vysoké $\sigma^{(m)}$, krátké volné dráhy, hodně rozptylů), pak betonem (střední $\sigma^{(m)}$) a nakonec vzduchem (malé $\sigma^{(m)}$, letí dlouho rovně). MC simulace odhadne, kolik procent fotonů projde skrz a s jakou energií.

Klíčové: účinný průřez = pravděpodobnost interakce, volná dráha = exponenciální rozdělení, $\nu$ určuje větvení. Pohlcení končí, rozptyl pokračuje, štěpení plní zásobník. Nehomogenita se řeší po úsecích.

Otázka 5

Modelování systémů hromadné obsluhy metodou Monte Carlo.

Teorie front. Cílem simulace je najít optimální nastavení systému, počet linek, organizaci front, prioritizaci, aby systém zvládal nápor, ale ne zbytečně přeplýtval kapacitou. Aplikace: call centra, banky, IT servery, fronty u kasy, sítě, výrobní linky.

Stavební kameny systému

Tok požadavků: posloupnost událostí, které do systému vstupují.
Linka (uzel): zařízení, které požadavek obsluhuje.
Doba obsluhy $\tau_z$: jak dlouho zpracování trvá.
Doba čekání $\tau_p$: jak dlouho je požadavek ochotný čekat.

Náhodné veličiny, jejichž rozdělení potřebujeme znát (z teorie nebo empiricky):

Časy příchodu událostí.
Doby obsluhy.

Popis systému

Vstupní tok: časy příchodů, priority, doby čekání.
Struktura: počet linek, kapacita zásobníku (fronty).
Režim: pravidla pro zařazení do front (FIFO, LIFO, prioritní, …).

Kategorie podle doby čekání

Typ	$\tau_p$	Chování
S odmítáním	$\tau_p = 0$	Není-li ihned volná linka, požadavek odchází neobsloužen.
Čekací	$\tau_p = +\infty$	Požadavek čeká, dokud na něj nepřijde řada.
Smíšený	$0 < \tau_p < +\infty$	Konečná tolerance, po vypršení odchází.

Jeden požadavek je popsán vektorem $(t_i, a_1, \dots, a_m)$, kde $i$ je pořadové číslo, $t_i$ čas příchodu a $a_j$ další atributy (např. priorita).

Vlastnosti toku požadavků

Označme $\eta_i = t_i - t_{i-1}$ interval mezi po sobě jdoucími příchody (pro $i=1$ je $\eta_1 = t_1$). Sledujeme pravděpodobnost $P(k, t_0, t)$, že v intervalu $\langle t_0, t_0 + t \rangle$ dorazí právě $k$ požadavků.

Homogenní tok: jediným parametrem je čas příchodu $t_i$ (žádné priority, kategorie).
Tok s omezenou závislostí: intervaly $\eta_i$ jsou nezávislé. Sdružená hustota se rozpadne:
$g_s(z_1, \dots, z_s) = \displaystyle\prod_{k=1}^{s} f_k(z_k)$
Stacionární tok: $P$ nezávisí na startu $t_0$, jen na délce $t$. Hustota intervalu $f$ je konstantní v čase. Intenzita toku:
$\lambda = \left(\displaystyle\int_0^\infty x\,f(x)\,dx\right)^{-1}$
Ordinární tok: $P(\ge 2$ příchodů v $\Delta t) \to 0$ rychleji než $\Delta t$. Tedy „dva najednou“ je zanedbatelné.

Poissonův tok: klasická aproximace

Intervaly mezi příchody mají exponenciální rozdělení $f(x) = \lambda e^{-\lambda x}$. Počet příchodů má Poissonovo rozdělení.

Stacionární tok (intenzita $\lambda$ konstantní):

$P(k, t) = \dfrac{(\lambda t)^k}{k!}\,e^{-\lambda t}$

Nestacionární tok (intenzita $\lambda(u)$ se mění, třeba špička v bance v 11:00):

$P(k, t_0, t) = \dfrac{\Lambda^k(t_0, t)}{k!}\,e^{-\Lambda(t_0,t)}, \qquad \Lambda(t_0,t) = \displaystyle\int_{t_0}^{t} \lambda(u)\,du$

Poissonův tok = „úplně náhodné příchody bez paměti“. Pravděpodobnost, že další zákazník přijde za 5 minut, nezávisí na tom, kdy přišel ten předchozí. Intervaly se generují přes inverzní transformaci: $\eta = -\frac{1}{\lambda}\ln \Upsilon$.

Dva přístupy k simulaci

1. Modely s diskrétním časem. Čas se rozseká na pevné kroky $\Delta t$. Na konci každého počítač zkontroluje stav systému (přišel požadavek? uvolnila se linka?).

Jednoduchá implementace.
Velká nevýhoda: pokud je $\Delta t$ malé, většina kroků je prázdných a marně se vyhodnocují.

2. Modely s diskrétními událostmi (event-driven). Čas se „přeskakuje“ rovnou na další skutečnou událost (příchod, dokončení obsluhy, porucha). Žádné prázdné kroky.

Řádově efektivnější.
Pro Monte Carlo simulace front standard.

Diskrétní čas = budík každou minutu kontroluje, jestli něco. Diskrétní události = budík zazvoní přesně ve chvíli, kdy se má něco stát. Druhé je téměř vždy lepší.

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

Co když chceme změřit pravděpodobnost zahlcení systému, ale ono nastane jednou za milion běhů? Naivní MC by potřeboval astronomický počet pokusů. Pomáhá splitting:

Princip

Sledování stavu: definujeme kontrolní okamžiky $T_1, T_2, \dots$ a kritickou hranici $N$ (např. počet požadavků v systému).
Rozštěpení: jakmile v okamžiku $T_i$ stav dosáhne $n \ge N$, simulaci v tomto bodě rozštěpíme na $M$ identických podsystémů. Každému přiřadíme váhu $Q = 1/M$.
Další větvení: v $T_{i+1}$ ověříme podmínku pro každý podsystém. Ty, co podmínku splní, znovu rozštěpíme, váhy se vynásobí ($1/M^2$, dále $1/M^3$, …). Větve, co nesplní, zahodíme.
Sloučení: výsledné rozdělení stavů spočítáme jako $p_r = Q_r / Q$, kde $Q_r$ je součet vah systémů ve stavu $r$, a $Q$ celkový součet.

Místo abychom marně čekali na vzácný jev, „klonujeme“ trajektorie, které už jsou na cestě k němu. Vážením je zajištěno, že odhad pravděpodobnosti zůstává nestranný, jen jsme dramaticky zvýšili počet pokusů v zajímavé oblasti.

Co se v simulaci vyhodnocuje

Průměrná délka fronty, max. délka fronty.
Průměrná doba čekání a obsluhy.
Procento zamítnutých / odchozích požadavků.
Vytížení linek (% času, kdy byly obsazené).
Pravděpodobnost vzácných stavů (přes splitting).

Kavárna se dvěma baristy. Intenzita příchodu $\lambda = 20$ zákazníků/hod (Poisson), doba obsluhy exponenciální se střední hodnotou 5 min. Otázka: stačí dva baristy v poledne, kdy $\lambda(u)$ vyskočí na 50?
Simulujeme event-driven: každá událost je „přišel zákazník“ nebo „barista skončil“. Sledujeme délku fronty, počet odcházejících. Spustíme 10 000 nezávislých běhů, srovnáme scénáře s 2, 3, 4 baristy.

Zapamatuj si rozdíl diskrétní čas vs. diskrétní události, kategorie podle čekání (odmítání / čekací / smíšený), a roli Poissonova rozdělení. Pro vzácné jevy je nástrojem splitting, větvení s vahami.

Materiál ke státnicím · Metody Monte Carlo · 5 otázek