Monte Carlo integrování

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]

Budeme odhadovat hodnotu integrálu ${\displaystyle I}$ funkce ${\displaystyle f:\Omega \rightarrow \mathbb {R} }$ definovanou jako

${\displaystyle I=\int \limits _{\Omega }f(x)\,\mathrm {d} x.}$

Primární estimátor ${\displaystyle F_{\mathrm {prim} }}$ integrálu ${\displaystyle I}$ definujeme vzorcem

${\displaystyle F_{\mathrm {prim} }={\frac {f(X)}{p(X)}},}$

kde ${\displaystyle X}$ je libovolná náhodná veličina s hustotou pravděpodobnosti ${\displaystyle p(x)}$, pro kterou platí

${\displaystyle \forall x\in \Omega :f(x)\neq 0\implies p(x)\neq 0.}$

Nestrannost primárního estimátoru plyne z jeho definice

${\displaystyle 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 se používá jako měřítko pro určení chyby estimátoru.

${\displaystyle 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 ${\displaystyle p(x)}$ bude podobná ${\displaystyle f(x)}$ bude výsledný rozptyl malý.

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á ${\displaystyle N}$ nezávislých náhodných veličin ${\displaystyle Y_{i}=f(X_{i})/p(X_{i})}$ a jako výsledek se vezme jejich průměr, tedy:

${\displaystyle F_{N}={\frac {1}{N}}\sum \limits _{i=1}^{N}{\frac {f(X_{i})}{p(X_{i})}}.}$

Tento estimátor je opět nestranný, jelikož platí:

${\displaystyle 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.

Pro rozptyl sekundárního estimátoru platí vztah

${\displaystyle V[F_{N}]=V\left[{\frac {1}{N}}\sum _{i}{Y_{i}}\right]={\frac {1}{N^{2}}}NV[Y_{i}]={\frac {1}{N}}V[F_{\mathrm {prim} }].}$

Rozptyl je tedy přímo úměrný rozptylu primárního estimátoru s koeficientem ${\displaystyle 1/N}$. Standardní odchylka je definována jako odmocnina z rozptylu, a proto je pouze ${\displaystyle {\sqrt {N}}}$-krát menší.

Mějme funkci ${\displaystyle f(x,y)={\sqrt {R^{2}-x^{2}-y^{2}}}}$ pro ${\displaystyle x^{2}+y^{2}\leq R^{2}}$, tedy předpis polokoule s poloměrem ${\displaystyle R}$. Budeme se snažit odhadnout hodnotu integrálu této funkce na množině ${\displaystyle \{(x,y),|x|<=R\land |y|<=R\}}$. Rozšíříme definiční obor funkce pro ${\displaystyle x,y\in \mathbb {R} }$. Takže

${\displaystyle {\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 ${\displaystyle 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<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;
}

Části vzorkovaného intervalu, kde má ${\displaystyle 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í. 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 ${\displaystyle f(x)}$, tedy bychom měli

${\displaystyle p(x)={\frac {f(x)}{N}},}$

kde ${\displaystyle N}$ je normalizační faktor definovaný jako

${\displaystyle N=\int \limits _{\Omega }f(x)\,\mathrm {d} x,}$

který zajišťuje, že integrál ${\displaystyle p(x)}$ přes celou doménu je ${\displaystyle 1}$ a tedy že ${\displaystyle p(x)}$ je hustota pravděpodobnosti. Rozptyl odhadu by pak byl nulový, jak je vidět z následujících úprav.

${\displaystyle {\begin{aligned}&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{aligned}}}$

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.

Jednou z možností jak snížit rozptyl odhadu je využití takzvané řídící funkce ${\displaystyle g}$. Důležité je aby tato funkce byla analyticky integrovatelná a zároveň aproximovala námi integrovanou funkci ${\displaystyle f}$.

${\displaystyle 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 ${\displaystyle g}$ aproximuje funkci ${\displaystyle 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.

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ší.

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.

Nejdříve mějme funkci

${\displaystyle L(z)={\begin{cases}1,&{\text{pro }}0\leq z_{1},\dots ,0\leq z_{s}\\0,&{\text{jinak,}}\end{cases}}}$

kde ${\displaystyle z\in \mathbb {R} ^{s}}$. Objem kvádru definovaného funkcí ${\displaystyle L}$ v ${\displaystyle A\subseteq \mathbb {R} ^{s}}$ je

${\displaystyle V(A)=\prod _{j=1}^{s}{v_{j}}.}$

Víme, že Monte Carlo odhad objemu tohoto kvádru je roven

${\displaystyle {\frac {1}{N}}\sum _{i=1}^{N}{L(z_{i})}={\frac {m(A)}{N}},}$

kde ${\displaystyle N}$ je počet vzorků, ${\displaystyle z_{i}}$ jsou jednotlivé vzorky a ${\displaystyle m(A)}$ značí počet vzorků, které spadly kvádru ${\displaystyle 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

${\displaystyle D^{*}(z_{1},\dots ,z_{N})=\sup _{A}{\left|{\frac {m(A)}{N}}-V(A)\right|}.}$

Diskrepance slouží jako míra uniformity dané množiny. Konverguje k nule pro ${\displaystyle 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í. Jedna ze sekvencí s nízkou diskrepancí je Van der Corputova a její rozšíření do prostorů s vyšší dimenzí od Haltona nebo Hammersleyho.