Metoda Monte Carlo – cvičení ⬩ Adamátorovy zápisky

Lineární kongruenční generátory

Nechť $p$ je prvočíslo, $x_{0} \in \hat{p - 1}$ a $a$ je generátor multiplikativní grupy modulo $p$ . Definujme rekurentní posloupnost $x_{n} ≔ a \cdot x_{n - 1} mod p$ . Potom tato posloupnost bude postupně procházet všechna čísla z $\hat{p - 1}$ v pseudonáhodném pořadí. Uděláme-li z ní náhodnou proměnnou $X$ , potom bude $X \sim U (\hat{p - 1})$ neboli pro dostatečně velké $p$ máme $Y ≔ \frac{X}{p} \sim U ((0, 1))$ . Ovšem zkusíme-li vyplnit čtverec tím, že střídavě generujeme souřadnice $x$ a $y$ , zjistíme, že ho to nevyplní, ale bude to na něm vytvářet jakési čáry.

Spočtěme si pro $Y$ střední hodnotu a rozptyl:

E Y = \int_{0}^{1} x f_{Y} (x) d x = \int_{0}^{1} x d x = \frac{1}{2},

Var Y = E (X - \frac{1}{2}) = \frac{1}{12} .

Ze zákona velkých čísel víme, že definujeme-li $Z ≔ \sum_{i = 1}^{m} Y_{i}$ pro dostatečně velké $m$ , bude přibližně $Z \sim N (\frac{m}{2}, \frac{m}{12})$ . Když z toho chceme udělat standardní normální rozdělení, znormujeme to:

Q ≔ \frac{Z - \frac{m}{2}}{\sqrt{\frac{m}{12}}} .

Ale pozor, zvolíme-li pro jednoduchost třeba $m = 12$ , tahle náhodná proměnná ve skutečnosti vždy bude v intervalu $(- 6, 6)$ . Tedy nemáme skutečně normální rozdělení.

Výpočet objemu

Pojďme si metodou Monte Carlo spočítat objem paraboloidu. Mějme bod $u = (x, y, z) \in {[- R, R]}^{2} \times [0, H]$ a paraboloid $A$ daný nerovnicí $z \geq \frac{H}{R^{2}} (x^{2} + y^{2})$ . Budeme-li bod brát rovnoměrně náhodně, potom objem paraboloidu můžeme spočítat jako

V = 4 R^{2} H \cdot p [u \in A] ≕ p .

Samozřejmě abychom tuto pravděpodobnost spočetli, musíme hodněkrát náhodně střílet. Budeme tedy mít náhodnou proměnnou $M \sim Bi (n, p)$ . Poté můžeme odhadnout $p \approx \frac{M}{N}$ . Skutečně platí $E M = N \cdot p$ . Ale s jakým rozptylem?

σ^{2} = N \cdot p \cdot (1 - p) \approx N \cdot \frac{M}{N} \cdot (1 - \frac{M}{N}) = M \cdot (1 - \frac{M}{N})

σ = \sqrt{M} \cdot \sqrt{1 - \frac{M}{N}}

γ = \frac{σ}{E M} = \frac{\sqrt{\frac{1 - M}{N}}}{\sqrt{M}} < \frac{1}{\sqrt{M}} = 𝒪 (M^{- \frac{1}{2}}) = 𝒪 (N^{- \frac{1}{2}})

Taky můžeme použít polární souřadnice. Zvolíme náhodně $φ \sim U (0, 2 π)$ . Ovšem u poloměru už nechceme rovnoměrné rozdělení, protože bodů vzdálenějších od středu je víc. Na přednášce si odvodíme, že je vhodné zvolit $r \sim \sqrt{U (0, 1)}$ . Potom máme nerovnici $z \geq \frac{H}{R^{2}} \cdot r$ a stačí spočítat $π R^{2} H \cdot p$ .

Ještě jiný přístup. Místo abychom náhodně stříleli body v trojrozměrném prostoru, je budeme střílet ve dvojrozměrném a pro každý spočteme výšku paraboloidu v daném bodě. Tuto výšku spočteme jako

z = (1 - \frac{r^{2}}{R^{2}}) \cdot H = (1 - \frac{x^{2} + y^{2}}{R^{2}}) \cdot H ≕ g (x, y),

E g (X, Y) \approx \frac{1}{N} \sum_{i = 1}^{N} g (x_{i}, y_{i}) ≕ \bar{g},

s_{x}^{2} = \frac{1}{N - 1} \sum_{i = 1}^{N} {(g (x_{i}, y_{i}) - \bar{g})}^{2} .

Z toho si už můžeme snadno spočíst objem:

V = π r^{2} \cdot E g (X, Y) .

Nerovnoměrná rozdělení

Chceme generovat náhodnou proměnnou s exponenciální hustotou pravděpodobnosti

f (x) = [x \geq 0] \cdot \frac{1}{μ} \cdot exp (- \frac{x}{μ}) .

Odpovídající distribuční funkce je

F (x) = \int_{- \infty}^{x} f = [x \geq 0] \cdot (1 - exp (- \frac{x}{μ})) .

Jak takovéto rozdělení nagenerujeme z náhodné proměnné s rovnoměrným rozdělením

r \sim R (0, 1)

? Stačí vyřešit rovnici

F (x) = r

. Vyjde nám

x = - μ \cdot ln (1 - r) = - μ \cdot ln r^{*}, r^{*} \sim R (0, 1) .

Tento přístup můžeme samozřejmě použít i pro jiná rozdělení. Mějme třeba Cauchyovo rozdělení

f (x) = \frac{1}{π (1 + x^{2})},

F (x) = \frac{1}{2} + \frac{arctan x}{π} .

Opět budeme řešit rovnici

F (x) = r

, čímž dostaneme

x = tan (\frac{π}{2} \cdot (2 r - 1)) .

Jako poslední příklad si zkusme tímto způsobem vygenerovat rovnoměrné rozdělení

x \sim R (a, b)

f (x) = \frac{[x \in [a, b]]}{b - a},

F (x) = {\begin{cases} 0 & x < a, \\ \frac{x - a}{b - a} & x \in [a, b], \\ 1 & x > b . \end{cases}

Řešením rovnice

F (x) = r

, přičemž uvažujeme jenom „prostou část“ uprostřed, dostaneme

x = a + (b - a) \cdot r .

To jsme samozřejmě mohli uhodnout, ale je dobré si ověřit, že systematický přístup sedí s intuicí.

Testování náhodné volby

Mějme nějaké náhodné rozdělení. Zvolíme si body

x_{1} < x_{2} < \dots < x_{n - 1}

, přičemž bereme

x_{0} ≔ - \infty, x_{n} ≔ \infty

. Tím si rozdělíme reálnou osu na šuplíčky. Pravděpodobnost, že se náhodná proměnná bude vyskytovat v

k

-tém šuplíčku, je

p_{k} ≔ \int_{x_{k - 1}}^{x_{k}} f = F (x_{k}) - F (x_{k + 1}) .

Provedeme-li

N

pozorování, dostaneme nějaká čísla

N_{1}, \dots, N_{n}

vyjadřující, kolikrát jsme se trefili do kterého šuplíčku. Střední hodnota bude

O_{k} = N \cdot p_{k} .

V praxi chceme šuplíčky zvolit tak, aby bylo

O_{k} \geq 5

. Podle nějaké Pearsonovy věty platí

χ^{2} ≔ \sum_{k = 1}^{n} \frac{{(N_{k} - O_{k})}^{2}}{O_{k}} \sim χ_{n - 1}^{2} .

Box-Müllerova transformace

Cílem je vygenerovat z rovnoměrného rozdělení normální rozdělení. Kdybychom chtěli použít obecnou metodu pro nerovnoměrná rozdělení, dopadli bychom blbě, protože jak víme, kumulativní distribuční funkce nemá elementární tvar:

F (x) = \int_{- \infty}^{x} f = \int_{- \infty}^{x} \frac{1}{\sqrt{2 π}} \cdot exp (- \frac{t^{2}}{2}) d t .

Použijeme oblíbený trik, že si to rozšíříme do dvou rozměrů. Mějme funkci

f_{2} (x, y) ≔ \frac{1}{\sqrt{2 π}} \cdot exp (- \frac{x^{2} + y^{2}}{2}) .

Přejdeme-li do polárních souřadnic, máme

φ \sim U ((o, 2 π))

. Pro distribuční funkci poloměru platí

F_{R} (r) = \iint_{x^{2} + y^{2} \leq r^{2}} \frac{1}{2 π} \cdot exp (- \frac{x^{2} + y^{2}}{2}) d (x, y) = \int_{0}^{r} \int_{0}^{2 π} \frac{1}{2 π} \cdot exp (- \frac{{\tilde{r}}^{2}}{2}) d φ d \tilde{r} = \dots = 1 - exp (- \frac{r^{2}}{2}) .

Tuhle funkci už dokážeme invertovat: máme-li rovnoměrnou proměnnou

u

, potom

F_{R} (r) = u

, tedy

r = \sqrt{- 2 \cdot ln (1 - u)} = \sqrt{- 2 \cdot ln u^{*}}

Poissonovo rozdělení

Máme exponenciální rozdělení

X \sim Ex (μ)

f (x) = \frac{1}{μ} \cdot exp (- \frac{x}{μ}) \cdot I (x \geq 0),

F (x) = (1 - exp (- \frac{x}{μ})) \cdot I (x \geq 0) .

Dále existuje binomické rozdělení

Y \sim Bi (N, p)

p_{y} (N, p) = (\binom{N}{y}) p^{y} {(1 - p)}^{N - y} .

Jeho limitou pro

N \to \infty, N \cdot p = λ

je Poissonovo rozdělení

X \sim Po (λ)

p_{x} (λ) = lim_{N \to \infty} p_{y} (N, \frac{λ}{N}) = \frac{λ^{x}}{x!} \cdot exp (- λ) .

Teď nás zajímá, jak Poissonovo rozdělení generovat. Mohli bychom opět pro náhodnou proměnnou

U \sim U (0, 1)

řešit rovnici

F (x) = u

. Nyní už ale

F

není prostá funkce, takže si musíme ze všech možných řešení zvolit to nejmenší. Pro zjednodušení můžeme počítat

F_{0} = p_{0} = exp (- λ), F_{x} = F_{x - 1} + p_{x} = F_{x - 1} + \frac{λ}{x} \cdot p_{x - 1} .

Další možnost je aproximovat ho pomocí binomického, kde můžeme počítat

F_{0} = p_{0} = {(1 - p)}^{N}, F_{y} = F_{y - 1} + p_{y} = F_{y - 1} + \frac{p \cdot (N - y + 1)}{(1 - p) \cdot y} \cdot p_{y - 1} .

Výpočet integrálu

Chceme spočítat

I ≔ \int_{a}^{b} f (x) d x

. Víme-li, že

Z \leq f (x) \leq Q

, můžeme vzít náhodné proměnné

X \sim U (a, b), Y \sim U (Z, Q)

a máme

P [Y \leq f (X)] = \frac{\int_{a}^{b} (f (x) - z) d x}{(b - a) \cdot (Q - z)} ≕ p_{HIT},

I = (b - a) \cdot (Z + (Q - Z) \cdot p_{HIT}) .

Střelíme-li náhodně

N

-krát a z toho se

K

-krát trefíme, můžeme aproximovat

p_{HIT} \approx \frac{K}{N}

.Pro tuto metodu je

σ = (Q - Z) \cdot \sqrt{\frac{p_{HIT} \cdot (1 - p_{HIT})}{N}} = 𝒪 (\frac{1}{\sqrt{N}})

Jiný způsob, jak by na to šlo jít, je vzít si náhodnou proměnnou

X \sim U (a, b)

a počítat

E f (X) = \frac{1}{b - a} \cdot \int_{a}^{b} f (x) d x = \frac{I}{b - a},

I = (b - a) \cdot E f (X) \approx \frac{b - a}{N} \cdot \sum_{i = 1}^{N} f (X_{i}), x_{i} \sim U (a, b) .

Pro tuto metodu je

σ = \frac{1}{\sqrt{N}} \cdot \sqrt{(b - a) \cdot \int_{a}^{b} {g (x)}^{2} d x - {(\int_{a}^{b} g (x) d x)}^{2}} = 𝒪 (\sqrt{N}) .

Rozptyl je menší než u předchozí metody.

Nevlastní integrál

Zkusme spočítat integrál

\int_{0}^{\infty} \frac{cos ω x}{1 + x^{2}} d x .

Potřebujeme zavést nějakou substituci, která to dostane do omezeného intervalu. Jedna možnost je

ξ ≔ \frac{x}{1 + x}

. Tím dostaneme

\int_{0}^{1} \frac{cos \frac{ω ξ}{1 - ξ}}{{(1 - ξ)}^{2} + ξ^{2}} d ξ .

Další možnost by byla

ξ ≔ exp (- x)

, čímž dostaneme

\int_{0}^{1} \frac{cos (ω \cdot ln ξ)}{(1 + {ln}^{2} ξ) \cdot ξ} d ξ .

Nebo taky

ξ ≔ \frac{2}{π} \cdot arctan x

\int_{0}^{1} \frac{π \cdot cos (ω \cdot tan \frac{π ξ}{2})}{(1 + {tan}^{2} \frac{π ξ}{2}) \cdot 2 \cdot {cos}^{2} \frac{π ξ}{2}} .

Vícerozměrné integrály

Ukážeme si, v čem metoda Monte Carlo vyniká oproti deterministickým metodám integrace. Chceme-li spočítat

\int_{a}^{b} f

dělením intervalu na

N

kousíčků velikosti

h ≔ \frac{a - b}{N}

a chceme přesnost

ε = 𝒪 (h^{k})

, kde

k

je řád konvergence metody, časová náročnost bude

T (N) = 𝒪 (N) = 𝒪 (h^{- 1}) = 𝒪 (ε^{- \frac{1}{k}}) .

Počítáme-li

d

-rozměrný integrál, kde

j

-tou dimenzi dělíme na

N_{j}

kousíčků, časová náročnost bude

T (N) = 𝒪 (\prod_{j = 1}^{d} N_{j}) = 𝒪 (\prod_{j = 1}^{d} h_{j}^{- 1}) = 𝒪 (h^{- d}) = 𝒪 (ε^{- \frac{d}{k}}) .

Kdybychom to na oblasti

O \subset ℝ^{d}

počítali metodou Monte Carlo, tak jak už jsme si odvozovali, přesnost bude

ε = \frac{σ \cdot V (O)}{\sqrt{T (N)}},

z čehož máme

T (N) = 𝒪 (ε^{- 2})

. Metoda Monte Carlo bude tedy rychlejší v případě, kdy

\frac{d}{k} > 2

neboli

d > 2 k