Simulace v jazyce R

Autor článku: Lenka Fiřtová

V tomto článku je popsán způsob, jak provést v R simulaci. Simulaci provádíme obvykle ve chvíli, kdy je příliš obtížné spočítat to, co nás zajímá – třeba průměrnou hodnotu nějaké veličiny ve složitější situaci nebo pravděpodobnost nějakého jevu.

Příklad 1

Zjistěte simulací, jaká je pravděpodobnost, že na kostce padne šestka.

Nejprve úvod do pravděpodobnosti: pravděpodobnost jevu (například jevu „na kostce padla šestka“) zjistíme tak, že vydělíme počet příznivých výsledků (ten je jeden: „padla šestka“) počtem všech možných výsledků (těch je šest: „padlo některé z čísel 1 až 6“). Například při hodu šestistěnnou hrací kostkou je tedy pravděpodobnost, že padne šestka, rovna 1/6.

Kdybychom to nevěděli a neuměli to spočítat, pak bychom mohli kostkou házet pořád dokola a postupně bychom viděli, že šestka padá v jedné šestině případů. Z toho bychom odvodili, že pravděpodobnost, že padne šestka, je rovna 1/6.

Abychom určili přesnou pravděpodobnost, museli bychom kostkou hodit nekonečněkrát. To samozřejmě není možné, hodíme-li ale kostkou třeba tisíckrát, získáme už poměrně přesný odhad pravděpodobnosti (čím vícekrát kostkou hodíme, tím přesnější tento odhad bude).

Samozřejmě doopravdy nebudeme házet dokola kostkou, ale hody si pouze nasimulujeme.

Nasimulujeme například tisíc hodů. Tisíckrát tedy náhodně vybereme jedno číslo od 1 do 6. Pak spočítáme, v kolika z této tisícovky hodů padla šestka.

Kód pro R může vypadat například takto:

hod = sample(c(1:6), 1000, replace = TRUE)
length(which(hod == 6))
[1] 0.162

První řádek kódu říká, že nasimulujeme 1000 hodů kostkou: z čísel 1 až 6 (1:6) jich vybereme tisíc s vracením (s opakováním). Argument replace musíme nastavit na TRUE, jinak bychom vybírali bez opakování. Výsledek uložíme do proměnné hod, což bude tedy řada tisíce číslic 1 až 6.

Teď se podívejme na druhý řádek. Příkaz which(hod == 6) nám řekne, které prvky vektoru hod jsou rovny 6 (tedy ve kterých hodech padla šestka). Příkaz length vrací délku vektoru. Řekne tedy, kolik bylo celkem takových hodů, ve kterých padla šestka.

V našem případě šlo o 16,2 % hodů. Správně by měla šestka padnout v 1/6 hodů, což se rovná 16,67 %. Vidíme, že simulace nám umožnila určit pravděpodobnost padnutí šestky poměrně přesně i s „pouhou“ tisícovkou hodů.

Příklad 2

V článku o while cyklech je následující příklad: uvažujte skupinu, kde je 45 mužů a 5 žen. Vybírejte postupně náhodně jednu osobu. Na kolikátý pokus jste vybrali ženu? (tutéž osobu můžete vybrat i opakovaně, tzn. vybranou osobu vracíte zpět do skupiny)

Tento příklad je řešen následovně, pro vysvětlení viz daný článek.

skupina = c(rep(0, 45), rep(1, 5))
kolikrat = 0
vybrany = 0
while(vybrany == 0) {
  kolikrat = kolikrat + 1
  vybrany = sample(skupina, 1)
}

Při jednom konkrétním pokusu můžeme vybrat ženu hned napoprvé, ale třeba i až na posté. Nyní nás tedy bude zajímat, na kolikátý pokus vybereme ženu v průměru. Použijeme opět simulaci. Celý výše uvedený kód spustíme tisíckrát a pokaždé si uložíme, na kolikátý pokus jsme vybrali ženu. Tyto hodnoty pak zprůměrujeme.

Založíme si novou proměnnou výsledek. Tam v každém běhu simulace uložíme, na kolikátý pokus jsme vybrali ženu.

vysledek = rep(NA, 1000)

Zbytek kódu bude vypadat následovně.

for(i in 1:1000) {
  kolikrat = 0
  vybrany = 0
  while(vybrany == 0){
    kolikrat = kolikrat + 1
    vybrany = sample(skupina, 1)
  }
 vysledek[i] = kolikrat
}

Vidíme, že jde o tentýž kód jako výše, pouze „zabalený“ do for cyklu, který se používá při simulacích velice často (pro for cykly viz článek).

V každém z tisíce běhů simulace (tedy v každé tzv. replikaci) pro i od 1 do 1000 proběhne kód uvedený výše. Výsledek se pak uloží na i-té místo proměnné výsledek. Na pátém místě proměnné výsledek tedy například bude, na kolikátý pokus jsme vybrali ženu v pátém běhu simulace (v páté replikaci).

Pak už jen stačí spočítat průměrnou hodnotu proměnné výsledek.

> mean(vysledek)
[1] 9.969

Jelikož jde o náhodnou veličinu (provádíme náhodné výběry ze skupiny), bude se výsledné číslo vždy mírně lišit, ale mělo by se pohybovat okolo 10. Ženu tedy v průměru vybereme na desátý pokus.

Číslo, které by mohlo být poměrně obtížné spočítat, jsme simulací zjistili velmi snadno.