Chcemy obliczyć całkę \(\di J = \int _a^b f(x)\,dx\).
Metoda I (uogólnienie metody obliczania \(\pi \)), przykład 12.1.
Załóżmy, że \(f\) jest ograniczona, dla uproszczenia, że \(0 \le f(x) \le c\), gdzie \(c \in \r \).
Niech \(U_1 = (U_{11},U_{12}), \dots , U_n = (U_{n1},U_{n2}) \) będą niezależnymi wektorami losowymi o rozkładzie jednostajnym na prostokącie \([a,b]\times [0,c]\). Określamy zmienne losowe:
\[ X_i = \left \{\begin {array}{ll} 1, \mbox { gdy } & f(U_{i1}) \le U_{i2} \\ 0, \mbox { gdy } & f(U_{i1}) > U_{i2} \end {array} \right . \]
Oczywiście zmienne losowe \(X_i\) są niezależne i mają rozkład \((0,1,p)\), gdzie \(\di p = \frac {J}{(b-a)c}\). Mamy więc \(m = E(X_i) = p\), \(\sigma ^2 = D^2(X_i) = p(1-p)\). Ponieważ \(\bar {X}_n\) jest nieobciążonym i silnie zgodnym estymatorem \(m\), to \(\bar {J}_n =(b-a)c\bar {X}_n\) jest nieobciążonym i silnie zgodnym estymatorem całki \(J\).
Znajdziemy dwustronny przedział ufności dla \(J\) na poziomie \(1 - \alpha \). Ponieważ dla \(m\) takim przedziałem jest \((\bar {X}_n - \ve , \bar {X}_n + \ve )\), gdzie \(\ve = \frac {\sigma _X}{\sqrt {n}} \Phi ^{-1}\left (1 - \frac {\alpha }{2}\right ). \), to dla \(J = (b-a)c m\) jest przedział: \((\bar {J}_n - \delta _1,\bar {J}_n + \delta _1)\), gdzie:
\[\di \delta _1 = (b-a)c \ve = (b-a)c \frac {\sqrt {\frac {J}{(b-a)c}(1 - \frac {J}{(b-a)c})}}{\sqrt {n}} \Phi ^{-1}\left (1 - \frac {\alpha }{2}\right ) = \]
\[\sqrt {J(b-a)c - J^2} \frac {1}{\sqrt {n}} \Phi ^{-1}\left (1 - \frac {\alpha }{2}\right ).\]
Metoda II. Zauważmy, że
\[\di J = \int _{a}^b f(x)\,dx = (b-a)\int _{\r }\frac {1}{b-a} I_{[a,b]}(x) f(x)\,dx = (b-a)E(f(X)),\]
gdzie \(X\) jest zmienną losową o rozkładzie \(U(a,b)\).
Stosując Przykład 12.10 do zmiennej \(Y = f(X)\) mamy wzór na estymator nadziei \(m_Y = E(Y)\) oraz umiemy znaleźć przedziały ufności dla \(m_Y\).
Estymatorem \(\hat {J}\) dla całki \(J\) jest więc \(\hat {J} = (b-a) \bar {Y}_n= \frac {b-a}{n} \sum _{i=1}^n Y_i\), gdzie \(Y_i = f(X_i)\), a \(X_1,X_2, \dots , X_n\) jest próbką prostą z rozkładu \(P_X = U(a,b)\).
Przedziałem ufności dla \(J\) jest więc: \(\di ((b-a)\bar {Y}_n - \delta _2, (b-a)\bar {Y}_n + \delta _2)\), gdzie
\[ \delta _2= (b-a)\frac {\sigma _Y}{\sqrt {n}} \Phi ^{-1}\left (1 - \frac {\alpha }{2}\right ), \]
gdzie \(\di \sigma _Y^2 = D^2(Y) = E(f(X)^2) - E(f(X))^2 = \frac {1}{b-a}\int _a^bf(x)^2\,dx - \frac {J^2}{(b-a)^2}\).
Tak więc:
\[ \delta _2 = \sqrt {(b-a) \int _a^b f(x)^2 \,dx - J^2} \frac {1}{\sqrt {n}} \Phi ^{-1}\left (1 - \frac {\alpha }{2}\right ). \]
Wzory na przedział ufności otrzymane w powyższych metodach zawierają wielkość \(J\), której nie znamy, więc nie mogą być bezpośrednio stosowane.
Zauważmy jednak, że Metoda II dostarcza potencjalnie lepszy (na ogół mniejszy) przedział ufności niż Metoda I. Mianowicie: \(\di (b-a) \int _a^b f(x)^2 \,dx \le J(b-a)c\). Czyli \(\delta _2 \le \delta _1\). Przy nieuważnym doborze \(c\) ta nierówność może być bardzo istotna.
Faktyczny przedział ufności można wyznaczyć dopiero po przeprowadzeniu eksperymentu.
W Metodzie I jest to przedział: \(\di ((b-a)c\bar {X}_n - \hat {\delta }_1,b-a)c\bar {X}_n + \hat {\delta }_1)\), \(\di \hat {\delta }_1 = (b-a)c \frac {\hat {\sigma }_X}{\sqrt {n}} \Phi ^{-1}\left (1 - \frac {\alpha }{2}\right )\), gdzie \(\di \hat {\sigma }_X^2 = \frac {1}{n}\sum _{i=1}^n(X_i - \bar {X}_n)^2\) = \(\di \frac {1}{n}\sum _{i=1}^n X_i^2 - (\bar {X}_n)^2 \).
W Metodzie 2 jest to przedział: \(\di ((b-a)\bar {Y}_n - \hat {\delta }_2, (b-a)\bar {Y}_n + \hat {\delta }_2)\), \(\di \hat {\delta }_2 = (b-a)\frac {\hat {\sigma }_Y}{\sqrt {n}} \Phi ^{-1}\left (1 - \frac {\alpha }{2}\right ), \) gdzie \(\di \hat {\sigma }_Y^2 = \frac {1}{n}\sum _{i=1}^n(Y_i - \bar {Y}_n)^2\) = \(\frac {1}{n}\sum _{i=1}^n Y_i^2 - (\bar {Y}_n)^2 \).
Przykład – 12.16
Policzmy całkę \(J = \int _0^3 x^2\,dx\) stosując obie metody.
Losujemy (Maple) \(n = 10 000\) punktów (z prostokąta \([0,3]\times [0,9]\) w Metodzie I, z odcinka \([0,3]\) w Metodzie II), ustalamy \(\alpha := 0.05\). W tym przypadku umiemy oczywiście wyliczyć bezpośrednio: \(J = 9\), a także \(\delta _1 = 0.249462688356251\) oraz \(\delta _2 = 0.157774057303750\).
W Metodzie I otrzymujemy \(\bar {J}_n = 8.934300000\) oraz
\(\hat {\delta }_1 = 0.249003673187557\).
W Metodzie II otrzymujemy \(\bar {J}_n = 8.91313296682201 \) oraz
\(\hat {\delta }_2 =0.157905405151449\).
Przykład – 12.17
Porównajmy obydwie metody dla bardziej skomplikowanej funkcji.
Niech \(f(x) = x+ \sin x+ \cos 2x\).
Interesuje na całka
\(\di J = \int _0^{4\pi } f(x)\,dx\).
Łatwo bezpośrednio obliczyć, że \(J = 8\pi ^2\).
Uruchamiamy 10 razy metodą 1 i 10 razy metodą 2. Najpierw porównujemy same przybliżenia \(J\), a później odpowiednie przedziały ufności (\(\alpha = 0.05\)). W każdej z metod przeprowadzamy 2000 losowań z rozkładu jednostajnego. Metoda 1 – kolor czerwony, Metoda 2 – kolor niebieski.
1. Opisanymi dwiema metodami oraz innymi metodami Monte Carlo można obliczać całki wielowymiarowe:
\[ \int _D f\, d\mu , \mbox { gdzie } D \subset \rn ,\ \ \mu \mbox { -- dana miara}. \]
2. W wielu przypadkach metody Monte Carlo są jedynymi metodami obliczania takich całek!
3. Oprócz dwóch omówionych powyżej istnieje więcej metod Monte Carlo obliczania całek. Istotnym problemem jest to, aby do danego problemu dobrać metodę, które dają szansę na mały przedział ufności.
4. Wiele problemów obliczeniowych można sprowadzić do odliczania całek.
5. Metody Monte Carlo stosowane są też w innych zagadnieniach. W szczególności metami Monte Carlo owiązywane są złożone zadania optymalizacyjne.