Metoda Monte Carlo (simulace v R-ku)

Autor: Lenka Fiřtová

V tomto článku si povíme něco málo o tom, jak provádět simulace.

Proč bychom vůbec chtěli provádět Monte Carlo simulace?

Monte Carlo simulace je simulace, která pracuje s náhodnými čísly (tedy spíše s čísly pseudonáhodnými: programy většinou vytvářejí “náhodná čísla” pomocí algoritmů tak, aby měla stejné vlastnosti jako čísla skutečně náhodná). Monte Carlo simulace je statistický experiment, v rámci kterého vygenerujeme velké množství (pseudo)náhodných čísel a pak z nich vyvozujeme závěry. Monte Carlo je vlastně řešením úlohy pomocí statistického experimentu.

Monte Carlo simulaci tedy využijeme všude, kde pracujeme s náhodou resp. s pravděpodobností. Například můžeme simulací zjistit, jaká je pravděpodobnost, že při hodu třemi kostkami padne jednou sudé číslo a dvakrát liché. Stačí si nasimulovat dostatek (čím více tím lépe, ideálně alespoň několik tisíc) takových hodů a podívat se, v kolika procentech případů nastala popsaná situace. Při znalosti “vzorečků” pravděpodobnosti lze zrovna tuto situaci snadno spočítat i bez simulace, ale některé situace obsahují výpočty natolik složité (či dokonce neproveditelné), že je lepší sáhnout po simulaci.

Příklad: odhad čísla π

Jedním z nejjednodušších příkladů použití Monte Carlo simulace je odhad čísla π. Jak na to? Představme si kruh vepsaný čtverci. Do tohoto čtverce náhodně “hážeme” body.

Uvažujme, že čtverec má stranu rovnou 2 a kruh jemu vepsaný má tedy poloměr roven 1 (r = 1). Pak obsah tohoto čtverce je roven 2r · 2r a obsah kruhu je roven π · r²  Podíl obsahu kruhu a obsahu čtverce je roven (π · r²)/(2r · 2r) = π / 4. 

Pokud do tohoto čtverce náhodně “hážeme” body, pak z nich určitý podíl spadne do kruhu (a zbytek mimo kruh). Jestliže body “hážeme” skutečně náhodně, pak podíl obsahu kruhu a obsahu čtverce se musí rovnat podílu bodů, které “spadly” do kruhu, a to vše se rovná π / 4 (viz předchozí odstavec). To znamená, že číslo π odpovídá čtyřnásobku podílu bodů, které spadly do kruhu.

Snáze si to představíme, podíváme-li se na následující obrázek:

Na obrázku je několik tisíc bodů: část z nich (černá) v kruhu, zbytek (modrá) mimo kruh. Řekli jsme, že čtverec má stranu o délce 2 (na jednotkách nezáleží). To znamená, že každý bod můžeme popsat dvěma souřadnicemi x, y například na škále od –1 do 1. Střed čtverce a současně střed kruhu má souřadnice [0; 0]. Jeden vygenrerovaný bod může mít například souřadnice [0,5; 0,75]. Leží takový bod uvnitř kruhu? To můžeme spočítat Pythagorovou větou. Vzdálenost tohoto bodu (označme c) od středu kruhu je totiž stejná jako délka přepony v trojúhelníku s odvěsnami 0,5 a 0,75. Takže platí: c² = 0,5² + 0,75², tzn. c = 0,901. To je méně než poloměr kruhu 1, takže tento bod leží v kruhu.

Použijeme R-ko a zkusíme vygenerovat 10 000 bodů: pro každý bod vybereme náhodně dvě čísla od –1 do 1 a to budou jeho souřadnice. Pythagorovou větou spočítáme, jestli bod leží uvnitř kruhu. Pak zjistíme, jaký podíl z 10 000 bodů ležel uvnitř kruhu.

Kód bude vypadat takto:

> r = 1  # poloměr kruhu a současně polovina strany čtverce

> k = 0  # do proměnné k si budeme započítávat body, které padly do kruhu: pokaždé, když bod padne do kruhu, navýšíme k o jedničku<

> for (i in 1:10000)   {   # 10 000krát zopakujeme následující postup:
           # z uniformního rozdělení funkcí runif vygenerujeme souřadnici x a souřadnici y
           x = runif(1,-1,1) # vygenerujeme jedno náhodné číslo z intervalu od –1 do 1
           y = runif(1,-1,1) # vygenerujeme jedno náhodné číslo z intervalu od –1 do 1
           if (sqrt(x*x + y*y) <= r)  { # použijeme Pythagorovu větu: 
           # pokud je vzdálenost bodu se souřadnicemi x, y od středu menší než poloměr r...
               k  = k + 1     # pak číslo k navýšíme o jedničku
               }
   }

Kolik tedy činí π? Jedná se o čtyřnásobek podílu bodů v kruhu (to je k) na všech bodech (to je 10 000).

Tedy π se rovná:

> 4*k/10000

Protože se jedná o simulaci, vyjde číslo pokaždé trochu jinak. Čím více bodů vygenerujeme, tím přesnější číslo bude. Mělo by se ale pohybovat okolo skutečné hodnoty π (3,14).

Vidíme, že jsme dokázali spočítat nějakou matematickou veličinu nikoli vzorečkem, nýbrž pomocí simulace.

Další příklad na simulace v R-ku viz zde.

Vysvětlení for cyklů viz zde.