Metoda Monte Carlo obliczania całek + implementacja

Zbiór wzorów, definicji i najczęściej poruszanych problemów z probabilistyki oraz statystyki matematycznej.
Awatar użytkownika
Emiel Regis
Użytkownik
Użytkownik
Posty: 1495
Rejestracja: 26 wrz 2005, o 17:01
Płeć: Mężczyzna
Lokalizacja: Kraków
Podziękował: 71 razy
Pomógł: 225 razy

Metoda Monte Carlo obliczania całek + implementacja

Post autor: Emiel Regis »

Metoda Monte Carlo obliczania całek


Często jest tak, iż wiemy, że istnieje całka oznaczona z funkcji \(\displaystyle{ f}\) jednak nie potrafimy jej analitycznie policzyć. W większości przypadków stosuje się wtedy metody numeryczne, jednak istnieje wobec nich stosunkowo prosta alternatywa. Otóż dla tak deterministycznego problemu jakim jest liczenie całek przychodzi z pomocą probabilistyka. Metoda, którą poniżej omówię jest szczególnie istotna, gdy funkcja, którą całkujemy jest bardzo nieregularna bądź też w przypadku całek wielokrotnych.

Jako ciekawostkę dodam, że pierwszy raz zastosowano na dużą skalę metody oparte o użycie liczb losowych w latach czterdziestych XX wieku. W Los Alamos w ramach "Projektu Manhattan" (nad budową bomby jądrowej pracowali wtedy m.in. John von Neumann, Stanisław Ulam, Richard Feynman) dotyczyły one rozpraszania i absorpcji neutronów.
Tam także powstała nazwa (Monte Carlo - kasyno w Monaco, nawiązanie do losowości gier hazardowych) jako określenie tego typu metod matematycznych.

Przybliżone obliczanie całki
\(\displaystyle{ \boxed{ \int_a^b f(x) dx }}\)

Na początku przedstawię najogólniejszą postać a później jej szczególne przypadki, czasami łatwiej może być rozpocząć czytanie właśnie od szczególnych przypadków, dlatego jeśli ktoś nie widzi poniżej zbieżności wynikającej z MPWL to niech najpierw poczyta dalszy fragment a następnie wróci zobaczyć uogólnienie.

Kluczowa rzecz w metodzie:

Zbieżność ciągów \(\displaystyle{ S_n}\), \(\displaystyle{ S _n'}\), \(\displaystyle{ S_n''}\) gwarantuje nam MPWL.


\(\displaystyle{ \hline}\)

\(\displaystyle{ (X_n)_{n \in \mathbb{N}} \hbox{ - ciąg niezależnych zmiennych losowych o nośniku }(a,b) \hbox{ oraz o gęstości }g}\)

\(\displaystyle{ S_n = \frac{1}{n} \sum_{i=1}^n \frac{f(X_i)}{g(X_i)} \stackrel{n \to \infty}{\longrightarrow} E \left(\frac{f(X_1)}{g(X_1)}\right) =
\int_a^b \frac{f(x)}{g(x)}g(x)dx = \int_a^b f(x) dx}\)



Widzimy więc, że dzięki MPWL powyższa suma zbiega do wartości szukanej całki.


Mogą się nasunąć przynajmniej dwa pytania:

1. Jak duże musi być n aby uzyskane przybliżenie nas satysfakcjonowało?
2. Skąd wziąć wartości zmiennych losowych?

Odpowiedź na pierwsze pytanie oczywiście zależy od tego jak dokładny ktoś chce mieć wynik. Jednak bez wnikania w niuanse można powiedzieć, że metoda ta jest bardzo dokładna już dla stosunkowo niedużych \(\displaystyle{ n}\) (niedużych z punktu widzenia obecnych komputerów).
Metody numeryczne obliczania całek także potrzebują komputerów wiec ich użycie w powyższej metodzie oczywiście nie jest żadnym nadzwyczajnym wymogiem. Warto także zaznaczyć, że metoda Monte Carlo jest szybsza od większości standardowych metod numerycznych.

Odnośnie drugiego pytania:
Nie jest to tematem tego artykułu dlatego tylko wspomnę, że istnieją generatory liczb (pseudo) losowych z przedziału \(\displaystyle{ (0,1)}\).
Mając te liczby potrafimy generować liczby także z innych rozkładów.

Gdy \(\displaystyle{ X_1, \ldots, X_n}\) są niezależnymi zmiennymi losowymi o rozkładzie jednostajnym na przedziale \(\displaystyle{ (0,1)}\) (czyli ich wartości stanowią nasze wygenerowane liczby), a dystrybuanta \(\displaystyle{ F}\) rozkładu, z którego chcemy generować liczby jest odwracalna to \(\displaystyle{ F^{-1}(X_1), \ldots,F^{-1}(X_n)}\) są niezależne i mają rozkład o dystrybuancie \(\displaystyle{ F}\).
Często przy generowaniu liczb z różnych rozkładów korzysta się także z innych metod niż odwracanie dystrybuanty. Jednak dla tematu, który omawiam najważniejsze jest, żeby mieć świadomość, iż da się generować liczby z różnych rozkładów.


\(\displaystyle{ \hline}\)

Przypadek I


\(\displaystyle{ (a,b) = (0,1)}\)

\(\displaystyle{ X_i \sim U(0,1) \Longrightarrow g(x) = 1_{(0,1)}(x)}\)

\(\displaystyle{ S_n' = \frac{1}{n} \sum_{i=1}^n f(X_i) \stackrel{n \to \infty}{\longrightarrow} Ef(X_1) =
\int_{\mathbb{R}} f(x)g(x)dx = \int_0^1 f(x) dx}\)



\(\displaystyle{ \hline}\)

Przypadek II


\(\displaystyle{ X_i \sim U(a,b) \Longrightarrow g(x) = \frac{1}{b-a} \cdot 1_{(a,b)}(x)}\)

\(\displaystyle{ S_n'' = \frac{1}{n} \sum_{i=1}^n f(X_i) \stackrel{n \to \infty}{\longrightarrow} Ef(X_1) =
\int_{\mathbb{R}} f(x)g(x)dx = \frac{1}{b-a} \int_a^b f(x) dx}\)


Przekształcając otrzymamy:

\(\displaystyle{ \frac{b-a}{n} \sum_{i=1}^n f(X_i) \stackrel{n \to \infty}{\longrightarrow} \int_a^b f(x) dx}\)

Można zauważyć, że metoda w punkcie drugim podobnie jak ta z samego początku służy do liczenia dowolnej całki. Nasuwa się wiec pytanie czemu w zależności od przedziału całkowania \(\displaystyle{ [a,b]}\) zawsze nie brać rozkładu jednostajnego na nim. Otóż czasami warto wziąść inne bardziej wymyślne rozkłady z tego względu że wariancja zmiennej \(\displaystyle{ S_n}\) może być mniejsza niż \(\displaystyle{ S_n''}\).


\(\displaystyle{ \hline}\)

Na koniec, aby empirycznie potwierdzić wiarygodność tego co piszę zaprezentuję konkretne wykorzystanie metody. Otóż przedstawię poniżej własną implementację metody Monte Carlo w programie do obliczeń symbolicznych Maple.

Definicja funkcji, którą będziemy całkować:
(wybrałem funkcję elementarną, żeby móc porównać wynik metody Monte Carlo z dokładnym)

Kod: Zaznacz cały

f:=x->ln(x);
Implementacja algorytmu:

Kod: Zaznacz cały

a := 5:
b := 27:
n := 100:
suma := 0:
for i to n do suma := suma+f(stats[random, uniform[a, b]](1)) end do:
wynik := (b-a)*suma/n;
Wynik:

Kod: Zaznacz cały

evalf(int(f(x), x = a .. b))
Funkcja evalf jest po to, żeby niewymierny wynik przybliżyć wymiernym.



Rezultaty wykonania programu.

Wartość całki:
\(\displaystyle{ \boxed{\int_a^b f(x) dx = 58.94040585}}\)

Wartość całki liczona metodą Monte Carlo:
\(\displaystyle{ \hbox{dla } n = 10}\):

\(\displaystyle{ 60.34541047, \ \ \ 61.28887530, \ \ \ 55.27822168, \ \ \ 61.75680023, \ \ \ 56.45305178}\)


\(\displaystyle{ \hbox{dla } n = 100}\):

\(\displaystyle{ 60.30928158, \ \ \ 58.54156748, \ \ \ 59.98247772, \ \ \ 58.65675627, \ \ \ 60.82712744}\)


\(\displaystyle{ \hbox{dla } n = 1000}\):

\(\displaystyle{ 59.01833876, \ \ \ 59.52496339, \ \ \ 58.94087063, \ \ \ 59.40483336, \ \ \ 59.08315534}\)


\(\displaystyle{ \hbox{dla } n = 10000}\);

\(\displaystyle{ 58.82425644, \ \ \ 58.94127837, \ \ \ 59.01296060, \ \ \ 58.96964135, \ \ \ 59.05471524}\)


\(\displaystyle{ \hbox{dla } n = 100000}\);

\(\displaystyle{ 58.96083314, \ \ \ 58.93575466, \ \ \ 58.92137904, \ \ \ 58.92459274, \ \ \ 58.98589176}\)


\(\displaystyle{ \hbox{dla } n = 1000000}\);

\(\displaystyle{ 58.94260775, \ \ \ 58.93068625, \ \ \ 58.93767684, \ \ \ 58.93710768, \ \ \ 58.95010141}\)


\(\displaystyle{ \hline}\)

Powyżej przedstawiłem tylko przykład, zainteresowanym polecam samodzielne testy metody dla innych funkcji.

Nie powinno nikogo dziwić, że czasami dla mniejszych \(\displaystyle{ n}\) otrzymamy dokładniejszy wynik niż dla większych - pamiętajmy, że tutaj mimo wszystko mamy do czynienia z probabilistyką. Natomiast możemy się spodziewać (i tak jest) tendencji wzrostu dokładności wyniku wraz ze wzrostem \(\displaystyle{ n}\).

Polecam także, aby losować punkty nie z rozkładu jednostajnego. Tutaj nie miałem wpływu na wariancję wyniku a można ją zmniejszyć. Mimo wszystko biorąc pod uwagę, że dla \(\displaystyle{ n}\) rzędu \(\displaystyle{ 10}\) czy \(\displaystyle{ 100}\) obliczenia można wykonać nawet na kalkulatorze to wyniki są przyzwoite. Oczywiście jeśli ktoś potrzebuje bardzo dokładnego wyniku to on wtedy prowadzi obliczenia na dużo mocniejszym komputerze, przez co może wziąć \(\displaystyle{ n}\) większe od miliona. Na moim prywatnym komputerze wykonanie programu dla największego \(\displaystyle{ n}\) trwało około \(\displaystyle{ 570s}\).


\(\displaystyle{ \hline}\)
ODPOWIEDZ