Metody Monte Carlo: wykorzystanie losowości do rozwiązywania problemów deterministycznych

Szacowanie π za pomocą losowych punktów

Klasyczna demonstracja Monte Carlo zaczyna się od geometrii. Rozważmy okrąg jednostkowy — o promieniu 1, wyśrodkowany w początku układu — wpisany w kwadrat 2×2. Pole okręgu wynosi π·r² = π. Pole kwadratu wynosi 4. Ich stosunek to dokładnie π/4.

Teraz wygeneruj N losowych punktów równomiernie rozłożonych w kwadracie, o współrzędnych (x, y), gdzie zarówno x, jak i y mieszczą się w zakresie od −1 do 1. Dla każdego punktu sprawdź, czy x² + y² ≤ 1. Jeśli tak, punkt leży wewnątrz okręgu; nazwijmy liczbę takich punktów M. Wówczas:

π ≈ 4 × M / N

Przy 1000 punktach możesz otrzymać π ≈ 3,1 — z grubsza, ale rozpoznawalnie. Przy 1 milionie punktów zazwyczaj osiągasz około 3,141 — dokładność do trzech miejsc po przecinku. Błąd oszacowania Monte Carlo skaluje się jako 1/√N: aby zyskać jedno dodatkowe miejsce dziesiętne dokładności, potrzeba 100 razy więcej próbek. Ta powolna zbieżność to cena prostoty.

Co czyni to godnym uwagi, to nie efektywność — istnieją znacznie szybsze sposoby obliczania π — ale zasada, którą to demonstruje: losowe próbkowanie może wydobyć dokładną informację matematyczną zarówno z geometrii, prawdopodobieństwa, jak i fizyki.

Całkowanie w wysokich wymiarach

Tradycyjne metody całkowania numerycznego, takie jak reguła trapezów czy reguła Simpsona, działają poprzez umieszczenie siatki punktów obliczeniowych nad domeną całkowania. W jednym wymiarze N punktów daje błąd, który zazwyczaj skaluje się jako 1/N² lub lepiej. Ale istnieje śmiertelny problem: w d wymiarach siatka z N punktami na oś wymaga N^d obliczeń łącznie.

Dla d = 10 wymiarów przy 100 punktach na oś to 10²⁰ obliczeń — daleko poza możliwościami jakiegokolwiek komputera. Dla d = 100 liczba wymaganych punktów siatki jest astronomicznie większa niż liczba atomów w obserwowalnym wszechświecie. To klątwa wymiarowości.

Całkowanie Monte Carlo jest na nią odporne. Wystarczy próbkować N losowych punktów w d-wymiarowej domenie i uśrednić wartości funkcji podcałkowej. Błąd zawsze wynosi 1/√N, niezależnie od liczby wymiarów. Monte Carlo staje się relatywnie bardziej efektywne niż jakakolwiek metoda siatkowa wraz ze wzrostem wymiarowości, dlatego dominuje w dziedzinach takich jak:

• Chemia kwantowa (obliczanie energii elektronowych cząsteczek w setkach wymiarów)
• Fizyka statystyczna (funkcje podziału nad ogromnymi przestrzeniami konfiguracyjnymi)
• Instrumenty pochodne finansowe (wycena opcji przy modelach skorelowanych wielu aktywów)
• Transport promieniowania (śledzenie milionów ścieżek fotonów lub neutronów)

Początki w Projekcie Manhattan

Nazwa „Monte Carlo” została ukuta przez Stanisława Ulama i Johna von Neumanna w latach 40. XX wieku, na cześć słynnej dzielnicy kasyn w Monako — nawiązanie do centralnej roli przypadku. Okazją był Projekt Manhattan, a konkretnym problemem była dyfuzja neutronów przez materiał rozszczepialny.

Gdy neutron przemieszcza się przez uran lub pluton, przechodzi przez złożony łańcuch zdarzeń rozpraszania, absorpcji i rozszczepienia, których prawdopodobieństwa zależą od właściwości materiału i energii neutronu. Napisanie równań analitycznych śledzących wszystkie możliwe ścieżki jest niewykonalne. Ale symulowanie losowej ścieżki pojedynczego neutronu jest proste — a symulowanie milionów z nich statystycznie odtwarza zagregowane zachowanie reakcji łańcuchowej.

Enrico Fermi wykonywał w rzeczywistości mentalne obliczenia Monte Carlo jeszcze wcześniej, wykorzystując swoją słynną zdolność do szacowania przekrojów czynnych neutronów przez analogię i mentalną randomizację. Ale to Ulam sformalizował metodę po wyzdrowieniu z choroby i spędzeniu czasu na grze w pasjansa — zastanawiając się, jaki odsetek gier można wygrać, nie rozgrywając ich wszystkich do końca.

Von Neumann natychmiast dostrzegł zastosowanie do problemów neutronowych i, co kluczowe, rozpoznał, że nowo zbudowany komputer ENIAC mógłby uruchamiać takie symulacje mechanicznie. Pierwsze poważne obliczenia Monte Carlo działały na ENIAC w 1947 roku, a metoda ta jest odtąd niezastąpiona w fizyce jądrowej, projektowaniu broni i inżynierii reaktorów.

Markov Chain Monte Carlo

Głębszy problem pojawia się, gdy chcesz próbkować z rozkładu prawdopodobieństwa, którego nie możesz bezpośrednio odwrócić ani znormalizować. Jest to powszechne w statystyce bayesowskiej, gdzie rozkład a posteriori parametrów modelu może być wysokowymiarową, nieznormalizowaną gęstością proporcjonalną do iloczynu wiarygodności i rozkładu a priori.

Markov Chain Monte Carlo (MCMC) rozwiązuje ten problem, konstruując łańcuch Markowa — sekwencję losowych stanów, w której każdy stan zależy tylko od poprzedniego — którego rozkład stacjonarny jest równy rozkładowi docelowemu. Po okresie „rozgrzewki”, podczas którego łańcuch osiąga równowagę, kolejne stany łańcucha są (skorelowanymi) próbkami z rozkładu docelowego.

Fundamentalnym algorytmem jest Metropolis-Hastings (1953, rozszerzony w 1970):

1. Zacznij w stanie x.
2. Zaproponuj nowy stan x' z rozkładu propozycji q(x'|x).
3. Oblicz współczynnik akceptacji α = [p(x') · q(x|x')] / [p(x) · q(x'|x)], gdzie p to gęstość docelowa.
4. Zaakceptuj x' z prawdopodobieństwem min(1, α); w przeciwnym razie pozostań w x.
5. Powtórz.

Geniusz tego algorytmu polega na tym, że stała normalizująca p redukuje się w stosunku α — trzeba oceniać jedynie nieznormalizowaną gęstość. MCMC jest dziś centralne dla wnioskowania bayesowskiego, przewidywania struktury białek oraz symulowania modelu Isinga w fizyce statystycznej, gdzie bezpośrednie próbkowanie z rozkładu Boltzmanna jest niemożliwe.

Monte Carlo w finansach

Model Blacka-Scholesa zakłada, że ceny akcji podążają za geometrycznym ruchem Browna — czystym modelem analitycznym z zamkniętą formułą wyceny opcji. Rzeczywiste portfele finansowe są znacznie bardziej złożone: dziesiątki skorelowanych aktywów, zmienność, która sama fluktuuje (zmienność stochastyczna), procesy skokowe, wypłaty zależne od ścieżki i ograniczenia regulacyjne.

Symulacja Monte Carlo obsługuje to wszystko naturalnie. Symuluj tysiące lub miliony możliwych trajektorii rynkowych na przestrzeni czasu życia opcji, oblicz wypłatę przy wygaśnięciu w każdym scenariuszu, zdyskontuj do dzisiaj i uśrednij. Wynikiem jest godziwa wartość opcji przy założeniach twojego modelu.

To podejście napędza również obliczenia Value at Risk (VaR): symuluj wartość portfela w tysiącach możliwych scenariuszy rynkowych na horyzoncie jedno- lub dziesięciodniowym i zgłoś poziom straty przekraczany zaledwie w 1% (lub 5%) scenariuszy. Po kryzysie finansowym z 2008 roku bardziej wyrafinowane ramy testów warunków skrajnych rozszerzyły to podejście na setki scenariuszy makroekonomicznych — zasadniczo wielkoskalowe Monte Carlo nad przestrzeniami stanów gospodarczych.

Prawo wielkich liczb i granice błędu

Zbieżność Monte Carlo jest gwarantowana przez prawo wielkich liczb: gdy N rośnie, średnia próbki zbiega niemal na pewno do prawdziwej średniej. Centralne twierdzenie graniczne mówi więcej: rozkład średniej próbki jest w przybliżeniu normalny, z odchyleniem standardowym σ/√N, gdzie σ to odchylenie standardowe próbkowanej wielkości.

Daje to dokładną, obliczalną granicę błędu — w przeciwieństwie do wielu deterministycznych metod numerycznych, gdzie analiza błędu wymaga założeń o gładkości. Zawsze możesz oszacować swój błąd Monte Carlo, uruchamiając symulację wielokrotnie i mierząc zmienność.

Aby zmniejszyć wariancję i przyspieszyć zbieżność, praktycy stosują techniki takie jak:

• Próbkowanie ważne (importance sampling): próbkuj intensywniej z regionów, które wnoszą największy wkład do całki.
• Zmienne antytetyczne: dla każdej losowej próbki użyj też jej „odbicia lustrzanego”, by zniwelować skorelowane błędy.
• Zmienne kontrolne: odejmij skorelowaną wielkość, której prawdziwa średnia jest znana.
• Quasi-Monte Carlo: zastąp sekwencje pseudolosowe sekwencjami o niskiej rozbieżności (Sobol, Halton), które wypełniają przestrzeń bardziej równomiernie. W praktyce zbieżność QMC może zbliżyć się do O(1/N) — dramatycznie szybciej niż standardowe 1/√N — szczególnie w umiarkowanych wymiarach.

Trwała siła Monte Carlo nie leży pomimo jej losowości, ale dzięki niej: losowość zapewnia nieobciążoną, niezależną od wymiaru ścieżkę do odpowiedzi, których metody deterministyczne nie są w stanie osiągnąć.