Całkowanie metodą Monte Carlo

Całkowanie metodą Monte Carlo

Plan wykładu:

1. Podstawowa metoda Monte Carlo

2. Metody MC o zwiększonej efektywności a) losowania ważonego

b) zmiennej kontrolnej c) losowania warstwowego d) obniżania krotności całki

3. Przykład zastosowania MC ddo rozwiązania równania dyfuzji i szacowania wartości całek

Przypomnienie wiadomości ze statystyki

i) funkcja gęstości prawodpodobieństwa (fgp) ii) wartość oczekiwana

iii) wariancja (odchylenie standardowe)

Funkcja gęstości prawdopodobieństwa (fgp) ta ma następujące własności

Przy jej pomocy można określić prawdopodobieństwo zdarzenia że zmienna losowa x przyjmie wartość pomiędzy x a x+dx:

Dla danej funkcji gęstości prawdopodobieństwa można określić dystrybuantę rozkładu

która jest funkcją prawostronnie ciągłą i niemalejącą

3 Przykład. Rozkład normalny N(μ,σ) (Gaussa)

fgp:

dystrybuanta:

- wartość oczekiwana zmiennej losowej x

- wariancja zmiennej losowej x

- odchylenie standardowe

Wartość oczekiwana m oraz odchylenie standardowe s są parametrami funkcji gęstości prawdopodobieństwa f(x).

Zazwyczaj chcemy określić wartość oczekiwaną zmiennej losowej z będącej funkcją np. innej zmiennej losowej x, czyli z=z(x):

ale … zazwyczaj g(z) jest nieznana.

Można jendak skorzystać z prawa przenoszenia prawdopodobieństwa

i powyższą całkę przekształcić do postaci

A wariancję

i odchylenie standardowe

liczymy w znany już sposób.

Jeśli ciąg liczb

stanowią zmienne losowe o funkcji gęstości prawdopodobieństwa f(x) to estymatorem wartości oczekiwanej m(z) zmiennej losowej z(x_i) jest średnia z próby

z wariancją

Uwaga:

(N-1) w mianowniku wynika z faktu że średnią wyliczamy z N wartości z(x_i) – znając jej wartość możemy wyliczyć dowolną z(x_i) dysponując N-1

pozostałymi wartościami. Liczba stopni swobody zmniejsza się o 1. W praktyce, dla dużych N jedynkę można pominąć.

Miarą „rozrzutu” zmiennych losowych z_i wokół wartości średniej jest odchylenie standardowe

Ale też jest zmienną losową, ponieważ konstruujemy ją ze zmiennych z_i (każda z nich ma identyczną wariancję) . Jakie jest odchylenie standardowe średniej?

Do jego estymacji możemy użyć s(z)

Wariancję możemy liczyć na dwa sposoby

1 – musimy zapamiętać wszystkie N liczb (odporność na błędy) 2 – na bieżąco liczymy wartości obu sum i zapamiętujemy tylko 2 liczby

(błędy mogą spowodować ujemną wartość różnicy - ale to rzadki przypadek )

Podstawowa metoda całkowania Monte Carlo

Interesuje nas wyznaczenie (a raczej estymacja) wartości oczekiwanej zmiennej losowej

która jest funkcją wektora zmiennych (losowych):

Rozkład prawdopodobieństwa zmiennej losowej z opisuje funkcja gęstości g(z)

a rozkład prawdopodobieństwa wektora x opisuje funkcja gęstości f(x)

9 Przy takich założeniach, zgodnie z CTG dystrybuanta wartości oczekiwanej ma rozkład normalny (niezależnie od rozkładu zmiennej losowej z)

metodę Monte Carlo szacowania wartości całek w wersji podstawowej definiują wzory:

a) wartość całki

Uwaga: x – jest wektorem, którego składowe są niezależnymi zmiennymi losowymi o określonych funkcjach gęstości prawdopodobieństwa

b) błąd oszacowania

Zazwyczaj obszarem całkowania jest określony podzbiór przestrzeni R^M.

W takim przypadku obliczaną całkę trzeba zapisać w nieco zmienionej postaci:

gdzie: jest funkcją przynależności do zbioru

„Kwadratura” Monte Carlo (metoda orzeł-reszka)

Uwagi:

a) w powyższym przypadku zakładamy, że fgp jest stała w obszarze W

b) wydajność metody zależy od stosunku wielkości obszaru V i obszaru W.

11 Przykład. Wyznaczyć pole powierzchni obiektu o nieregularnym kształcie.

Uwagi:

- brzeg obszaru V (tu powierzchnia S) może być opisywany zestawem skomplikowanych równań

- Obszar Ω powienien mieć regularny kształt (np. n-wymiarowy prostopadłościan lub kula) aby w łatwy sposób wygenerować punkty równomiernie go pokrywające

- metoda mało efektywna w d=3 i więcej wymiarach (jest to odpowiednik metody eliminacji von Neumanna)

12 Przykład

Należy obliczyć numerycznie wartość całki

a) metoda trapezów

b) „kwadratura” Monte Carlo

gdzie:

- objętość obszaru regularnego obliczeniowego

(kostka 5-wymiarowa)

13 Rys. Wykres błędu oszacowania wartości całki w zależności od liczby

węzłów (trapezy)/losowań (MC).

C – wartość dokładna całki

I – wartość całki wyznaczona numerycznie

Przykład

Dzielnik napięcia powinien zapewniać tłumienie o wartości 0.5 z dokładnością 2%.

Opory r₁ i r₂ mają rozrzuty produkcyjne które można reprezentować za pomocą niezależnych zmiennych

o funkcjach gęstości prawdopodobieństwa

Wyznaczyć estymatę uzysku produkcyjnego η, czyli średniego odsetka układów sprawnych.

Tłumienie napięciowe dzielnika

15 Tłumienie jest realizacją zmiennej losowej:

Rozkład tej zmiennej opisuje fgp:

zależna od

Warunkiem sprawności układu (jednej z wielu realizowanych możliwości) jest :

Wykorzystujemy metodę MC do estymacji wartości oczekiwanej:

k jest funkcją wektora losowego

dlatego uzysk produkcyjny można wyrazić wzorem na średnią wartość funkcji przynależności:

gdzie:

jest iloczynem ze względu na niezależność zmiennych losowych r₁ i r ₂. Estymatę uzysku można obliczać jako średnią arytmetyczną

gdzie:

są niezależnymi realizacjami wektora losowego r

17 Algorytm wyznaczenia uzysku:

1) Wylosuj parę liczb: r₁ i r₂, zwiększ N o 1

2) Jeśli obliczone k mieści się w obszarze V wówczas zwiększ N_s o 1 3) Uzysk oblicz jako wartość ułamka

Przykład. Wyznaczyć minimalną liczbę N próbek wystarczającą do

wyznaczenia estymaty uzysku z trzysigmowym błędem względnym:

Dla δ =0.1%,1%,10%.

Obliczamy wariancję estymatora:

Błąd względny:

Przekształcając go można otrzymać wyrażenie na minimalną liczbę próbek potrzebną do uzyskania wymaganej dokładności:

19 Metody zwiększania efektywności metody Monte Carlo

Dokładność wyznaczenia całki metodą MC zależy od liczby próbek N oraz wariancji zmiennej losowej:

Wydajność metody można zwiększyć ustalając N i dokonując takiej transformacji funkcji podcałkowej, aby nowa zmienna losowa (funkcja) miała mniejszą wariancję.

Co oznacza wydajność? - to dokładniejszy wynik (mniejszy błąd) uzyskany dla tej samej liczby losowań zmiennej niezależnej.

20 a) Metoda losowania ważonego

Pierwotna postać całki

Pod całkę wprowadzamy nową fgp g_x(x)

i wprowadzamy nową zmienną losową y

wówczas całkę I możemy obliczyć losując wektor x z rozkładem g_x(x) i sumując uzyskane wartości y_n

Zmienna losowa y ma taką samą wartość oczekiwaną jak zmienna losowa G oraz wariancję zależną od nowej fgp g (x).

21 Najmniejszą wartość wariancja osiąga dla:

Jeżeli G(x) jest funkcją nieujemną, wówczas minimalna wariancja estymatora ważonego jest równa 0 - należałoby jednak w takim przypadku znać wartość całki w mianowniku. Zazwyczaj nie jest to możliwe, dlatego funkcję G(x)

zastępuje się inną G₁(x), której całka może być łatwo obliczona. Minimalizacja wariancji w takim przypadku zależy od jakości zastosowanego przybliżenia.

b) Metoda zmiennej kontrolnej Metoda polega na dekompozycji całki

gdzie:

- jest aproksymacją funkcji G(x) umożliwiającą łatwe obliczenie pierwszego wyrazu po prawej stronie (analitycznie lub numerycznie).

Wariancja zmiennej losowej

ma znacznie mniejszą wartość niż pierwotna zmienna G(x).

23 c) Losowanie warstwowe

W metodzie tej obszar całkowania V dzieli się na K rozłącznych podobszarów:

Całkę I oblicza się jako sumę całek w podobszarach

gdzie: k=1,2,3,...,K

Całki I_k można obliczać za pomocą podstawowej wersji metody MC (ze zmodyfikowaną fgp)

Próbki

są realizacjami wektora losowego x o fgp

d) Metoda obniżania krotności całki

Obniżenia krotności całki można dokonać gdy jest możliwa dekompozycja wektora oryginalnych zmiennych losowych:

oraz obszaru

że zachodzi

Zmienna losowa

ma zazwyczaj mniejszą wariancję niż G(x) co pozwala dość łatwo obliczyć całkę zewnętrzną.

Metoda jest skuteczna jeśli potrafimy dość dokładnie i szybko obliczyć całkę

25

Metoda MC wymaga zastosowania generatora liczb pseudolosowych o zadanym rozkładzie gęstości prawdopodobieństwa. Generatory (a raczej ciągi generowanych liczb) muszą spełniać określone warunki (korelacja, okres, fgp itp.).

Zastosowania metody Monte Carlo

Ogólnie: w problemach nauki/techniki, w których zawodzą metody deterministyczne a) sumulacja komputerowa probabilistycznego modelu matematycznego/fizycznego

(np. model Isinga, rozpraszanie neutronów w reaktorze, symulacja zmian wartości akcji w czasie).

b) obliczanie wartości całek wielokrotnych

(obliczanie objętości, momentów bezwładności itp. obiektów o nieregularnym kształcie, całki oddziaływania w chemii/fizyce kwantowej)

c) optymalizacja

(np. minimalizacja czasu oczekiwania pacjenta w kolejce do lekarza, problem komiwojażera,

poszukiwanie minimum energii układu – stan stacjonarny)

d) rozwiązywanie równań różniczkowych zależnych i niezależnych od czasu (np. równanie Poissona metodą błądzenia przypadkowego,

dyfuzja w obszarze o skomplikowanej geometrii)

26

Przykład. Równanie dyfuzji Musimy podać warunek początkowy

Jeśli jako warunek początkowy zadamy deltę Diraca

to dalszą ewolucję u(x,t) opisuje formuła

Załóżmy teraz że u(x,t) ma opisywać zbiór N cząstek, które dla t=0

były skupione w niewielkim obszarze, a dla t>0 dyfundują w prawo i w lewo.

Rozkład gęstości w chwili t=Dt

można wygenerować stochastycznie przesuwając każdą cząstkę o Dx, przy czym Dx otrzymujemy z rozkładu normalnego (np. generator Boxa-Mullera):

27 Rozkład wędrowców (dyfuzja składnika) w wybranych chwilach czasowych

(jeśli wektor przesunięcia wychodzi poza ściankę – to zostaje ona od niej odbita)

Przykład. Obliczanie momentu bezwładności kuli

Obszar kuli:

def. momentu bezwładności:

Wariancja:

g - ciężar właściwy materiału kuli