Metoda Monte Carlo – cvičení
Lineární kongruenční generátory
Nechť je prvočíslo, a je generátor multiplikativní grupy modulo . Definujme rekurentní posloupnost . Potom tato posloupnost bude postupně procházet všechna čísla z v pseudonáhodném pořadí. Uděláme-li z ní náhodnou proměnnou , potom bude neboli pro dostatečně velké máme . Ovšem zkusíme-li vyplnit čtverec tím, že střídavě generujeme souřadnice a , zjistíme, že ho to nevyplní, ale bude to na něm vytvářet jakési čáry.
Spočtěme si pro střední hodnotu a rozptyl:
Ze zákona velkých čísel víme, že definujeme-li pro dostatečně velké , bude přibližně . Když z toho chceme udělat standardní normální rozdělení, znormujeme to:
Ale pozor, zvolíme-li pro jednoduchost třeba , tahle náhodná proměnná ve skutečnosti vždy bude v intervalu . 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 a paraboloid daný nerovnicí . Budeme-li bod brát rovnoměrně náhodně, potom objem paraboloidu můžeme spočítat jako
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 . Poté můžeme odhadnout . Skutečně platí . Ale s jakým rozptylem?
Taky můžeme použít polární souřadnice. Zvolíme náhodně . 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 . Potom máme nerovnici a stačí spočítat .
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 toho si už můžeme snadno spočíst objem:
Nerovnoměrná rozdělení
Chceme generovat náhodnou proměnnou s exponenciální hustotou pravděpodobnosti
Odpovídající distribuční funkce je
Jak takovéto rozdělení nagenerujeme z náhodné proměnné s rovnoměrným rozdělením ? Stačí vyřešit rovnici . Vyjde nám
Tento přístup můžeme samozřejmě použít i pro jiná rozdělení. Mějme třeba Cauchyovo rozdělení
Opět budeme řešit rovnici , čímž dostaneme
Jako poslední příklad si zkusme tímto způsobem vygenerovat rovnoměrné rozdělení .
Řešením rovnice , přičemž uvažujeme jenom „prostou část“ uprostřed, dostaneme
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 , přičemž bereme . Tím si rozdělíme reálnou osu na šuplíčky. Pravděpodobnost, že se náhodná proměnná bude vyskytovat v -tém šuplíčku, je
Provedeme-li pozorování, dostaneme nějaká čísla vyjadřující, kolikrát jsme se trefili do kterého šuplíčku. Střední hodnota bude
V praxi chceme šuplíčky zvolit tak, aby bylo . Podle nějaké Pearsonovy věty platí
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:
Použijeme oblíbený trik, že si to rozšíříme do dvou rozměrů. Mějme funkci
Přejdeme-li do polárních souřadnic, máme . Pro distribuční funkci poloměru platí
Tuhle funkci už dokážeme invertovat: máme-li rovnoměrnou proměnnou , potom , tedy .
Poissonovo rozdělení
Máme exponenciální rozdělení ,
Dále existuje binomické rozdělení ,
Jeho limitou pro je Poissonovo rozdělení ,
Teď nás zajímá, jak Poissonovo rozdělení generovat. Mohli bychom opět pro náhodnou proměnnou řešit rovnici . Nyní už ale 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
Další možnost je aproximovat ho pomocí binomického, kde můžeme počítat
Výpočet integrálu
Chceme spočítat . Víme, li, že , můžeme vzít náhodné proměnné a máme
Střelíme-li náhodně -krát a z toho se -krát trefíme, můžeme aproximovat .
Pro tuto metodu je
Jiný způsob, jak by na to šlo jít, je vzít si náhodnou proměnnou a počítat
Pro tuto metodu je
Rozptyl je menší než u předchozí metody.
Nevlastní integrál
Zkusme spočítat integrál
Potřebujeme zavést nějakou substituci, která to dostane do omezeného intervalu. Jedna možnost je . Tím dostaneme
Další možnost by byla , čímž dostaneme
Nebo taky :
Vícerozměrné integrály
Ukážeme si, v čem metoda Monte Carlo vyniká oproti deterministickým metodám integrace. Chceme-li spočítat dělením intervalu na kousíčků velikosti a chceme přesnost , kde je řád konvergence metody, časová náročnost bude
Počítáme-li -rozměrný integrál, kde -tou dimenzi dělíme na kousíčků, časová náročnost bude
Kdybychom to na oblasti počítali metodou Monte Carlo, tak jak už jsme si odvozovali, přesnost bude
z čehož máme . Metoda Monte Carlo bude tedy rychlejší v případě, kdy neboli .