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]

Ilustrace Monte Carlo integrace. Použití náhodných vzorků pro určení toho jak velkou část čtverce zabírá kruh

ÚlohaEditovat

Budeme odhadovat hodnotu integrálu   funkce   definovanou jako

 

EstimátoryEditovat

Primární estimátorEditovat

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átoruEditovat

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átoruEditovat

Rozptyl se používá jako měřítko pro určení chyby estimátoru.

 

Pokud   bude podobná   bude výsledný rozptyl malý.

Sekundární estimátorEditovat

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átoruEditovat

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, 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 odhaduEditovat

Vzorkování podle důležitosti (Importance sampling)Editovat

Čá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)Editovat

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)Editovat

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 CarloEditovat

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 diskrepanceEditovat

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í diskrepanceEditovat

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.

ReferenceEditovat

  1. Press et al, 2007, Chap. 4.
  2. Press et al, 2007, Chap. 7.
  3. Newman, 1999, Chap. 2.