Lekce 6 - Generování (pseudo)náhodných čísel

V minulé lekci, LU rozklad, vlastní čísla a definitnost matic, jsme se naučili dělat LU rozklad, ukázali si jeho aplikace a dozvěděli se, co jsou to vlastní čísla a definitnost matic.

Ať už jde o šifrování nebo generování úrovní ve hrách, potřebujeme v aplikacích často používat náhodná čísla nebo alespoň čísla, která se jako náhodná tváří. Čistě náhodná čísla klasický deterministický počítač generovat nedokáže, protože se řídí přesnými pravidly a neumí zkrátka jen tak "něco střelit". To lze v praxi obejít generováním tzv. pseudonáhodných čísel, která jsou sice počítána čistě předvídatelným způsobem, nicméně mají charakter náhodného šumu, což nám většinou postačuje. Pro případy, kdy jsou striktně vyžadována čistě náhodná čísla, mají některé procesory obvody, které taková čísla umí generovat na základě zákonů kvantové mechaniky. Tyto obvody ale nejsou příliš rychlé a běžně nejsou potřeba, proto se zde budeme zabývat čísly pseudonáhodnými a odteď kdykoli budu psát o náhodných číslech, myslím implicitně čísla pseudonáhodná.

Kongruentní generátor

Základem generování náhodných čísel je procedura generující hodnoty v určitém intervalu. Typicky se používá tzv. kongruentní generátor (existují alternativy jako třeba Mersenne twister). Kongruentní generátor funguje na základě rovnic:

x₀ = Z

x_{n + 1} = (A * x_n + C) mod M

Vidíme, že jde o vzorec pro posloupnost hodnot, kdy první hodnotu (x₀), tzv. seed, zvolíme (Z) a každou další spočítáme tak, že vezmeme aktuální hodnotu, vynásobíme ji zvolenou konstantou (A), přičteme jinou zvolenou konstantu (C) a s výsledkem provedeme operaci modulo M (zbytek po celočíselném dělení M). Konstanta M je tzv. perioda. Jde o důležité číslo, protože nejpozději po m krocích se generovaná posloupnost začne opakovat. Musíme tedy s tímto faktem počítat a volit M vhodně vysoké.

Při nevhodně zvolených parametrech můžeme u generovaných čísel vypozorovat určitou souvislost (např. když je zakreslíme jako body roviny, můžeme dostat přímku), proto existují různé statistické testy pro určování kvality generátoru. Vhodné hodnoty, které využívá např. standardní knihovna jazyka C jsou M = 2³², A = 1103515245, C = 12345.

Je zřejmé, že generátor bude vždycky vytvářet stejnou posloupnost v závislosti na počáteční hodnotě (seedu). Pro každé spuštění aplikace tedy dostaneme vždycky stejný výsledek, což se někdy, např. při simulacích, hodí a jindy, třeba ve hrách, ne. Pro dosažení různosti posloupností se proto typicky jako seed používá hodnota založená na aktuálním čase – jelikož program spouštíme v různých časech, generovaná posloupnost bude pokaždé jiná.

Jakmile máme číslo vygenerované, je vhodné jej převést do určitého intervalu, např. <0,1) (bez koncové jedničky, ta nemůže být vygenerována). To uděláme samozřejmě tak, že jej v plovoucí řádové čárce vydělíme hodnotou M.

Důvod, proč jedničku ve výstupu nechceme, vyplývá z požadavku použít tuto funkci pro generování celých čísel. Pokud chceme funkci použít pro generování celých čísel v intervalu <0,N>, tak jednoduše vynásobíme výstup generátoru (tedy číslo 0 až 1, bez koncové jedničky) hodnotou N + 1 a číslo zaokrouhlíme odseknutím desetinné části.

Pro úplnost uveďme příklad kongruentního generátoru v jazyce C:

#include <stdio.h>

unsigned int x = 10341; // nahodila hodnota - seed

double muj_rand()  // vraci nahodnou hodnotu v rozsahu <0,1>
  {
    x = (1103515245 * x + 12345) % 4294967296;
    return x / ((double) 4294967296); // prevedeni do intervalu <0,1)
  }

int main()
  {
    int i;
    for (i = 0; i < 100; i++)
      printf("%lf\n",muj_rand());

    getchar();
    return 0;
  }

Generování různých pravděpodobnostních rozložení

Generátor sám o sobě produkuje čísla v rovnoměrném rozložení, tzn. každé číslo z intervalu má stejnou pravděpodobnost být vygenerováno. To ovšem ne vždy potřebujeme. Někdy chceme generovat rozložení normální, exponenciální či nějaké vlastní. V takovém případě klasicky vygenerujeme hodnotu v běžném, rovnoměrném rozložení a použijeme některou z metod pro převod na jiná rozložení:

• inverzní transformace – Pro použití inverzní transformace potřebujeme znát distribuční funkci pravděpodobnosti daného rozložení, což je funkce, která říká, jaká je pravděpodobnost, že náhodná veličina nebude větší než daná hodnota, a která se značí F(x). Pokud k této funkci dokážeme sestrojit inverzní funkci F^-1(x), pak stačí vygenerovat náhodné číslo v intervalu <0,1> a předat jej jako argument této funkci. Číslo, které nám vyjde, je číslo v požadovaném rozložení. Dodejme, že inverzní funkce k funkci g(x) je funkce, které předáme číslo a ona nám řekne, jaký argument musíme dát funkci g(x), abychom z ní toto číslo dostali (zkrátka si osy x a y prohodí role).

Příklad inverzní transformace

• vylučovací metoda – Jde o poměrně přímočarou metodu. Pokud známe funkci rozložení hustoty pravděpodobnosti p(x), můžeme s použitím rovnoměrného rozložení generovat náhodně body ve dvourozměrném prostoru, tzn. souřadnice x a y, tak dlouho, dokud vygenerovaný bod nepadne pod graf funkce, což znamená, že číslo x je číslem vygenerovaným v daném rozložení. Hodnoty pro x samozřejmě generujeme v intervalu rozsahu funkce hustoty pravděpodobnosti a hodnoty y v intervalu oboru jejích hodnot, jinak bychom se výsledku nemuseli dočkat. Metoda není vhodná pro případy, kdy je velká pravděpodobnost, že vygenerovaný bod padne mimo, což znamená nové generování bodu a zpomalení algoritmu (takže např. funkce s definičním oborem všech reálných čísel jako je normální rozložení). Uveďme krátký příklad pro ilustraci:

double hustota_pravdepodobnosti(double x)

  // definicni obor funkce je zde od 0,2 po 1,5
  {
    ... // zde je definice funcke hustoty pravdepodobnosti
  }

double moje_rozlozeni()
  {
    double x, y;

    while (1)
      {
        x = muj_rand() * 1.3 + 0.2; // generuje x v rozsahu 0,2 az 1,5
        y = muj_rand() * 1.1; // generuje cislo v rozsahu 0 az 1,1

        if (y <= hustota_pravdepodobnosti(x))
          return x;
      }
  }

• centrální limitní věta pro normální rozdělení – Tato matematická věta říká, že součet několika hodnot v libovolném rozložení konverguje k rozložení normálnímu. Můžeme ji proto využít tehdy, když potřebujeme vygenerovat přibližně normální rozložení za pomoci rovnoměrného. Princip se dá demonstrovat na případu, kdy generujeme čísla 2 kostkami – je jenom jeden způsob, jak hodit číslo 2 (1 na obou kostkách), ale je mnoho způsobů, jak hodit číslo 7 (1 + 6, 2 + 5, 3 + 4 atd.), číslo 7 má tedy vyšší pravděpodobnost výskytu. Příkladem jednoduché funkce může být tato:

double normalni_rozlozeni()
  {
    int soucet, i;

    for (i = 0; i < 10; i++) // 10x rovnomerne rozlozeni
      soucet += muj_rand() * 11; // cislo 0 az 10

    return soucet / ((double) 100);
  }

• kombinace různých metod (pro složitější rozložení)

Závěrem

Vidíme, že generování náhodných čísel není triviálním problémem. Při mnohých vědeckých simulacích může být navíc tato akce kritická a proto bychom si měli být vědomi alespoň základních principů fungování náhodných generátorů, které by tato lekce měla pokrývat.

V další lekci, Algoritmus pro numerické integrování - Obdélníková metoda, si popíšeme algoritmus pro numerický výpočet integrálu, tedy výpočet obsahu pod křivkou. Tato část je zaměřena na Obdélníkovou formuli.