Monte Carlo integrování
![](http://upload.wikimedia.org/wikipedia/commons/thumb/b/b0/MonteCarloIntegrationCircle.png/220px-MonteCarloIntegrationCircle.png)
V matematice je Monte Carlo integrování postup numerického odhadu hodnoty integrálu funkce pomocí náhodného vzorkování. Jedná se o speciální případ Monte Carlo metody. Zatímco jiné algoritmy obvykle vyhodnocují integrand v pravidelné mřížce,[1] tak Monte Carlo algoritmus volí tyto body náhodně.[2] Tato metoda poskytuje lepší výsledky v prostorech s vyšší dimenzí, než klasické kvadraturní vzorce.[3]
Úloha
Budeme odhadovat hodnotu integrálu funkce definovanou jako
Estimátory
Primární estimátor
Primární estimátor integrálu definujeme vzorcem
kde je libovolná náhodná veličina s hustotou pravděpodobnosti , pro kterou platí
Nestrannost primárního estimátoru
Nestrannost primárního estimátoru plyne z jeho definice
To znamená, že střední hodnota primárního estimátoru je hodnota, kterou odhadujeme.
Rozptyl primárního estimátoru
Rozptyl se používá jako měřítko pro určení chyby estimátoru.
Pokud bude podobná bude výsledný rozptyl malý.
Sekundární estimátor
Jelikož použití pouze jednoho vzorku v primárním estimátoru nezaručuje dostatečně malý rozptyl odhadu, používá se estimátor sekundární. Ten využívá nezávislých náhodných veličin a jako výsledek se vezme jejich průměr, tedy:
Tento estimátor je opět nestranný, jelikož platí:
Z odvození nestrannosti vidíme, že hustota pravděpodobnosti obecně nemusí být identická pro všechny vzorky.
Rozptyl sekundárního estimátoru
Pro rozptyl sekundárního estimátoru platí vztah
Rozptyl je tedy přímo úměrný rozptylu primárního estimátoru s koeficientem . Standardní odchylka je definována jako odmocnina z rozptylu, a proto je pouze -krát menší.
Mějme funkci pro , tedy předpis polokoule s poloměrem . Budeme se snažit odhadnout hodnotu integrálu této funkce na množině . Rozšíříme definiční obor funkce pro . Takže
Budeme předpokládat, že hustota pravděpodobnosti je konstantní (jedná se o uniformní rozdělení). Toto je jednoduchá ukázka kódu napsaném v programovacím jazyce C++, který by daný problém mohl řešit.
#include<math.h> #include<random> double rand01() { return (double)rand() / RAND_MAX; } const double R = 1; double f(double x, double y) { double fx = R * R - x * x - y * y; return fx < 0 ? 0 : sqrt(fx); } void integral(int n, double &I) { double r, x; r = 0; for(int i= 0; i < n; i++) { x = f(R * (2 * rand01() - 1), R * (2 * rand01() - 1)); r += x; } // 1/(4 * R * R) je pravdepodobnost vyberu vzorku z uniformniho rozdeleni. // Timto prvkem muzeme vysledek delit az nakonci pouze diky vyuziti uniformniho rozdeleni, // nebot kazdy prvek ma tu samou pravdepodobnost vyskytu. I = 4 * R * R * r / n; }
Metody pro snížení rozptylu odhadu
Vzorkování podle důležitosti (Importance sampling)
Části vzorkovaného intervalu, kde má větší hodnotu, jsou důležitější, protože vzorky z těchto oblastí více ovlivňují výsledek. Vzorkování podle důležitosti umisťuje vzorky přednostně do takových oblastí. Toho se docílí tím, že pdf, ze které se vybírají vzorky, bude podobná integrandu. Použitím vzorkování podle důležitosti lze snížit rozptyl při zachování nestrannosti.
Pokud bychom použili pdf přímo úměrnou , tedy bychom měli
kde je normalizační faktor definovaný jako
který zajišťuje, že integrál přes celou doménu je a tedy že je hustota pravděpodobnosti. Rozptyl odhadu by pak byl nulový, jak je vidět z následujících úprav.
Bohužel v praxi se taková pdf nedá použít, protože pro její konstrukci bychom potřebovali znát hodnotu integrálu, který se snažíme spočítat.
Odhad integrálu pomocí řídící funkce (Control Variate)
Jednou z možností jak snížit rozptyl odhadu je využití takzvané řídící funkce . Důležité je aby tato funkce byla analyticky integrovatelná a zároveň aproximovala námi integrovanou funkci .
Pravý člen posledního výrazu umíme zintegrovat analyticky a levý člen integrujeme numericky jako doposud. Výhodou je to, že tím jak funkce aproximuje funkci , tak jejich rozdíl bude mít menší rozptyl. Tato metoda je lépe použitelná v případě, kdy se analyticky integrovatelná funkce vyskytuje v integrované funkci jako aditivní člen.
Vzorkování po částech (Stratified Sampling)
Při hledání vzorků pro odhad integrálu dochází v popsaných metodách často k náhodnému shlukování prvků. Tyto shluky silně přispívají k velikosti rozptylu. Popíšeme dvě metody snažící se o potlačení těchto shluků. První z nich je právě vzorkování po částech a jak název napovídá, tak prostor, přes který integrujeme, rozdělíme do disjunktních částí a odhad integrálu budeme počítat na těchto podmnožinách. Výsledkem bude součet odhadů integrálů na těchto podmnožinách. Tento přístup potlačuje shlukování vzorků a dá se ukázat, že rozptyl takového odhadu bude menší nebo roven rozptylu sekundárního estimátoru se stejným počtem vzorků. Tento postup je velmi účinný pro nízkou dimenzi prostoru, v kterém integrujeme. Nejjednodušší možností jak prostor rozdělit je uniformní rozklad. Lepších výsledků lze ale dosáhnout, pokud budeme dělení volit tak, aby rozptyl na jednotlivých částech prostoru byl co nejmenší.
Metody Quasi Monte Carlo
Místo náhodného vzorkování se používají čistě deterministické metody volby vzorků. Prezentované výsledky platí i pro QMC metody, ale v důkazech nelze využít statistiky. Determinismem se snažíme odstranit některé vlastnosti náhodných vzorků, které nám vadily a to především shlukování. Z toho důvodu zavádíme pojem diskrepance.
Definice diskrepance
Nejdříve mějme funkci
kde . Objem kvádru definovaného funkcí v je
Víme, že Monte Carlo odhad objemu tohoto kvádru je roven
kde je počet vzorků, jsou jednotlivé vzorky a značí počet vzorků, které spadly kvádru . Diskrepance množiny bodů je pak definovaná jako maximální možná chyba odhadu objemu přes všechny možné kvádry v daném prostoru. Tedy diskrepance je
Využití diskrepance
Diskrepance slouží jako míra uniformity dané množiny. Konverguje k nule pro .
Existuje teoreticky (Koksma–Hlawka nerovnost) odhad horní hranice chyby MC odhadu integrálu funkce, který je omezený součinem právě diskrepance a variací integrované funkce. To je důvod proč se snažíme najít vzorkovací sekvence s co nejnižší diskrepancí. Jedna ze sekvencí s nízkou diskrepancí je Van der Corputova a její rozšíření do prostorů s vyšší dimenzí od Haltona nebo Hammersleyho.