STANOVENÍ HODNOTY ČÍSLA $\pi$ METODOU MONTE CARLO

Abstrakt. Náhodné generování bodů v rovině umožňuje simulovat chování mnoha fyzikálních systémů. Článek pojednává o metodě Monte Carlo a její aplikaci při stanovení hodnoty čísla $\pi$.

Představme si následující problém: Máme čtverec o straně délky $a$, do kterého je vepsán kruh. Jaká je pravděpodobnost (jev $\mathsf{A}$), že náhodně zvolený bod leží uvnitř vepsaného kruhu? Předpokládejme, že každý bod může být vybrán se stejnou pravděpodobností. Velikost množiny příznivých výsledků a všech výsledků bude charakterizována svým obsahem. Potom užitím geometrické definice pravděpodobnosti dostáváme $$P(\mathsf{A})=\frac{m(\mathsf{A})}{m(\Omega)}=\frac{\frac{\pi a^2}{4}}{a^2}=\frac{\pi}{4}.$$

Budeme-li náhodný pokus: „vybrat náhodně bod ze čtverce o straně délky $a$“ opakovat $n$ krát, můžeme zaznamenávat, kolikrát se bod nacházel uvnitř kruhu $n(\mathsf{A})$ a vyjadřovat poměr $n(\mathsf{A})/n$, tedy tzv. relativní četnost výskytu jevu $\mathsf{A}$. A co víc, dostatečně velkým opakováním náhodného pokusu můžeme zjišťovat hodnotu matematické konstanty $\pi$. Platí totiž, že $$\pi\sim 4\frac{n(\mathsf{A})}{n},$$

neboť $\displaystyle \lim_{n\to \infty} \frac{n(\mathsf{A})}{n}=P(\mathsf{A})$.

Popsaná metoda se nazývá Monte Carlo právě proto, že se využívá náhodných procesů ke hledání řešení konkrétních problémů (město Monte Carlo je známé svými casiny a hazardem). Touto metodou lze v praxi rychle určovat například obsahy složitých geometrických obrazců (výpočty integrálů, apod.).

Pro simulování náhodného výběru bodu ze čtverce o straně délky $a$ na počítači dobře poslouží generátor náhodných čísel. Postupů, kterými počítače generují náhodná čísla, je mnoho. Pro naše účely dobře poslouží standardní generátory, které nalezneme již v tabulkových procesorech nebo, jako v našem případě, v programovacím prostředí Free Pascalu. Zvolme délku strany $a=1$. Generátor náhodných čísel vybere pro každý pokus dvě čísla z intervalu $\langle 0,1\rangle$, kterým přiřadí význam souřadnic $x$ a $y$ bodu v jednotkovém čtverci.

Monte Carlo simulace použitá ke hledání čísla pí.

Ukázka výsledku náhodného generování $n$ bodů v jednotkovém čtverci, pro $n=100$, $n=1000$, $n=2000$, $n=3000$, $n=4000$ a $n=5000$.

Ve druhé části článku se budeme zabývat důležitou otázkou, kterou je přesnost výpočtu čísla $\pi$ metodou Monte Carlo. Je potřeba dodat, že následující text vyžaduje od čtenáře základní orientaci v teorii pravděpodobnosti, tedy znalost pojmů, jako je například náhodná veličina, střední hodnota, rozptyl, rozdělení pravděpodobnosti, centrální limitní věta a další. Není tedy primárně určen pro středoškolské studenty, ale i pro ně může být rozhodně přínosný.

Abychom mohli provést kvalifikovaný odhad chyby, se kterou se hodnota čísla $\pi$ náhodnými pokusy získává, potřebujeme do celého problému proniknout poněkud hlouběji. Začneme tím, že si zavedeme náhodnou veličinu $M$, která bude vyjadřovat, kolik testovacích bodů z daného počtu pokusů $n$ padne do kruhu o průměru $d=1$. Náhodná veličina $M$ tak může nabývat všech hodnot od $0$ do $n$. Označme jako $p$ pravděpodobnost, při které náhodně zvolený bod padne do daného obrazce (v našem případě kruhu; lze zobecnit i na jiné útvary). Pravděpodobnost, že náhodně vybraný bod bude ležet vně obrazec, je potom $1-p$. Náhodná veličina $M$ má potom binomické rozdělení pravděpodobnosti $M\sim\mathsf{Bi}(n,p)$, neboť pravděpodobnost jevu, že do obrazce padne právě $k$ náhodně vybraných bodů je dána předpisem $$P[M=k]=\binom{n}{k}p^k(1-p)^{n-k}.$$ Je obecně známo, že střední hodnota binomického rozdělení $M$ je $\mathsf{E}M=np$ a jeho rozptyl $\mathsf{D}$ je $\mathsf{D}M=np(1-p)$.

Pro naše účely bude výhodné namísto s náhodnou veličinou $M$ pracovat s náhodnou veličinou $\frac{M}{n}$, která vyjadřuje relativní četnost úspěchu jevu, při kterém se náhodně vybraný bod nachází uvnitř obrazce. Pro veličinu $\frac{M}{n}$ platí $\lim_{n\to\infty}\frac{M}{n}=p$. Využitím obecných vlastností střední hodnoty a rozptylu dostáváme \begin{align*} \mathsf{E}\left(\frac{M}{n}\right)&=\frac{1}{n}\mathsf{E}M=\frac{1}{n}\cdot np=p,\\ \mathsf{D}\left(\frac{M}{n}\right)&=\frac{1}{n^2}\mathsf{D}M=\frac{np(1-p)}{n^2}=\frac{p(1-p)}{n}. \end{align*}

Pro libovolnou náhodnou veličinu $X$ se střední hodnotou $\mathsf{E}X$ a rozptylem $\mathsf{D}X$ je pravděpodobnost, že absolutní hodnota $|X-\mathsf{E}X|$ bude menší než libovolné $\varepsilon>0$ omezena tzv. Čebyševovou nerovností II. typu, kterou nejčastěji píšeme ve tvaru $$P\left[|X-\mathsf{E}X|<\varepsilon\right]\geq 1-\frac{\mathsf{D}X}{\varepsilon^2}.$$ Hodnotu parametru $\varepsilon$ můžeme v našem případě interpretovat jako velikost chyby, se kterou se liší vypočtená (změřená) střední hodnota od skutečné. Aplikací Čebyševovy nerovnosti na náhodnou veličinu $\frac{M}{n}$ dostáváme $$P\left[\left|\frac{M}{n}-p\right|<\varepsilon\right]\geq 1-\frac{p(1-p)}{n\varepsilon^2}.$$ Jelikož na pravé straně nerovnosti se nachází pravděpodobnost $p$, kterou obecně neznáme, využijeme následující odhad. Funkce $f(p)=p(1-p)$ je kvadratická funkce, která nabývá svého maxima v bodě $[\frac{1}{2};\frac{1}{4}]$, a je tudíž omezená shora. Můžeme tedy provést horní odhad ve tvaru $p(1-p)\leq\frac{1}{4}$. Celkem po dosazení platí $$P\left[\left|\frac{M}{n}-p\right|<\varepsilon\right]\geq 1-\frac{p(1-p)}{n\varepsilon^2}\geq 1-\frac{1}{4n\varepsilon^2}.$$

Pro odhad pravděpodobnosti na levé straně zavedeme parametr $\alpha$, který představuje hladinu významnosti. Za hladinu významnosti volíme nejčastěji hodnotu $\alpha=0,05$, tedy $5\;\%$. Požadujeme, aby pravděpodobnost, že se od sebe naměřená a skutečná velikost střední hodnoty liší o méně než $\varepsilon$ byla $$P\left[\left|\frac{M}{n}-p\right|<\varepsilon\right]=1-\alpha.$$

Jednoduchými úpravami postupně dostáváme \begin{align*} 1-\alpha&\geq 1-\frac{1}{4n\varepsilon^2}\\ \alpha &\leq\frac{1}{4n\varepsilon^2}\\ \varepsilon^2 &\leq\frac{1}{4n\alpha},\\ \end{align*}

což po odmocnění dává $$\varepsilon\leq\frac{1}{2\sqrt{n\alpha}}.$$

Jiný odhad můžeme nalézt za použití tzv. Moivre-Laplaceovy centrální limitní věty. Nechť $X\sim\mathsf{Bi}(n,p)$, potom $$\lim_{n\to\infty}P\left[a<\frac{X-np}{\sqrt{np(1-p)}}< b\right]=\Phi(b)-\Phi(a),$$

kde $\Phi(a)$, resp. $\Phi(b)$ je hodnota distribuční funkce v bodě $a$, resp. $b$ normovaného normálního rozdělení $\mathsf{N}(0,1)$. Moivre-Laplaceova centrální limitní věta popisuje konvergenci $\frac{X-np}{\sqrt{np(1-p)}}\to\mathsf{N}(0,1)$.

Hodnoty $a$ a $b$ budeme volit symetricky od střední hodnoty $0$ tak, aby pravděpodobnost $P[a< Y< b]=1-\alpha$, kde $Y\sim\mathsf{N}(0,1)$. To odpovídá volbě $a=z_{\frac{\alpha}{2}}$ a $b=z_{1-\frac{\alpha}{2}}$, kde funkce $z$ je tzv. kvantilová funkce normovaného normálního rozdělení nebo též inverzní funkce k distribuční funkci $\Phi$. Potom totiž platí $$\lim_{n\to\infty}P\left[z_\frac{\alpha}{2}<\frac{X-np}{\sqrt{np(1-p)}}< z_{1-\frac{\alpha}{2}}\right]=\Phi(z_{1-\frac{\alpha}{2}})-\Phi(z_\frac{\alpha}{2})=\left( 1-\alpha+\frac{\alpha}{2}\right)-\frac{\alpha}{2}=1-\alpha.$$ Nyní lze upravit levou stranu, resp. výraz pro pravděpodobnost. Položme $X=M$ a využijeme skutečnost, že $z_\frac{\alpha}{2}=-z_{1-\frac{\alpha}{2}}$. Postupně dostáváme \begin{multline*} P\left[z_\frac{\alpha}{2}<\frac{M-np}{\sqrt{np(1-p)}}< z_{1-\frac{\alpha}{2}}\right]= P\left[-z_{1-\frac{\alpha}{2}}<\frac{M-np}{\sqrt{np(1-p)}}< z_{1-\frac{\alpha}{2}}\right]=\\ =P\left[\left|\frac{M-np}{\sqrt{np(1-p)}}\right|< z_{1-\frac{\alpha}{2}}\right] =P\left[\left|\frac{\frac{M-np}{n}}{\frac{\sqrt{np(1-p)}}{n}}\right|< z_{1-\frac{\alpha}{2}}\right]=\\ =P\left[\left|\frac{M-np}{n}\right|< \sqrt{\frac{np(1-p)}{n^2}}z_{1-\frac{\alpha}{2}}\right]=P\left[\left|\frac{M}{n}-p\right|< \sqrt{\frac{p(1-p)}{n}}z_{1-\frac{\alpha}{2}}\right]. \end{multline*} Označíme-li $\left|\frac{M}{n}-p\right|$ jako $\varepsilon$, jež má význam velikosti odchylky hledané pravděpodobnosti od statisticky zjištěné, lze za použití již uvedeného odhadu $p(1-p)\leq \frac{1}{4}$ psát $$\varepsilon<\frac{z_{1-\frac{\alpha}{2}}}{2\sqrt{n}},$$

čímž nalézáme druhý předpis pro odhad chyby užitím Moivre-Laplaceovy centrální limitní věty.

Shrnutím závěrů z obou postupů můžeme konstatovat, že velikost chyby $\varepsilon$ při metodě Monte Carlo se řádově řídí předpisem $\varepsilon\sim\frac{1}{\sqrt{n}}$. V praxi tento závěr můžeme interpretovat tak, že k dosažení přesnosti $5$ desetinných míst ($\varepsilon=0,00001$), by mělo postačovat $10^{10}$ bodů, což je 10 miliard bodů. To je hodnota, kerá je na běžném počítači s běžným systémovým vybavením nedosažitelná. Z tohoto důvodu nemůžeme metodu Monte Carlo považovat za způsob, jakým dosahovat vysoké přesnosti hodnoty čísla $\pi$. V praxi však její uplatnění nalézáme tam, kde je vyžadována spíše relativně veliká rychlost nalezení řešení, než jeho vysoká přesnost.

Výše uvedený graf znázorňuje závislost přesnosti výpočtu čísla $\pi$ na počtu náhodně vybraných bodů $n$ (osa $x$ je v logaritmické škále). Sytě je vyznačena skutečná hodnota čísla $\pi$. Znázorněny jsou i odhady průběhů chyby, které byly výše odvozeny jak pomocí Čebyševovy nerovnosti, tak Moivre-Laplaceovy centrální limitní věty. Vzhledem k tomu, že chyba se může od čísla $\pi$ odklonit ke kladnému i zápornému směru, jsou v grafu vyznačeny křivky ve tvaru $\pi\pm\varepsilon$. Hladina významnosti je 5 %, kvantil normálního rozdělení $z_{1-\frac{\alpha}{2}}$ je $z_{1-\frac{0,05}{2}}=1,96$. V grafu je zobrazeno rozložení experimentálně získané série dat (posloupnost $a_n$), která byla vygenerována s ekvidistantním krokem $200$ bodů. Průměrování různých sérií nebylo provedeno, tzn. jedná se o ukázku právě jedné série dat. Z grafu je patrné, že přesnost počítačem vygenerovaných hodnot čísla $\pi$ metodou Monte Carlo kopíruje trend teoreticky odvozených předpokladů o velikosti chyby takto získané hodnoty.


Monte Carlo simulace ke stažení
Monte Carlo ikona Samospustitelný program ve formátu exe. Jedná se o konzolovou aplikaci, kterou jsem naprogramoval v jazyce Pascal. Po zadání počtu pokusů se náhodně vybere požadovaný počet bodů a vypočítá se hledaná hodnota čísla $\pi$. S novými operačními systémy lze počítat s menší podporou vašeho počítače, někdy je nutné spustit program v režimu kompatibility Windows XP.
© Michal Řepík