Duża część zastosowań metod MC wiąże się z obliczaniem całek lub sum. Typowe zadanie polega na obliczeniu wartości oczekiwanej
![]() |
gdzie jest zmienną losową o rozkładzie prawdopodobieństwa
na przestrzeni
, zaś
. Zazwyczaj
jest podzbiorem wielowymiarowej przestrzeni euklidesowej
lub jest zbiorem skończonym ale bardzo licznym (na przykład zbiorem pewnych obiektów kombinatorycznych). Jeśli
i rozkład
jest opisany przez
gęstość
to całka określająca wartość oczekiwaną jest zwykłą całką Lebesgue'a.
W tym przypadku zatem możemy napisać
![]() |
Równie ważny w zastosowaniach jest przypadek dyskretnej przestrzeni . Jeśli
, to
![]() |
W dalszym ciągu tego rozdziału utożsamiamy rozkład z funkcją
.
Zwróćmy uwagę, że przedstawiamy tu
jako wartość oczekiwaną tylko dla wygody oznaczeń;
w istocie każda całka, suma, (a także prawdopodobieństwo zdarzenia losowego) jest
wartością oczekiwaną.
Na pierwszy rzut oka nie widać czym polega problem! Sumowanie wykonuje każdy kalkulator, całki się sprawnie oblicza numerycznie. Ale nie zawsze. Metody MC przydają się w sytuacjach gdy spotyka się z następującymi trudnościami (lub przynajmniej którąś z nich).
Przestrzeń jest ogromna. To znaczy wymiar
jest bardzo duży lub skończona przestrzeń zawiera astronomicznie dużą liczbę punktów.
Rozkład prawdopodobieństwa jest ,,skupiony” w małej części ogromnej przestrzeni
.
Nie ma podstaw do zakładania, że funkcja jest w jakimkolwiek sensie ,,gładka”
(co jest warunkiem stosowania standardowych numerycznych metod całkowania).
Gęstość rozkładu
jest znana tylko z dokładnością do stałej normującej. Innymi słowy, umiemy obliczać
ale nie znamy stałej
. Czsem zadanie polega właśnie na obliczeni tej stałej (
jest nazwane ,,funkcją podziału” lub sumą statystyczną).
Prosta metoda MC (po angielsku nazywana bardziej brutalnie: Crude Monte Carlo,
czyli CMC) nasuwa się sama. Należy wygenerować niezależnych zmiennych losowych
o jednakowym rozkładzie
i za estymator wartości oczekiwanej
wziąć średnią z próbki,
![]() |
Mocne Prawo Wielkich Liczb gwarantuje, że prawie na pewno, gdy
.
W terminologii statystycznej,
jest mocno zgodnym estymatorem obliczanej wielkości. Zadanie wydaje się rozwiązane. Są jednak dwa zasadnicze kłopoty.
Estymator skonstruowany metodą CMC może się zbliżać do
przerażająco wolno.
Co robić, jeśli nie umiemy losować z rozkładu ?
Zadziwiająco skutecznym sposobem na oba przedstawione wyżej kłopoty
jest losowanie istotne (Importance Sampling, w skrócie IS).
Przypuśćmy, że umiemy losować z rozkładu o gęstości . Zauważmy, że
![]() |
gdzie . Piszemy tu wzory dla całek ale oczywiście dla sum jest tak samo.
Niech
będą niezależnymi zmiennymi losowymi o jednakowym rozkładzie
,
![]() |
(8.1) |
gdzie
![]() |
traktujemy jako wagi wylosowanych punktów .
Z tego co wyżej powiedzieliśmy wynika, że
oraz
prawie na pewno, gdy
. Mamy więc estymator nieobciążony
i zgodny. Milcząco założyliśmy, że
w każdym wylosowanym punkcie, czyli, że
. Z tym zwykle nie ma wielkiego kłopotu.
Ponadto musimy założyć, że umiemy obliczać funkcję
w każdym wylosowanym punkcie.
Jeśli znamy tylko
a nie znamy stałej
to jest kłopot. Możemy tylko obliczyć
. Używamy zatem nieco innej postaci estymatora
IS, mianowicie
![]() |
(8.2) |
Ten estymator wygląda dziwnie, bo w mianowniku występuje estymator jedynki
ale chodzi o to, że we wzorze (8.2) w przeciwieństwie do (8.1)
możemy zastąpić przez
, bo nieznany czynnik
skraca się.
Możemy jeszcze zapisać nasz estymator w zgrabnej formie
![]() |
gdzie są unormowanymi wagami. Trzeba jednak
pamiętać, że to ,,unormowanie” polega na podzieleniu wag przez zmienną losową. W rezultacie
estymator
jest obciążony. Jest jednak mocno zgodny, bo licznik
we wzorze (8.2) po podzieleniu przez
zmierza do
a mianownik do 1,
prawie na pewno.
Estymator
jest znacznie częściej używany niż
. Oprócz uniezależnienia się od stałej normującej, jest jeszcze inny powód. Okazuje się, że (pomimo obciążenia)
estymator
może być bardziej efektywny. Sens tego pojęcia wyjaśnimy w
następnym podrozdziale. Zarówno dobór przykładów jak i ogólny tok rozważań jest zapożyczony z monografii Liu [15].
Naturalne jest żądanie, aby estymator był mocno zgodny, prawie na pewno przy
. Znaczy to, że zbliżymy się dowolnie blisko aproksymowanej wielkości, gdy tylko dostatecznie długo przedłużamy symulacje. Chcielibyśmy jednak wiedzieć coś konkretniejszego,
oszacować jak szybko błąd aproksymacji
maleje do zera i jakie
wystarczy
do osiągnięcia wystarczającej dokładności. Przedstawimy teraz
najprostsze i najczęściej stosowane w praktyce podejście, oparte na tzw. asymptotycznej
wariancji i konstrukcji asymptotycznych przedziałów ufności.
W typowej sytuacji estymator
ma następującą własność, nazywaną
asymptotyczną normalnością:
![]() |
w sensie zbieżności według rozkładu. W konsekwencji, dla ustalonej liczby ,
![]() |
gdzie jest dystrybuantą rozkładu
. Często nie znamy asymptotycznej
wariancji
ale umiemy skonstruować jej zgodny estymator
.
Dla ustalonej ,,małej” dodatniej liczby
łatwo dobrać kwantyl rozkładu normalnego
tak, żeby
. Typowo
i
. Otrzymujemy
![]() |
Ten wzór interpretuje się w praktyce tak: estymator ma błąd nie przekraczający
z prawdopodobieństwem około
. W żargonie statystycznym
jest nazywane asymptotycznym poziomem ufności.
Chociaż opisame powyżej podejście ma swoje słabe strony, to prowadzi do prostego kryterium
porównywania estymatorów. Przypuśćmy, że mamy dwa estymatory i
, oba asymptotycznie normalne, o asymptotycznej wariancji
i
odpowiednio. Przypuśćmy dalej, że dla obliczenia pierwszego z tych estymatorów generujemy
punktów, zaś dla drugiego
. Błędy obu estymatorów na tym mym poziomie istotności są ograniczone przez, odpowiednio,
i
. Przyrównując te wyrażenia do siebie dochodzimy do wniosku, że oba estymatory osiągają podobną dokładność, jeśli
. Liczbę
![]() |
nazywamy względną efektywnością (asymptotyczną). Czasami dobrze jest wybrać
za ,,naturalny punkt odniesienia” estymator CMC i zdefiniować efektywność
estymatora o asymptotycznej wariancji
jako
![]() |
Mówi się też, że jeśli wygenerujemy próbkę punktów i obliczymy
to ,,efektywna liczność próbki” (ESS, czyli effective sample size) jest
. Tyle
bowiem należałoby wygenerować punktów stosując CMC, żeby osiągnąć podobną dokładność.
Asymptotyczna normalność prostego estymatora CMC wynika wprost z Centralnego Twierdzenia Granicznego (CTG). Istotnie, jeśli generujemy niezależnie
o jednakowym rozkładzie
, to
![]() |
gdzie
![]() |
Zupełnie podobnie, dla losowania istotnego w formie (8.1) otrzymujemy asymptotyczną normalność i przy tym
![]() |
Z tego wzorku widać, że estymator może mieć wariancję zero jeśli
. Niestety, żeby obliczyć
potrzebna jest znajomość
współczynnika proporcjonalności, który jest równy…
. Nie wszystko jednak stracone.
Pozostaje ważna reguła heurystyczna:
Gęstość należy tak dobierać, aby jej ,,kształt” był zbliżony do funkcji
.
Dla drugiej wersji losowania istotnego, (8.2), asymptotyczna normalność wynika z następujących rozważań:
![]() |
gdzie
![]() |
Wyrażenie w kwadratowym nawiasie może być ujemne, jeśli jest duża dodatnia korelacja
zmiennych i
. W tej sytuacji estymator IS2 jest lepszy od IS1. Okazuje się
więc rzecz na pozór paradoksalna: dzielenie przez estymator jedynki może poprawić estymator.
Poza tym oba estymatory IS1 i IS2 mogą mieć mniejszą (asymptotyczną) wariancję, niż CMC. Jeśli
efektywność jest większa niż 100%, to używa się czasem określenia ,,estymator superefektywny”.
Jeśli interesuje nas obliczenie wartości oczekiwanej dla wielu różnych funkcji , to
warto wprowadzić uniwersalny, niezależny of
wskaźnik efektywności. Niezłym takim wskaźnikiem jest
.
Istotnie, ,,odchylenie
” pomiędzy gęstością instrumentalną
i docelową
, jest określone poniższym wzorem:
![]() |
Niestety, to ,,odchylenie” nie jest odległością, bo nie spełnia warunku symetrii.
Niemniej widać, że małe wartości świadczą o bliskości
i
.
Zauważmy, jeszcze, że
jest wariancją asymptotyczną estymatora
IS1 dla funkcji
.
Zgodnym estymatorem wielkości jest
![]() |
Dzielenie przez (a nie
) wynika z tradycji. Dzielenie przez
powoduje, że
estymator nie wymaga znajomości czynnika normującego gęstość
.
Jest to w istocie estymator typu IS2. Estymator typu IS1 jest równy po prostu
, ale można go stosować tylko gdy znamy czynnik normujący.
Wprowadźmy teraz bardziej formalnie pojęcie ważonej próbki. Niech będzie
parą zmiennych losowych o wartościach w
. Jeśli rozkład brzegowy
ma gęstość
zaś docelowa gęstość jest postaci
to żądamy aby
![]() |
dla prawie wszystkich . Mówimy wtedy, że próbka jest dopasowana do rozkładu
. Idea jest taka jak w losowaniu istotnym IS2: dopasowanie gwarantuje, że dla każdej funkcji
mamy
. Nie jest konieczne, aby
było deterministyczną funkcją
. Pozostawiamy sobie możliwość generowania wag losowo.
Stała normująca jest na ogół nieznana, a więc
waga
jest określona z dokładnością do proporcjonalności. W poprzednich
podrozdziałach dla podkreślenia opatrywaliśmy takie nieunormowane wagi znaczkiem ,,prim”
ale teraz to pomijamy.
Gdy generujemy ciąg ważonych zmiennych losowych
, to
oczywiście wymaga się aby
, gdzie
stała
musi być jednakowa dla wszystkich
.
Jeśli założy się niezależność par
, to zgodność i asymptotyczną normalność estymatora
danego równaniem (8.2) wykazuje się tak samo jak poprzednio. Asymptotyczna wariancja jest równa
; czynnik
pojawia się dlatego, że operujemy nieunormowanymi wagami.
Jest jednak dużo ciekawych zastosowań, w których
punkty są zależne. Wtedy nawet zgodność estymatora
nie
jest automatyczna. Asymptotyczna normalność może zachodzić z graniczną
wariancją zupełnie inną niż w przypadku niezależnym lub nie zachodzić wcale.
Interesujące jest powiązanie idei eliminacji z ważeniem próbek. Czysta metoda eliminacji pracuje w sytuacji, gdy gęstość ,,instrumentalna” majoryzuje funkcję proporcjonalną do gęstości docelowej, czyli . Jeśli ten warunek nie jest spełniony, to można
naprawić odpowiednio ważąc wylosowaną próbkę. Dokładniej, algorytm ważonej eliminacji
(Rejection Control, w skrócie RC) jest następujący.
repeat |
Gen ![]() |
![]() |
if ![]() |
begin |
![]() |
Gen ![]() |
end; |
until (![]() ![]() |
![]() |
Zauważmy, że prawdopodobieństwo akceptacji w tym algorytmie jest równe
, gdzie
![]() |
Rozkład na wyjściu algorytmu ma zatem gęstość
. Waga
jest więc ,,dopasowana” do rozkładu
w sensie zdefiniowanym w poprzednim podrozdziale, bo
![]() |
Zauważmy jeszcze, jak należy modyfikować wagi w RC, jeśli na wejściu mamy próbkę ważoną,
dopasowaną do .
Jeżeli punkt
, pochodzący z rozkładu
, ma na wejściu wagę
, to
na wyjściu zaakceptowany punkt
ma rozkład
i powinien mieć wagę
. Widać stąd, że
![]() |
gdzie jest napisanym wyżej prawdopodobieństwem akceptacji w RC.
Po angielsku nazywają się Self Avoiding Walks, w skrócie SAW. Niech będzie
-wymiarową kratą całkowitoliczbową.
Mówimy, że ciąg
punktów kraty jest SAW-em jeśli
każde dwa kolejne punkty i
sąsiadują ze sobą, czyli różnią się
o
na dokładnie jednej współrzędnej,
żadne dwa punkty nie zajmują tego samego miejsca, czyli dla
.
Zbiór wszystkich SAW-ów o ogniwach w
oznaczymy
, a
dla
i
ustalonych w skrócie SAW. Przykład
widać na rysunku.
Jak policzyć SAW-y, czyli obliczyć ?
Zainteresujmy się teraz ,,losowo wybranym SAW-em”. Rozumiemy przez to zmienną losową
o rozkładzie jednostajnym w zbiorze
, czyli taką, że
dla każdego
. W skrócie,
.
Niech
oznacza odległość euklidesową końca SAW-a, czyli punktu
, od początku, czyli punktu
. Na przykład dla lańcuszka widocznego na rysunku mamy
.
Można zadać sobie pytanie, jaka jest średni kwadrat takiej odległości, czyli
Jak obliczyć ,
przy założeniu, że
?
Można zastosować prostą metodę MC, czyli eliminację. Niech oznacza zbiór
wszystkich ,,błądzeń”, czyli ciągów
niekonieczne spełniających
warunek ,,
dla
”. Oczywiście
i metoda generowania ,,losowego błądzenia” (z rozkładu
) jest bardzo prosta: kolejno losujemy pojedyncze kroki, wybierając zawsze jedną z
możliwości.
Żeby otrzymać ,,losowy SAW”, stosujemy eliminację. Ten sposób pozwala w zasadzie estymować
(przez uśrednianie długości zaakceptowanych błądzeń) oraz
(przez zanotowanie frakcji akceptowanych błądzeń).
Niestety, metoda jest bardzo nieefektywna, bo dla dużych
prawdopodobieństwo akceptacji szybko zbliża sie do zera.
Metoda ,,wzrostu” zaproponowana przez Rosenbluthów polega na losowaniu kolejnych kroków błądzenia spośród ,,dopuszczalnych punktów”, to znaczy punktów wcześniej nie odwiedzonych.
W każdym kroku, z wyjątkiem pierwszego mamy co najwyżej możliwości. W błądzeniu
widocznym na rysunku kolejne kroki wybieraliśmy spośród:
![]() |
możliwych. Nasz SAW został zatem wylosowany z prawdopodobieństwem
![]() |
Powiedzmy ogólniej, że przy budowaniu SAW-a mamy kolejno
![]() |
możliwości (nie jest przy tym wykluczone, że w pewnym kroku nie mamy żadnej możliwości, ).
Używając terminologii i oznaczeń związanych z losowaniem istotnym powiemy, że
![]() |
jest gęstością instrumentalną dla , gęstość docelowa jest stała, równa
, zatem funkcja podziału jest po prostu liczbą SAW-ów:
.
Wagi przypisujemy zgodnie ze wzorem
(jeśli wygenerowanie
SAW-a się nie udało,
dla pewnego
, to waga jest zero).
Niech teraz
będą niezależnymi błądzeniami losowanymi
metodą Rosenbluthów. Zgodnie z ogólnymi zasadami losowania istotnego,
![]() |
![]() |
Należy przy tym uwzględniać w tych wzorach błądzenia o wadze zero, czyli ,,nieudane SAW-y”.
Rozpatrzymy najprostszy model procesu opisującego straty i przychody w ubezpieczeniowej ,,teorii ryzyka”. Niech będą zmiennymi losowymi oznaczającymi straty netto
(straty
przychody) towarzystwa ubezpieczeniowego w kolejnych okresach czasu. Założymy
(co jest dużym uproszczeniem), że te zmienne są niezależne i mają jednakowy rozkład o gęstości
. Tak zwana ,,nadwyżka ubezpieczyciela” na koniec
-go roku jest równa
![]() |
gdzie jest rezerwą początkową.
Interesuje nas prawdopodobieństwo zdarzenia, polegającego na tym, że
dla pewnego
. Mówimy wtedy (znowu w dużym uproszczeniu) o ,,ruinie ubezpieczyciela”. Wygodnie jest przyjąć następujące oznaczenia i konwencje.
Zmienna losowa
![]() |
oznacza czas oczekiwania na ruinę, przy czym jeśli ruina nigdy nie nastąpi to ten czas uznajemy za nieskończony. Przy takiej umowie, prawdopodobieństwo ruiny możemy zapisać jako .
Wskaźnik
przy symbolu wartości oczekiwanej przypomina, że chodzi tu o ,,oryginalny” proces,
dla którego
.
Obliczenie analitycznie jest możliwe tylko w bardzo specjalnych przykładach. Motody numeryczne istnieją, ale też nie są łatwe. Pokażemy sposób obliczania
metodą Monte Carlo, który stanowi jeden z najpiękniejszych, klasycznych, przykładów losowania istotnego.
Przyjmiemy bardzo rozsądne założenie, że
. Funkcję tworzącą momenty, która odpowiada gęstości
określamy wzorem
![]() |
Założymy, że ta funkcja przyjmuje wartości skończone przynajmniej w pewnym otoczeniu zera
i istnieje takie , że
![]() |
Liczba jest nazywana współczynnikiem dopasowania i odgrywa tu kluczową rolę.
Metoda wykładniczej zamiany miary jest specjalnym przypadkiem losowania istotnego. Zeby określić instrumentalny rozkład prawdopodobieństwa, połóżmy
![]() |
Z definicji współczynnika dopasowania wynika, że jest gęstością prawdopodobieństwa, to znaczy
.
Generuje się ciąg
jednakowo rozłożonych zmiennych losowych o gęstości
.
Dla utworzonego w ten sposób ,,instrumentalnego procesu” używać będziemy symboli
i
. Zauważmy, że
![]() |
ponieważ funkcja tworząca momenty jest wypukła,
i
. Po zamianie miary, proces
na mocy Prawa Wielkich Liczb
zmierza prawie na pewno do minus nieskończoności, a zatem ruina następuje z prawdopodobieństwem jeden,
. Pokażemy, jak wyrazić prstwo ruiny dla oryginalnego procesu w
terminach procesu instrumentalnego. Niech
![]() |
Innymi słowy, zdarzenie zachodzi gdy
.
Mamy zatem
![]() |
Weźmy sumę powyższych równości dla i skorzystajmy z faktu, że
. Dochodzimy do wzoru
![]() |
(8.3) |
Ten fakt jest podstawą algorytmu Monte Carlo:
![]() ![]() |
for ![]() ![]() |
begin |
![]() ![]() |
repeat |
Gen ![]() ![]() |
until ![]() |
![]() ![]() |
end |
![]() |
Algorytm jest prosty i efektywny. Trochę to zadziwiające, że w celu obliczenia prawdopodobieństwa
ruiny generuje się proces dla którego ruina jest pewna. Po chwili zastanowienia można jednak zauważyć,
że wykładnicza zamiana miary realizuje podstawową ideę losowania istotnego: rozkład instrumentalny
,,naśladuje” proces docelowy ograniczony do zdarzenia .
Ciekawe, że wykładnicza zamiana miary nie tylko jest techniką Monte Carlo, ale jest też
techniką dowodzenia twierdzeń! Aby się o tym przekonać, zauważmy, że ,,po drodze”
udowodniliśmy nierówność (wynika to z podstawowego wzoru
(8.3) gdyż
).
Jest to sławna nierówność Lundberga i wcale nie jest ona oczywista.
Treść automatycznie generowana z plików źródłowych LaTeXa za pomocą oprogramowania wykorzystującego LaTeXML.
strona główna | webmaster | o portalu | pomoc
© Wydział Matematyki, Informatyki i
Mechaniki UW, 2009-2010. Niniejsze materiały są udostępnione bezpłatnie na licencji Creative Commons Uznanie autorstwa-Użycie niekomercyjne-Bez utworów zależnych 3.0 Polska.
Projekt współfinansowany przez Unię Europejską w ramach
Europejskiego Funduszu Społecznego.
Projekt współfinansowany przez Ministerstwo Nauki i Szkolnictwa Wyższego i przez Uniwersytet Warszawski.