Całkowanie metodą Monte Carlo

(1)

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

(2)

2 Podstawowa metoda Monte Carlo (MC)

Podstawowym zadaniem metody MC jest estymacja wartości oczekiwanej pewnej zmiennej losowej. Wiele problemów można sprowadzić do tego typu zagadnienia.

Założenie: jest zmienną losową o

funkcji gęstości prawdopodobieństwa (fgp) i wartości oczekiwanej:

Jeżeli

są niezależnymi zmiennymi losowymi o f.g.p.

to ich średnia arytmetyczna jest również zmienną losową:

będącą estymatorem wartości oczekiwanej zmiennej losowej

y

¹

_y

= E fyg =

Z

₁

¡1

yf

_y

(y)dy

fy

ⁿ

g = fy

ⁿ

jn = 1; 2; : : : ; Ng

f

_y

(y)

^

¹

_y

= 1 N

X

N n=1

y

_n

y

Zgodnie z założeniami centralnego twierdzenia granicznego przy

średnia arytmetyczna dąży do zmiennej losowej o rozkładzie normalnym z wartością oczekiwaną i wariancją:

Mając do dyspozycji N niezależnych realizacji zmiennej losowej

o pewnej fgp, estymatę wartości oczekiwanej wyznacza się ze wzoru:

a nieobciążoną estymatę wariancji tej zmiennej przy pomocy wyrażenia:

N ! 1

¹

_y

V ar f^ ¹

_y

g = 1 N

²

X

N n=1

V ar fy

ⁿ

g = ¾

_y²

N b

¹

y

fy

ⁿ

g = fy

ⁿ

jn = 1; 2; : : : ; Ng

^

¹

_y

= 1 N

X

N n=1

y

_n

^

¾

_y²

= 1 N ¡ 1

X

N n=1

(y

_n

¡ ^¹

^y

)

²

= 1

N ¡ 1

µ X

^N

n=1

y

²_n

¡ 1 N

µ X

^N

n=1

y

_n

¶

2

¶

¹

_y

(3)

3 Realizacją zmiennej losowej jest estymata

dlatego za miarę dokładności można przyjąć odchylenie standardowe estymatora zmiennej losowe

Obliczając wartość odchylenia standardowego zgodnie z (a) należy dwukrotnie odwoływać się do wartości każdej próbki. W przypadku (b)

odwołanie jest jednokrotne, ale błędy zaokrągleń przy sumowaniu mogą spowodować pojawienie się wartości ujemnych pod pierwiastkiem.

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

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

o funkcji gęstości prawdopodobienstwa

^

¹

y

¹ ^

y

^

¹

y

Podstawowa metoda Monte Carlo Wykorzystujemy przejście

Estymację wartości oczekiwanej zmiennej losowej y można wykonać wg wzoru

oraz wariancję tej zmiennej losowej:

Trzy powyższe wyrażenia są podstawą metody Monte Carlo przybliżonego obliczania całek.

E fyg =

Z

₁

¡1

yf

_y

(x)dy = Z

RRR^M

G(x x x)f

_x

(x x x)dx

^

¹

_y

= 1 N

X

N n=1

G(x x x

_n

)

^

¾

_y²

= 1 N ¡ 1

X

N n=1

µ

G(x x x

_n

) ¡ 1 N

X

N m=1

G(x x x

_m

)

¶

2

y = G(x x x)

x x x = £

x

₁

; : : : ; x

_M

¤

T

f

x

(x x x)

¾

_y_^

= ¾ ^

_y

p N

a

= v u

u t 1

N (N ¡ 1) X

N n=1

(y

_n

¡ ^ ¹

_y

)

²

=

b

= v u u t 1

N ¡ 1 µ 1

N

X

N n=1

y

_n²

¡ µ 1

N

X

N n=1

y

_n

¶

2

¶

(4)

4 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

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:

½ R

^M

r

₁

r

₂

f

_r₁

(r

₁

) f

_r₂

(r

₂

)

k = r

₁

r

₁

+ r

₂

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:

k = r

₁

r

₁

+ r

₂

f

k

(k) f

_r₁

(r

₁

) f

_r₂

(r

₂

)

k 2

= [0:49; 0:51]

rrr = [r

₁

; r

₂

]

^T

´ =

Z

₁

¡1

11 1

(k)f

_k

(k)dk I =

Z

G(x x x)f

x

(x x x)dx x x = Z

R^M

111 (x x x)G(x x x)f

x

(x x x)dx x x

1 11

(x x x) =

½ 1 dla x x x 2

0 dla x x x = 2

1 11

(x x x)

(5)

5 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 realizacjiami wektora losowego r

´ = Z

R²

111 (k(rrr))f

r

(rrr)dr

f

_r_r_r

= f

_r₁

(r

₁

)f

_r₂

(r

₂

)

^

´ = 1 N

X

N n=1

111 (k(rrr

_n

))

rrr

₁

; rrr

₂

rrr

₃

; : : :

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

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:

´ = N

_s

N

± = 3¾

_´_^

^

´

¾

_´²_^

= 1

N (N ¡ 1)

µ X

^N

n=1

¡ 11 1

(k(rrr

_n

)) ¢

2

¡ 1 N

µ X

^N

n=1

111 (k(rrr

_n

))

¶

2

¶

= ´(1 ^ ¡ ^´)

N ¡ 1

(6)

6

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:

Rys. Zależność minimalnej liczby próbek od założnego uzysku

± = 3

s 1 ¡ ^´

(N ¡ 1)^ ´

N = 1 ¡ ^´

^

´

µ 3

±

¶

2

Metody zwiększania efektywności metody Monte Carlo

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

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

I = Z

G(x x x)f

x

(x x x)dx x x = Z

R^M

1 11

(x x x)G(x x x)f

x

(x x x)dx x x

y = 111

(x x x)G(x x x)

(7)

7 a) Metoda losowania ważonego

Zakładamy że jest fgp dodatnio określoną dla

Całkę estymujemy:

Zmienna losowa z ma taką samą wartość oczekiwaną jak zmienna losowa y oraz wariancję zależną od fgp:

Wariancję etsymatora całki można zmniejszyć odpowiednio dobierając fgp.

I = 1 N

X

N n=1

z(x x x

_n

) z = 111

G(x x x)f

xxx

=g

xxx

I = E fzg = Z

½

G(x x x) f

_x_x_x

(x x x)

g

_xxx

(x x x) g

_xxx

(x x x)dx x x

¾ x x x 2

g

_xxx

(x x x) g

_x_x_x

(x x x)

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.

g

_xxx

(x x x) = 11 1

jG(xxx)jf

^x^x^x

(x x x) R

jG(xxx)jf

^x^x^x

(x x x)dx x x

(8)

8 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ą wariancję niż G(x).

c) Losowanie warstwowe

W metodzie tej obszar całkowania 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

Próbki

są realizacjami wektora losowego x o fgp

G(x ^ x x)

z = G(x x x) ¡ ^ G(x x x)

1

;

2

; : : : ;

k

¹(

_k

) = Z

k

f

_x_x_x

(x x x)dx x x

I ^

_k

= ¹(

_k

) N

_k

N_k

X

n=1

11 1

_k

(x x x

^(k)_n

)G(x x x

^(k)_n

)

fxxx

^(k)n

jn = 1; 2; : : : ; N

^k

g

f

_xxx;k

(x x x) = 111

_k

f

_x_x_x

(x x x)

¹(

_k

) I =

Z

G(x ^ x x)f

_xxx

dx x x + Z

h G(x x x) ¡ ^ G(x x x) i

f

_xxx

dx x x

I

_k

= Z

k

G(x x x)f

_x_x_x

(x x x)dx x x

= ¹(

_k

) Z

k

G(x x x)f

_x_x_x

=¹(

_k

)dx x x

(9)

9 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

=

u

£

^v

x x x

^T

= [u u u

^T

vvv

^T

]

f

xxx

(x x x) = f

uuu

(u u u)f

vvv

(vvv) u u u 2

^u

vvv 2

^v

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ę wewnętrzną (analitycznie lub numerycznie).

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

a) sumulacja komputerowa probabilistycznego modelu matematycznego/fizycznego

(kwantowa dyfuzyjna metoda MC).

b) Obliczanie wartości całek wielokrotnych

(obliczanie objętości, momentów bezwładności itp. obiektów o nieregularnym kształcie)

c) Optymalizacja (minimalizacja czasu

oczekiwania pacjenta w kolejce do lekarza) d) Rozwiązywanie równań różniczkowych (rów.

Poissona metodą błądzenia przypadkowego ze stałym lub zmiennym krokiem)

I = Z

u

½ Z

v

G(x x x(u u u; vvv))f

_v

(vvv)dvvv

¾

f

_u

(u u u)du u u

z = Z

v

G(x x x(u u u; vvv))f

_v

(vvv)dvvv