Metoda Monte Carlo

Z Wikipedie, otevřené encyklopedie

Monte Carlo je třída algoritmů pro simulaci systémů. Jde o stochastické metody používající náhodná čísla. Typicky využívány pro výpočet integrálů, zejména vícerozměrných, kde běžné metody nejsou efektivní. Metoda Monte Carlo má široké využití od simulaci experimentů přes počítání určitých integrálů až třeba řešení diferenciálních rovnic. Základní myšlenka této metody je velice jednoduchá, chceme určit střední hodnotu veličiny, která je výsledkem náhodného děje. Vytvoří se počítačový model toho děje a po proběhnutí dostatečného množství simulací se mohou data zpracovat klasickými statistickými metodami, třeba určit průměr a směrodatnou odchylku.

Rozlišují se dvě varianty metody Monte Carlo, analogový a neanalogový model.

Obsah

1 Analogový model
2 Neanalogový model
- 2.1 Integrace metodou Monte Carlo
- 2.2 Příklad
3 Přesnost metody Monte Carlo
4 Další použití metody
5 Související články

[editovat] Analogový model

Musíme umět modelovat celou situaci na počítači, to znamená například znát všechna pravděpodobnostní rozložení zkoumaných jevů a fyzikální zákonitosti, kterými se řídí. Provedením této simulace získáme výsledek,realizaci jakési náhodné veličiny $ξ$ . Tuto simulaci spustíme n-krát a získáme soubor historií x₁ ...x_n . Odhad střední hodnoty $ξ$ se určí :

$\xi \approx \frac{1}{N}\sum_{i=1}^N x_i$

a směrodatná odchylka se určí jednoduše jako

$\sigma = \sqrt{ \frac{\sum_{i=1}^n (\bar{x} - x_i)^2}{n-1}}$

[editovat] Neanalogový model

Tak se nazývá případ, při výpočtu nepoužíváme model reálného děje. Například výpočet určitého integrálu, případně obsahu ohraničeného útvaru.

[editovat] Integrace metodou Monte Carlo

Hlavní výhodou integrace touto metodou je její jednoduchost. Na druhou stranu, přesnost výpočtu je velice nízká obzvláště pokud má integrovaná funkce úzká maxima. Pokud se potřebuje přesný výsledek je výhodnější použít jiné metody.

Pokud určíme hodnotu funkce f v N náhodných bodech ležících v integrační oblasti V pak platí vztah

$\int_V f(\mathbf{x})dV \approx Vf(\mathbf{x})\pm\frac{\sigma V}{ \sqrt{n}}$

$\sigma = \sqrt{\frac{\sum_{i=1}^n (\overline{y} - y_i)^2}{n}}$

kde x označuje náhodný bod v oblasti V, $\bar{y}$ označuje průměr funkčních hodnot a $y i$ jednotlivé výsledky.

Integrování oblasti V v oblasti V'

Pokud budeme tento integrál počítat metodou Monte Carlo musíme oblast V uzavřít do co nejmenší oblasti V'. Případně oblast V rozdělit na více částí a integrovat je samostatně. Nyní zavedeme funkce $\tilde{f}$ takto

$\tilde{f}(\mathbf{x})=\begin{cases} 0 & \mathbf{x} \notin V \\ f(\mathbf{x})& \mathbf{x} \in V \end{cases}$

potom vygenerujeme N náhodných bodů ležících v oblasti V' a přibližnou hodnotu integrálů určíme jako

$I \approx \frac{V'}{N}\sum_{i=1}^N f'(\mathbf{x})$

a směrodatná odchylka se určí analogicky jako v předešlém případě

[editovat] Příklad

Je dána funkce $f(x,y) = \sqrt{R^2-x^2-y^2}$ tedy předpis polokoule. Oblast V bude kružnice s poloměrem R a V' čtverec o straně 2R který ji těsně obepíná. funkce $\tilde{f}$ bude definována takto

$\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}$

$\sigma = \sqrt{\frac{1}{n} \sum_{i=1}^n (x_i - \overline{x})^2} = \sqrt{ \frac{1}{n} \left(\sum_{i=1}^n x_i^2 - \bar{x}\sum_{i=1}^n x_i \right)}$

Během simulace si stačí zaznamenat pouze hodnoty $r = \sum_{i=1}^N x_i$ a $t = \sum_{i=1}^N x_i^2$

Zde je jednoduchý příklad jak by to mohlo vypadat v praxi

#include<math.h>
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 &u)
{
        double r,t,x;
        r = t = 0;
        for(int i= 0;i < n;i++){
                x = f(2*R*rand()/RAND_MAX-R,2*R*rand()/RAND_MAX-R);
                r+= x;
                t+= x*x;
        }
        I = 4*R*R*r/n;
        u = I*sqrt(t/(r*r)-1.0/n);
}

[editovat] Přesnost metody Monte Carlo

K odhadu chyby výsledku získaného metodou Monte Carlo se většinou používá střední kvadratická chyba aritmetického průměru. Chyba výsledku získaného pomocí n historií je úměrná $1/\sqrt{n}$ . Takže aby se zlepšil výsledek o jeden řád musí se počet historií zvýšit alespoň o dva řády. Abychom získali výsledek s přesností na 6 desetinných míst, což odpovídá přesnosti jiných metod musíme získat 10¹² historií.