Cesko.wiki Monte Carlo integrování
Author
Albert FloresIlustrace Monte Carlo integrace. Použití náhodných vzorků pro určení toho jak velkou část čtverce zabírá kruh 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, tak Monte Carlo algoritmus volí tyto body náhodně. Tato metoda poskytuje lepší výsledky v prostorech s vyšší dimenzí, než klasické kvadraturní vzorce.
Úloha
Budeme odhadovat hodnotu integrálu I funkce f:\Omega\rightarrow \mathbb{R} definovanou jako :I=\int\limits_{\Omega} f(x)\, \mathrm{d}x.
Estimátory
Primární estimátor
Primární estimátor F_{\mathrm{prim}} integrálu I definujeme vzorcem :F_{\mathrm{prim}}=\frac{f(X)}{p(X)}, kde X je libovolná náhodná veličina s hustotou pravděpodobnosti p(x), pro kterou platí :\forall x \in \Omega: f(x) \neq 0 \implies p(x) \neq 0.
Nestrannost primárního estimátoru
Nestrannost primárního estimátoru plyne z jeho definice :E[F_{\mathrm{prim}}]=\int\limits_\Omega \frac{f(x)}{p(x)}p(x)\,\mathrm{d}x=I. 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. :V[F_{\mathrm{prim}}]=E[F_{\mathrm{prim}}^2]-E[F_{\mathrm{prim}}]^2=\int\limits_\Omega\frac{f(x)^2}{p(x)}\,\mathrm{d}x-I^2 Pokud p(x) bude podobná f(x) 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á N nezávislých náhodných veličin Y_i=f(X_i)/p(X_i) a jako výsledek se vezme jejich průměr, tedy: :F_N=\frac{1}{N}\sum\limits_{i=1}^{N} \frac{f(X_i)}{p(X_i)}. +more Tento estimátor je opět nestranný, jelikož platí: :E[F_N]=E\left[\frac{1}{N}\sum{Y_i}\right]=\frac{1}{N}\sum{E[Y_N]}=\frac{1}{N}NI=I. 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 :V[F_N]=V\left[\frac{1}{N} \sum_i{Y_i}\right ] = \frac{1}{N^2} N V[Y_i] = \frac{1}{N}V[F_{\mathrm{prim}}]. Rozptyl je tedy přímo úměrný rozptylu primárního estimátoru s koeficientem 1/N. +more Standardní odchylka je definována jako odmocnina z rozptylu, a proto je pouze \sqrt{N}-krát menší.
Mějme funkci f(x,y) = \sqrt{R^2-x^2-y^2} pro x^2+y^2 \leq R^2, tedy předpis polokoule s poloměrem R. Budeme se snažit odhadnout hodnotu integrálu této funkce na množině \{(x,y),|x|. +more Rozšíříme definiční obor funkce pro x,y\in\mathbb{R}. Takže :\tilde{f}(x,y)=\begin{cases} 0 & x^2+y^2 > R^2 \\ f(x,y)& x^2+y^2 \leq R^2 \end{cases}.
Budeme předpokládat, že hustota pravděpodobnosti p(x) 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 #include
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
Metody pro snížení rozptylu odhadu
Vzorkování podle důležitosti (Importance sampling)
Části vzorkovaného intervalu, kde má f 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í. +more 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 f(x), tedy bychom měli :p(x)=\frac{f(x)}{N}, kde N je normalizační faktor definovaný jako :N=\int\limits_\Omega f(x)\, \mathrm{d}x, který zajišťuje, že integrál p(x) přes celou doménu je 1 a tedy že p(x) je hustota pravděpodobnosti. Rozptyl odhadu by pak byl nulový, jak je vidět z následujících úprav. +more : \begin{align} &V[F_{\mathrm{prim}}]=\int\limits_\Omega\frac{f(x)^2}{p(x)}\,\mathrm{d}x-I^2=\int\limits_\Omega f(x)\cdot N\,\mathrm{d}x-I^2=& \\ &=N\int\limits_\Omega f(x) \,\mathrm{d}x-I^2 = \left(\int\limits_\Omega f(x) \,\mathrm{d}x\right)^2-I^2=0& \end{align}.
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 g. Důležité je aby tato funkce byla analyticky integrovatelná a zároveň aproximovala námi integrovanou funkci f. +more :I=\int_{\Omega}f(x) \mathrm{d}x = \int_{\Omega}\left(f(x)-g(x)\right)\mathrm{d}x + \int_{\Omega}g(x) \mathrm{d}x 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 g aproximuje funkci f, 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. +more 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. +more 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 : L(z) = \begin{cases} 1, & \text{pro }0\leq z_1, \dots, 0\leq z_s \\ 0, & \text{jinak,} \end{cases}
kde z\in \mathbb{R}^s. Objem kvádru definovaného funkcí L v A\subseteq \mathbb{R}^s je :V(A)=\prod_{j=1}^s{v_j}. +more Víme, že Monte Carlo odhad objemu tohoto kvádru je roven :\frac{1}{N}\sum_{i=1}^N{L(z_i)}=\frac{m(A)}{N}, kde N je počet vzorků, z_i jsou jednotlivé vzorky a m(A) značí počet vzorků, které spadly kvádru L. 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 :D^*(z_1,\dots,z_N)=\sup_A{\left|\frac{m(A)}{N}-V(A)\right|}.
Využití diskrepance
Diskrepance slouží jako míra uniformity dané množiny. Konverguje k nule pro N\rightarrow \infty.
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í. +more Jedna ze sekvencí s nízkou diskrepancí je Van der Corputova a její rozšíření do prostorů s vyšší dimenzí od Haltona nebo Hammersleyho.