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)

Funkcja gęstości prawdopodobieństwa zmiennej losowej x.

Funkcja ta ma następujące własności:

Przy jej pomocy można określić

prawdopodobieństwo zdarzenia że zmienna 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ą.

^

x2[a;b]

f (x) ¸ 0

Z

b a

f (x)dx = 1

P fx · x

ⁱ

· x + dxg = f (x)dx

P fx

¹

· x · x

²

g = Z

b

a

f (x)dx = F (x

₂

) ¡ F (x

¹

)

Przykład. Rozkład normalny (Gaussa)

f (x) = 1 2 p

¼ exp µ

¡ (x ¡ ¹)

²

2¾

²

¶

F (x) = Z

x

¡1

f (x

⁰

)dx

⁰

(3)

- wartość oczekiwana zmiennej losowej x

- wariancja zmiennej losowej x

- odchylenie standardowe

¾

²

(x) = hx

²

i ¡ hxi

²

¾

²

(x) = h[x ¡ hxi]

²

i Z

_b

a

[x ¡ hxi]

²

dx

= Z

b a

£ x

²

¡ 2xhxi + hxi

²

¤

= Z

b a

x

²

f (x)dx ¡ 2hxi Z

b

a

xf (x)dx

+ hxi Z

b

a

f (x)dx

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

W podobny sposób możemy określić wartość

oczekiwaną funkcji, której argumentem jest zmienna losowa x o funkcji gęstości prawdopodobieństwa f(x)

i analogicznie jej wariancję

oraz odchylenie standardowe

hxi = E(x) = ¹(x) = Z

b

a

xf (x)dx

¾

²

(z) = hz

²

i ¡ hzi

²

¾(x) = p

hx

²

i ¡ hxi

²

¾(z) = p

hz

²

i ¡ hzi

²

hzi = ¹(z) =

Z

₁

¡1

zg(z)dz = Z

b

a

z(x)f (x)dx

(4)

Jeśli ciąg liczb

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

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ąć.

fx

ⁱ

g = fx

ⁱ

jn = 1; 2; : : : ; Ng

¹

z = 1 N

X

N i=1

z(x

_i

)

s

²

(z) = 1 N ¡ 1

X

N i=1

[z(x

_i

) ¡ ¹ z]

²

s

²

(z) = 1 N ¡ 1

X

N i=1

(z

_i

¡ ¹ z)

²

= 1

N ¡ 1

µ X

^N

i=1

z

_i²

¡ 1 N

µ X

^N

i=1

z

_i

¶

2

¶

s = p

s

²

(z)

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)

¹ z

¾

²

(¹ z) = ¾

²

Ã 1 N

X

N i=1

z

_i

!

= 1 N

²

X

N i=1

¾

²

(z) = 1

N ¾

²

(z)

s(¹ z) = s(z)

p N

(5)

5 Podstawowa metoda 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)

x x x = [x

1

; x

2

; : : : ; x

m

] z = z(x x x)

Z

₁

¡1

g(z(x x x))dz = 1

Przy takich założeniach, zgodnie z CTG

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

Z

V

f (x x x)dx x x = 1

hzi =

Z

₁

¡1

zg(z)dz = Z

V

z(x x x)f (x x x)dx x x

N

lim

!1

P 8 <

:

j¹ z ¡ hzij

¾(z)p N

· ¸ 9 =

; = 1 p 2¼

Z

¸

¡¸

e

^¡^u2²

du

I = Z

V

z(x x x)f (x x x)dx x x ¼ 1 N

X

N i=1

z(x x x)

¾(I) =

sZ

V

(z ¡ hzi)

²

f (x x x)dx x x

¼ s(z)

p N

(6)

6 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 funkcja gęstości prawdopodobieństwa jest stała w obszarze Ω b) wydajność metody zależy od stosunku wielkości obszaru V i obszaru Ω.

1 11

_V

(x x x)

1 11

_V

(x x x) =

½ 1 dla x x x 2 V 0 dla x x x = 2 V I =

Z

V

z(x x x)f (x x x)dx x x = Z

1 11

_V

(x x x)z(x x x)f (x x x)dx x x

V ½

I = Z

V

z(x x x)dx x x = Z

1

_V

(x x x)z(x x x)dx x x ¼ N

X

N i=1

1

_V

(x x x)z(x x x)

Przykład

Wyznaczyć pole powierzchni obiektu o nieregularnym kształcie.

S = N

X

N i=1

1

_V

= n N S =

Z

V

1d

²

rrr = Z

1

V

d

²

rrr

(7)

7

I

_trap

= h

⁵

X

n

i=0

w

_i

X

n j=0

w

_j

X

n k=0

w

_k

X

n

l=0

w

_l

X

n m=0

w

_m

g(x

_i

; x

_j

; x

_k

; x

_l

; x

_m

)

h = 1

n w

_i;j;k;l;m

=

½

₁

2

i; j; k; l; m = 0; n

1 i; j; k; l; m = 1; 2; : : : ; n ¡ 1

Przykład

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

a) metoda trapezów

b) Kwadratura Monte Carlo

gdzie:

jest wektorem, którego składowe są zmiennymi losowymi

X X X = [X

1

; X

2

; X

3

; X

4

; X

5

] I

_{M C}

= h

⁵

N

X

N i=1

g(X X X) I =

Z

1 0

dx

₁

Z

1

0

dx

₂

Z

1

0

dx

₃

Z

1

0

dx

₄

Z

1

0

dx

₅

g(x

₁

; x

₂

; x

₃

; x

₄

; x

₅

)

Z

₁

0

f (y)dy = h

"

f (y

₀

) + f (y

_n

)

2 +

n

X

¡1 i=1

f (y

_i

)

# y

0

= 0

y

n

= 1

(8)

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

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

(9)

9 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:

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 2 V

V = [0:49; 0:51]

r

₁

r

₂

f

_r₁

(r

₁

) f

_r₂

(r

₂

)

k = r

₁

r

₁

+ r

₂

f

k

(k) f

r₁

(r

1

) f

r₂

(r

2

)

´ =

Z

₁

¡1

11 1

_V

(k)f

_k

(k)dk

rrr = [r

₁

; r

₂

]

^T

(10)

10 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

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 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:

´ = N

_s

N

± = 3¾

_´_^

^

´

^

´ = 1 N

X

N n=1

111

_V

(k(rrr

_n

))

¾

_´²_^

= 1

N (N ¡ 1)

µ X

^N

n=1

¡ 11 1

_V

(k(rrr

_n

)) ¢

2

¡ 1 N

µ X

^N

n=1

1 11

V

(k(rrr

_n

))

¶

2

¶

= ´(1 ^ ¡ ^´) N ¡ 1

´ = Z

R²

111

V

(k(rrr))f

r

(rrr)dr

f

rrr

= f

r₁

(r

1

)f

r₂

(r

2

)

(11)

11

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)dx x x = Z

R^M

111

V

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

z = 1 11

_V

(x x x)G(x x x)

(12)

12 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.

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.

x x x 2 V g

_x_x_x

(x x x)

y = 111

V

G(x x x)f

xxx

=g

xxx

I = 1 N

X

N n=1

y(x x x

_n

)

g

_xxx

(x x x)

g

_x_x_x

(x x x) = 11 1

_V

jG(xxx)jf

^x^x^x

(x x x) R

V

jG(xxx)jf

^x^x^x

(x x x)dx x x

I = E(z) = Z

V

½

G(x x x) f

_x

(x x x)

g

_x

(x x x) g

_x

(x x x)dx x x

¾

(13)

13 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 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

Próbki

są realizacjami wektora losowego x o fgp

G(x ^ x x)

fxxx

^(k)n

jn = 1; 2; : : : ; N

^k

g I =

Z

V

G(x ^ x x)f

_x_x_x

dx x x + Z

V

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

f

_x_x_x

dx x x

V

1

; V

2

; : : : ; V

k

I

_k

= Z

Vk

G(x x x)f

_x_x_x

(x x x)dx x x

= ¹(V

_k

) Z

Vk

G(x x x)f

_x_x_x

=¹(V

_k

)dx x x

¹(V

_k

) = Z

Vk

f

_x_x_x

(x x x)dx x x

I ^

_k

= ¹(V

_k

) N

_k

N_k

X

n=1

111

_V_k

(x x x

^(k)_n

)G(x x x

^(k)_n

)

f

_xxx;k

(x x x) = 111

_V_k

f

_x_x_x

(x x x)

¹(V

_k

)

y = G(x x x) ¡ ^ G(x x x)

(14)

14 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ę 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)

V = V

u

£ V

^v

f

xxx

(x x x) = f

uuu

(u u u)f

vvv

(vvv) u u u 2 V

^u

vvv 2 V

^v

I =

Z

Vu

½ Z

Vv

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

_v

(vvv)dvvv

¾

f

_u

(u u u)du u u

z = Z

Vv

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

_v

(vvv)dvvv

x x x

^T

= [u u u

^T

vvv

^T

]