Transcript Monte Carlo

Metody Monte Carlo
Literatura




R.Wit, Metody Monte Carlo Wykłady, Wydawnictwo
Politechniki Częstochowskiej, Częstochowa 2004
N.P. Buslenko, D.I. Golenko,…,Metoda Monte
Carlo, Pwn, Warszawa 1967
S.M. Jermakow, Metoda Monte Carlo i zagadnienia
pokrewne, PWN, Warszawa 1976
D.W. Heermann, Podstawy symulacji
komputerowych w fizyce, WNT, Warszawa 1997
definicja



Nicolas Metropolis, Stanisław Ulam The Monte
Carlo Methods, 1949
Metoda Monte Carlo (MC) – dowolna technika
rozwiązująca problem wykorzystując liczby
losowe
Kluczowa rola – programowy generator liczb
losowych




Ruletka
Kostka do gry
Moneta
Losowania kul z urny
Definicja J.H. Haltona , 1970
MC polega na przedstawieniu
rozwiązania postawionego problemu w
postaci parametru pewnej hipotetycznej
populacji i używaniu losowej sekwencji
liczb do tworzenia próbki tej populacji,
na podstawie której można dokonać
statystycznego oszacowania wartości
badanego parametru
F – rozwiązanie problemu
liczba, zbiór liczb, wartość logiczna – decyzja
Oszacowanie MC wyniku F:
Fˆ  f (r1 , r2 ,...., rn )
{r1,r2,….,rn} - liczby losowe użyte w rozwiązaniu
Zastosowanie MC

Problemy deterministyczne




(wprowadzenie losowości)
Całkowanie
Znajdowanie pól i objętości
Obliczanie liczby π
bezpośrednia symulacja – badanie rzeczywistej
populacji
Problemy probabilistyczne/stochastyczne

Symulacje procesów zależnych od zmiennych
losowych
bezpośrednia symulacja
MC – całkowanie
problem deterministyczny -> problem stochastyczny
f ( x)
x  a; b
b
I

a
f ( x)dx
I
 f ( x) 
ba
ba
I
n
Wartość
średnia
n
 f (x )
i
i 1
Oszacowanie całki: losowy wybór n punktów
xi z <a;b> o rozkładzie równomiernym
Prawo wielkich liczb

prawo wielkich liczb:
xi - zmienne losowe
μ(x) – funkcja gęstości prawdopodobieństwa
ε >0
.

I

  ( x)dx  1

 f ( x) ( x)dx

1
lim P{I   
n 
n
n
 f (x )  I  }  1
i
i 1
niezależnie od wyboru szerokości przedziału wokół wartości
oczekiwanej, prawdopodobieństwo dla dużych n będzie dowolnie
bliskie 1.
Mocne prawo wielkich liczb
1
P{ lim
n  n
n
 f (x )  I}  1
i
i 1
przy wystarczająco dużej próbce losowej
oszacowanie MC zbiega się do poprawnej
wartości
Prawo wielkich liczb

Częstość występowania zdarzenia w n próbach jest
zbieżna do prawdopodobieństwa tego zdarzenia,
gdy n dąży do nieskończoności

zwiększając liczbę doświadczeń opartych na
zdarzeniach losowych, możemy oczekiwać rozkładu
wyników coraz lepiej odpowiadającego rozkładowi
prawdopodobieństw zdarzeń
przykład:
przeprowadzając wielką liczbę rzutów symetryczną
monetą, możemy oczekiwać że stosunek liczby
"wyrzuconych" orłów do liczby wszystkich rzutów będzie
bliski 0,5 (wartości prawdop.); tym większe są na to szanse
im większa jest liczba rzutów
Wyznaczenie liczby π

Pole figury :
x2+y2≤1
k=0; P=4;
for (i=1; i<=n; i++)
{
losuj punkt (x, y) z kwadratu opisanego na
kole: (-1, -1); (-1, 1); (1, 1); (1, -1);
if (x2+y2 <=1) k++;
}
pi=k/n*P;
Prelokacja –
metoda stochastyczna




Siatka 2D
Węzeł zapełniony z prawdopodobieństwem p є [0; 1]
pusty z prawdopodobieństwem (1-p)
Klaster – zbiór zapełnionych węzłów, z których każdy
jest połączony z przynajmniej jednym sąsiednim
(4 sąsiadów) zapełnionym węzłem
Dla p>=pc (próg prelokacji) wystąpi klaster
nieskończony (łączący każdy lub min. 1 brzeg siatki z
brzegiem przeciwnym)
prelokacja
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
określenie progu prelokacji –
symulacja MC


Analityczne wyrażenie na pc – dla siatki 2D i
nieskończeniewymiarowej
Metoda stochastyczna ??

Dla każdego p wielokrotne uruchomienie podprogramu
for(i=0;i<n;i++)
for(j=0;j<n;j++)
R=random(); // losuj liczbę R є[0,1]o rozkładzie równomiernym
if(R<p) L[i][j]=1;
else L[i][j]=0;
sprawdź, czy istnieje klaster prelokacji;

Uśrednienie wyników
problemy symulacji MC

Zależność wyników od rozmiaru sieci




Mniejsza sieć – mniejsza wartość pc
Skończone małe sieci – wartość progu relokacji –
rozmyta, trudna do oszacowania
Wybór odpowiedniej liczby próbek
Generator losowych liczb
Liczby Fibonacciego
fib(0)=1
fib(1)=1
fib(n)=fib(n-1) + fib(n-2) dla n>=2
fib(0)
fib(2)
procedure fib(x)
if(x<2)
then (wynik:=1)
else
(wynik:=fib(x-1)+fib(x-2)
fib(1)
fib(4)
fib(1)
fib(3)
fib(0)
fib(2)
fib(1)
Liczby Fibonacciego


Pierwsza próba modelowania wzrostu populacji
1202 Leonardo z Pizy „Książka o liczydle”




Para królików dojrzała po urodzeniu
Co kwartał – para potomstwa
Króliki nie umierają i nie uciekają
Proces rozmnażania zaczyna się od jednej pary

1  1  5   1  5  
Fk 

,




5  2   2  


k
Wzór Bineta - 1843
k
k 1
złoty podział odcinka
Fk
lim

k  Fk 1
x
1

1 x x
5 1
2
Liczby Fibonacciego




Zastosowanie – generowanie liczb losowych
Nierealistyczne założenia modelu
Nierealistyczne przewidywania w czasie
Użyteczny aparat matematyczny
Równanie logistyczne


Model Verhulsta 1845
Model populacji – liczebność
xn1  r  xn  (1  xn )
gdzie: xn – unormowana liczebność populacji w kroku n (0;1>
r – współczynnik przyjmujący wartości z przedziału <0; 4>,
opisuje liczebność kolonii pewnych osobników unormowaną w
stosunku do maksymalnej liczby elementów, mogących rozwijać
się w zadanym środowisku.
Równanie logistyczne
xn1  r  xn  (1  xn )
y  r  x  (1  x)
r<=1 x*=0 – punkt stabilny rozwiązania; iteracje
kończą się w x*=0
bifurkacje
 1<r<=3
x*=1-1/r - punkt stabilny
 r>3 – dwa punkty stabilne
 4 punkty stabilne
 8…
chaos
Równanie logistyczne – jako generator liczb losowych

równanie logistyczne
r=2
r=3
r=4
y=x
1,2
1
y
0,8
0,6
0,4
0,2
0
0
0,5
1
x
1,5
http://www.kasprzyk.demon.co.uk/www/Logistic.html
Symulowane odprężanie
Przykładem stochastycznej metody Monte-Carlo algorytm symulowanego odprężania

oparty na doświadczeniach, związanych z
oziębianiem materiałów


Powolne ochładzanie może doprowadzić do uzyskania
struktury kryształu, o minimalnej wartości energetycznej
Zgodnie z zasadami mechaniki statystycznej, układ
zmierzając do stanu minimalnej energii, poprzez jej
fluktuację, zależną od temperatury, może zwiększać swoją
energię o Ep z prawdopodobieństwem p, określonym
wzorem:
Symulowane odprężanie
 ΔE p
p  e k T
gdzie:
T – temperatura układu,
k - stała Boltzmanna,
Symulowane odprężanie




Przy wyższych temperaturach układ „przegląda” ogromną liczbę
stanów pokrywającą równomiernie dziedzinę rozwiązań.
Wraz z oziębianiem układu, przejścia do stanu o wyższej energii
są coraz rzadsze i rozwiązania poszukuje się w coraz węższych
obszarach rozwiązań
Akceptacja zmian, które przynoszą „gorsze rozwiązanie”,
umożliwia algorytmowi wyjście z lokalnych minimów.
Dzięki tym własnościom schemat Metropolisa został
wykorzystany w poszukiwaniu minimum globalnego
wielowymiarowych funkcji kryterium przez Kirkpatrick’a. Ze
względu na sposób przeszukiwania, obejmujący duży obszar
rozwiązań, jest on szczególnie efektywny w rozwiązywaniu
zadań kombinatorycznych. Na nim opiera się także działanie
maszyn Boltzmanna.
Schemat Metropolisa
procedure SM;{
i=początkowe rozwiązanie;
T=T0 ;
(wysoka wartość „temperatury początkowej”)
oblicz Ei; (wartość „energii”);
repeat
{
wygeneruj nowe rozwiązanie j (operator sąsiedztwa)
oblicz Ej;
if (Ej < Ei)
then
i=j;
else
if (random[0,1) < exp (- (Ej - Ei)/T))
then i=j;
T=T-T; }
until warunek stopu;}
Problemy statystyczne
symulacji MC




Generowanie ciągu liczb losowych o
zadanym rozkładzie opisującym zmienną
losową
statystyczna analiza modelu
statystyczna analiza wyników
Podstawowe znaczenie - liczby losowe z
rozkładu jednostajnego – można uzyskać
liczby losowe o dowolnym rozkładzie
Generatory liczb losowych





Liczby statystycznie niezależne
O rozkładzie jednostajnym
W ciągu o określonej długości – brak
powtarzania
Duża szybkość generowania
Minimalne wymagania na PAO
Kongruencja – przystawalność
liniowy generator liczb
pseudolosowych




a ≡ b (mod m) => a-b=k*m; kєZ
ni+1 =[ a · ni + b ] (mod m); n,a,b,m – liczby całkowite >=0
 a - mnożnik – duża liczba pierwsza; a mod 8 = 5
m/100 < a < m - √m
w systemie dwójkowym nieoczywisty szablon
 b – stała dodatkowa; b=0 – generator multiplikatywny
b/m ≈ ½ - √3/6
 n0- wartość początkowa – liczba nieparzysta
 m=231-1
ni < m
ni/m - <0;1>
Generator Fibonacciego

ni=(ni-1 + ni-2) mod m - addytywny



okres < 2m-1 – zbyt mało
Prosty
Przesunięta metoda Fibonacciego

ni=(ni-r + ni-s + c) mod m;
r>=s – przesunięcia
c – bit przeniesienia
c= 0 gdy ni-r + ni-s <= m
c=1 gdy ni-r + ni-s > m


Testy sprawdzające losowość
wygenerowanych liczb
Marshall – „fałszywy charakter
losowości wygenerowanej sekwencji
nie wpływa na obliczenia MC
Zadanie optymalizacyjne


min(x,y)F(x,y)
G(x,y)=y – h(x)=0
Szukane – min F(x,y) – funkcja celu;
y=h(x) – ograniczenia
F(x,y)=(x-4)2 + y2
G(x,y)=(x-2)2 + y2 -1 = 0
Szukane: minimum odległości punktu (4, 0) od punktu (x,y) znajdującego się na
okręgu o promieniu 1, o środku w punkcie (2, 0)
Rozwiązanie: (3,0)
(2,0)
(4,0)
Algorytm MC






Wylosuj liczbę losową yi z rozkładem równomiernym
na przedziale (-1; 1)
yi=2(zi-0.5); 0<zi<1
Oblicz xi(1) i xi(2) z G(x,y)=0 /równanie okręgu
Wybierz parę (xi,yi), dla której wartość funkcji celu
F(x,y) jest mniejsza
Wygeneruj yi+1, oblicz xi+1(1) i xi+1(2), wybierz lepsza
parę
Jeśli F(xi,yi) < F(xi+1,yi+1) pozostaw do dalszych
obliczeń parę (xi,yi), w przeciwnym wypadku parę
(xi+1,yi+1)
Przejdź do kroku pierwszego
Szukanie minimum funkcji
algorytm MC
START
N=0
Losuj x0
Wynik = x0
M = f(x0)
Losuj x
N++
N
N >= Nmax
T
T
f(x) < M
M = f(x)
Wynik = x
N
STOP
Zagadnienia optymalizacyjne


Wielomodalna funkcja celu
Deterministyczne procedury





często odnajdują minima lokalne
Są dobrze opisane
Wymagane założenia (ciągłość,..)
Zawarte w bibliotekach numerycznych
Metody wielostartowe


Różne punkty startowe
Deterministyczna procedura
Zagadnienia optymalizacyjne
Strategie poszukiwań
 Punkt startowy -> rozwiązanie; kolejny punkt startowy ->
rozwiązanie
 Jeśli rozwiązania podobne -> rozrzut punktów wokół znalezionego
minimum
 Jeśli różne -> kolejny punkt startowy,…zmiana procedury
optymalizacyjnej,…
 Na początku równomierne rozrzucenie n punktów startowych po
przestrzeni poszukiwań
 n procesów minimalizacyjnych
 Kosztowne
 Na początku równomierne rozrzucenie n punktów startowych po
przestrzeni poszukiwań
 Z każdego punktu – przypadkowy kierunek poszukiwań
(przeszukiwanie jednowymiarowe)
 Efektem – klastry ( zbiory punktów, które prowadzą do wspólnego
minimum)
 Środki ciężkości klastów – punkty startowe dla dokładnych
procedur minimalizacyjnych
podsumowanie

Zakres stosowania MC – bardzo szeroki




Modelowanie procesów rzeczywistych
Eksperymenty komputerowe
Modele deterministyczne i probabilistyczne
Wymagają najlepszych generatorów liczb
losowych

koniec
Gra w życie
automat komórkowy
1011 –
liczba neuronów w organizmie
liczba galaktyk we wszechświecie
6 * 1023
liczba Avogadra
Otoczenie - duże zespoły wzajemnie oddziałujących
elementów zmierzających do stanu równowagi
 Uproszczenia w symulacjach komputerowych



charakter oddziaływań
Ograniczenie oddziaływań do sąsiadów
Gra w życie J.H. Conway



2D; periodyczne warunki brzegowe
Elementy populacji – osobnicy – w węzłach siatki
Reguły przetrwania, śmierci, generowania nowych
osobników
Gra w życie - reguły
Dla każdego elementu populacji:



if ( 2<=liczba_sąsiadów<=3 ) element przeżywa 1
generację
if ( liczba_sąsiadów >=4 ) element umiera
//przeludnienie
if ( liczba_sąsiadów <=1 ) element umiera //izolacja
Dla każdego pustego pola:

if ( liczba_sąsiadów = 3 ) tworzy się nowy element
populacji w tym polu
Gra w życie reguły
Życie i śmierć zachodzą równocześnie
 Przesłanki


Nie istnieje konfiguracja pierwotna, dla której
można udowodnić, że rośnie ona w sposób
nieograniczony
Powinny istnieć konfiguracje pierwotne




prowadzące do wzrostu bez granic
znikające
tworzące stabilną konfigurację
wchodzące w nieskończone oscylacje
Gra w życie reguły zastosowania

Model formowania opinii społecznej






„warszawski”
„wrocławski”
Symulacja rozchodzenia się choroby zakaźnej
Model Isinga
CA + GA
Badanie gęstości upakowania kulek w polach
Model formowania opinii
społecznej







Stany: tak, nie si,j = -1 lub 1
Wartość początkowa – np. 80% populacji – tak; 20% - nie
Opinie rozrzucone losowo
Przeprowadzenie rund dyskusyjnych – wymiana z innymi
członkami populacji (sąsiedzi –odległość emocjonalna), siła
przekonywania
Fi,j – siła przekonywania
Dynamika układu – opisana regułą większości
S=f0,0s0,0(t) + f0,1s0,1(t) + f1,0s1,0(t) + f0,-1s0,-1(t) + f-1, 0s-1,0(t)
0,1
S0,0 = +1, jeśli S>=0
0,0
-1,0
0,-1
1,0
S0,0 = -1, jeśli S<0
Model formowania opinii
społecznej



Inny sposób wyboru sąsiadów , np. losowy, dla zadanego
zasięgu, proporcjonalny do odległości
Szumy
Znając dynamikę zmian pojedynczego osobnika –
obserwacja zmian rozkładu opinii

Tworzenie grup wokół przywódców (osoby o mocnym
wpływie)

Tworzenie wałów ochronnych (słabsi osobnicy za murem
osobników silniejszych)

Rozwój „grup oporu”
Model wrocławski



Bazuje na obserwacji zachowań stadnych
Jedna silnie skorelowana (to samo zdanie na
pewien temat) para potrafi narzucić swoje
zdanie sąsiadom
Jeśli para ma różne zdania – otoczenie nie
zmienia poglądów
Rozchodzenie się choroby
zakaźnej




Obszar N x N; rozmieszczenie osobników w polach
Periodyczne warunki brzegowe lub nie
Szczepienie – wśród losowo wybranej grupy
1 losowo wybrany niezaszczepiony osobnik –
źródłem chorob






V=liczba_osób_zaszczepionych/liczność_populacji
Nz – liczba osób niezaszczepionych
Zk- liczba osób zakażonych (spośród Nz)
I – wskaźnik infekcji I=Zk/Nz < 1
I(V) – funkcja malejąca
Istnieje Vc – wskaźnik infekcji gwałtownie maleje
Model Isinga






Badanie magnetycznych własności ciał
N spinów w węzłach siatki 2D
Spiny oddziałują z sąsiadami i zewnętrznym
polem magnetycznym
Stany spinów: dół -1; góra 1
Cel – minimum energii układu
Metoda Metropolisa
prelokacja



Symulacja pożaru lasu
 Szybkość i kierunek wiatru
 Wilgotność powietrza i poszycia
 Odległość między drzewami
 Istnienie i rozmiary przecinek
 Rozmieszczenie ognisk zapalnych
Pudełko z kulkami przewodnikami i izolatorami umieszczonymi w
dwóch przeciwnych ściankach; po przyłożeniu napięcia – prąd
popłynie jeśli utworzy się prelokujący klaster (przy powolnym
wzroście liczby przewodników - gwałtowne przejście do stanu
przewodnictwa)
Poszukiwanie ropy naftowej, wody – cechy porowatych skał
Algorytm prelokacji węzłowej




Przeglądanie tablicy zajętości wierszami, nadając
etykiety elementom zajętym;



Tablica zajętości
 Wypełniona losowo z prawdopodobieństwem p „1” lub „0”
Wektor pamięci ME
Etykiety – kolejne liczby całkowite (numer klastra)
elementowi, który ma sąsiadów (po lewej stronie i powyżej
(i-1,j); (i,j-1) przyporządkowana jest najniższa z etykiet
sąsiadów (y)
Ustawienie wektora pamięci ME(x)=y; x,y – etykiety; x>y
Uzgadnianie kilkustopniowe – zastąpienie etykiet z
wektora ME ( od największej wartości ) x->y
Etykietowanie elementów sąsiedztwo
1
1
2
3
1
1
1
1
1
3
2
1
1
1
1
1
1
Element
badany
sąsiedzi
1
1
1
prelokacja
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
3
1
1
1
7
2
4
4
7
4
7
2
2
5
3
3
8
3
4
3
3
5
3
3
3
3
8
7
7
1
2
3
ME(4)=3
3 -> 2
1
2
3
4
1
ME(3)=2
ME(6)=4
4 -> 3
1
3
3
3
3
3
3
4
4 4
4
7
6
6 -> 4
1
1
3
3
3
3
3
2
2
2
3
3
1
2
5
2
2
3
8
7
7
2
2
2
2
2
2
2
2
5
2
8
Kropla spadająca na wietrze



Kropla przesuwa się między węzłami sieci pod
wpływem siły ciężkości i wiatru
Cel – określenie średniej wartości dryftu (xk-x0)
Parametry symulacji:





x0, y0 – punkt początkowy
δ x = δy – odległości między współrzędnymi węzłów
p1, p2, p3, p4 - prawdopodobieństwa ruchu w 4 kierunkach
p 1+ p 2+ p 3+ p 4 = 1
Liczba symulacji
Kryterium stopu pojedynczej symulacji – yk= 0
Kropla spadająca na wietrze
Reguły poruszania się po sieci (xi,yi) :
1.
Wylosuj r rzeczywiste є(0,1)
2.
if (r <= p1) (xi,yj) -> (xi+1,yi)
3.
if (r >= p1 && r < p1+p2 ) (xi,yi) -> (xi,yi-1) ↑
4.
if (r >= p1+p2 && r < p1+p2 +p3 ) (xi,yj) -> (xi-1,yi) ←
5.
if (r >= p1+p2+p3 && r < 1 ) (xi,yj) -> (xi,yi+1) ↓
→
Kropla spadająca na wietrze
y0
→
↑
←
↓
yk =0
x0
xk
Ruchy Browna


Brown 1827, Einstein 1905, Smoluchowski
1906
Ruchy małych cząstek zawiesiny w sieci


Wyznaczenie położenia cząstki zawiesiny w chwili
t=0, Δt, 2Δt, 3Δt,….
oraz Δxi, Δyi
Badanie średniej wartości kwadratu przesunięcia
wraz z upływem czasu (liczba kroków)
1
(x) 
N
2

i
(x i ) 2
Ruchy Browna


śr_wartość_kwadratu_x(liczba kroków)zależność liniowa
śr_wartość_kwadratu_y(liczba kroków)zależność liniowa
1200
x
1000
y
800
600
400
200
0
0
500
1000
1500
2000
2500
Ruchy Browna

Nachylenie prostej


Dla gazów - 2D; D – współczynnik dyfuzji
D=kT/(6πηr)




k - stała Boltzmanna
T – temperatura
η współczynnik lepkości
r – promień cząstki
podsumowanie

Zakres stosowania MC – bardzo szeroki




Modelowanie procesów rzeczywistych
Eksperymenty komputerowe
Modele deterministyczne i probabilistyczne
Wymagają najlepszych generatorów liczb
losowych