Z Skrypty dla studentów Ekonofizyki UPGOW
Spis treści |
Próby i schemat Bernoulliego
Próbą Bernoulliego nazywamy dowolne doświadczenie losowe, w którym pytam tylko o dwa możliwe wyniki, będące zdarzeniami przeciwnymi. Prawdopodobieństwo tego, że w schemacie Bernoulliego o n próbach otrzymamy k razy sukces jest jest dane przez rozkład dwumianowy:
$$ p_n(k) = {n \choose k} \cdot p^k \cdot q^{n-k} = \frac{n!}{k! (n-k)!} \cdot p^k \cdot q^{n-k} $$ Przeprowadźmy symulację komputerową takiego procesu dla schematu Bernoulliego składającego się z sześciu prób z prawdopodobieństwem sukcesu \(p=0.6\). Dysponując jednorodnym generatorem liczb losowych z przedziału (0,1), łatwo możemy wygenerować rezultat próby Bernouliego. Korzystamy z faktu, że prawdopodobieństwo P zdarzenia że taka liczba losowa z przedziałuy (0,1) jest wieksza od \(0.6\) wynosi dokładnie \(p=0.6\). Ponieważ mamy sześć prób generujemy jednym poleceniem wektor:
rand(6,1)
a następnie wykonujemy na nim operację logiczną
rand(6,1)<0.6
i w wyniku otrzymujemy wektor zer i jedynek, przy czym jedynka odpowiada sukcesowi z prawdopodopieństwem 0.6 w jednej próbie. Eksperyment powtarzamy wiele razy, więc możemy wykorzystać drugi wskaźnik:
octave:87> rand(6,10)<0.5 ans = 0 1 0 0 1 0 1 0 1 0 0 0 1 1 0 1 0 1 0 0 0 0 1 1 1 0 1 1 0 1 1 1 1 0 1 0 1 1 0 0 0 0 0 1 0 1 0 0 1 1 1 1 0 1 1 1 1 0 0 0
Każda kolumna powyższej macierzy odpowiada jednemu eksperymentowi składającemu się z 6 zdarzeń
Chcemy wyliczyć prawdopodobieństwo z jakim sukces zajdzie dokładnie 3 razy, czyli częstość z jaką kolumna powyższej macierzy będzie zawierała dokładnie trzy jedynki. Można to zrobić sumując wszystkie rzędy tej macierzy wektora:
octave:88> sum( rand(6,10)<0.5 ) ans = 1 2 5 1 2 2 0 4 3 2
a następnie wykonując operację logiczną:
octave:89> sum( rand(6,10)<0.5 )==3 ans = 0 0 0 0 0 0 0 0 1 1
By obliczyć częstość musimy znowu zsumować powyższe "jedynki" i podzielić przez 10
octave:105> sum((sum( rand(6,10)<0.5 )==3))/10 ans = 0.20000
Otrzymany wynik jest średnia z 10 eksperymentów i jego dokładność nie jest duża. Dokładną wartość prawdopodobieństwa możemy oszacować z rozkładu dwumianowego, który jest w systemie Octave/Matlab zdefiniowany jako funkcja binopdf:
octave:106> binopdf(3,6,0.6 ) ans = 0.27648
By otrzymać lepszy wynik możemy powtórzyć nie 10 ale np. 1000000 razy nasz eksperyment:
octave:107> sum((sum( rand(6,1000000)<0.6 )==3))/1000000 ans = 0.27658
i okazuje się zę wynik wysymulowany zgadza się z dokładnym do trzech cyfr znaczących. Na współczesnym komputerze, taka operacja nie powinna zająć wiecej niż jedna sekunda.
Dysponując powyższymi narzędziami można narysować rozkład dwumianowy jako średnia z np. 1000 eksperymentów. Poniższy program ilustruje taką procedurę:
clear all; N=1000 p=0.6 n=61 tic; for k=0:n; eksp(k+1)=sum((sum( rand(n,N)<p )==k))/N; endfor toc plot(0:n,eksp,"*", 0:n,binopdf(0:n,n,p ),"o")
Proces Poissona
Rozważamy przedział liczbowy \([0, T]\). Z przedziału tego wybieram losowo jeden punkt, jedną liczbę. Ponieważ wszystkie liczby są "równo rozłożone", więc prawdopodobieństwo tego, że punkt ten jest w przedziale \((t_1, t_2)\subset [0, T] \) wynosi
$$P(A)= p = \frac{t_2 -t_1}{T}$$
DEFINICJA
Procesem Poissona \(N(t)\) nazywamy proces stochastyczny o następujących wlasnościach:
- Przestrzenią stanów jest zbiór liczb całkowitych nieujemnych, \(X=\{k\}_0^{\infty}\; = \{0, 1, 2, \dots \}\)
- \(N(0) = 0 \; \) (proces startujący z zera)
- \(N(t_2) - N(t_1)\; \) jest liczbą punktów w przedziale \((t_1, t_2)\)
- \(N(t)\) ma stacjonarne i niezależne przyrosty na nieprzekrywających się przedziałach o rozkładzie prawdopodobieństwa
Generowanie procesu Poissona
Ponieważ realizacja procesu Poissona jest dana przez niezależne punkty jednorodnie położone na odcinku \([0, T] \) to można wygenerować tę liczbe korzystając z rozkładu Poissona a następnie wygenerować położenie tych punktów na osi czasu zgodnie z rozkładem jednorodnym. Poniższy kod realizuje ten algorytm. Wykorzystana jest przy tym funkcja poissrnd wbudowana w system, która generuje punkty losowe z rozkładu Poissona \( e^{-\lambda} \; \frac{\lambda ^k}{k!}\)
clear all T=15; mu=1.3; N=poissrnd(T*mu); t=sort(T*rand(1,N)); x=1:(N); t=reshape([t;t],[2*length(t),1])(2:2*length(t)); x=reshape([x;x],[2*length(x),1])(1:(2*length(x)-1)); plot(t,x,"o-")
Proszę zwrócić uwagę na drugą i trzecią linie od końca skryptu, która służy wyłącznie do utworzenia wykresu schodkowego i dodaje "wewnętrzne" rogi schodków:
octave:218> t=[1,2,3]; octave:219> x=[5,6,7]; octave:220> t=reshape([t;t],[2*length(t),1])(2:2*length(t)) t = 1 2 2 3 3 octave:221> x=reshape([x;x],[2*length(x),1])(1:(2*length(x)-1)) x = 5 5 6 6 7 octave:223> [t,x] ans = 1 5 2 5 2 6 3 6 3 7
Inną i częściej spotykaną w praktyce możliwością generacji procesu Poissona jest generacja kolejnych punktów na osi czasu. Kolejne zdarzenia na osi czasu są z definicji procesu Poissona od siebie niezależne. Rozkład czasów miedzy dwoma kolejnymi punktami jest eksponencjalny \( p(t_i) = -\mu \; e^{-\mu t_i}\). Możemy więc wygenerować liczby losowe o rozkładzie eksponencjalnym i ich sekwencja zdefiniuje nam process Poissona.
clear all mu=1.3; N=19; t=exprnd (mu,1,N); for i=2:N t(i)+=t(i-1); endfor x=1:(N); t=reshape([t;t],[2*length(t),1])(2:2*length(t)); x=reshape([x;x],[2*length(x),1])(1:(2*length(x)-1)); plot(t,x,"o-")
Funkcja korelacyjna procesu Poissona
Proces Poissona jest procesem skorelowanym o funkcji korelacji danej wzorem:
$$R(t_2, t_1) = \langle N(t_2) N(t_1)\rangle = \mu^2 \;t_2 \;t_1 + \mu \; \mbox{min}(t_2, t_1) $$
Funkcję korelacyjną można wyznaczyć bezpośrednio z definicji, dysponując wieloma realizacjami procesu losowego.
$$\langle N(t_2) N(t_1)\rangle = \frac{1}{N_r} \sum_{i=1}^{N_r} N(t_2) N(t_1) $$
Dla danych czasów \(t_1,t_2\) obliczamy wartość danej realizacji procesu Poissona
clear all T=8; mu=1.3; Nr=2220; t1=4; t2=6.5; corr=0; for ii=1:Nr N=poissrnd(T*mu); t=[0,sort(T*rand(1,N)) ]; corr+=(max(find(t<t1)) - 1) * ( max(find(t<t2)) -1 ); endfor corr/=Nr; corr mu^2*t1*t2+mu*min(t1,t2)
Proces urodzin i śmierci
W procesie urodzin, liczba osobników nie maleje. W rzeczywistości zachodzą też procesy śmierci, czyli ubytek osobników. Proces ten można uwględnić w relacji w następujący sposób:
$$N(t) = N(0) + \sum_{i} \xi_i \theta (t-t_i)\;$$
gdzie zmienne losowe \(\xi_i=\{1, -1\}\) są niezależne między sobą i są o identycznych rozkładach prawdopodobieństa:
$$P(\xi_i = 1) = p, \; \; \; P(\xi_i = -1) = q, \; \; \; p+q=1$$
clear all mu=1.3; N=60; p=0.5; t=exprnd (mu,1,N); x(1)=0; for i=2:N t(i)+=t(i-1); x(i)=x(i-1)+((rand(1)>0.5)*2-1); endfor t=reshape([t;t],[2*length(t),1])(2:2*length(t)); x=reshape([x;x],[2*length(x),1])(1:(2*length(x)-1)); plot(t,x,"o-")
Błądzenie przypadkowe
Błądzenie przypadkowego to proces losowy zdefiniowany w następujący sposób:
Rozważmy nieskończoną jednowymiarową sieć (łańcuch) o strukturze periodycznej, o okresie \(L\). Węzły sieci oznaczymy liczbami całkowitymi \(\{\dots, -2, -1, 0, 1, 2, \dots \}\). Odległość między węzłami wynosi \(L\). Niech cząstka w chwili początkowej \(t=0\) znajduje się w węźle oznaczonym umownie \(r=0\). Cząstka co pewien ustalony czas \(T\) wykonuje krok albo w prawo (zdarzenie \(A_1\)) albo w lewo (zdarzenie \(A_2\)). Niech prawdopodobieństwo kroku w prawo wynosi \(p\), a kroku w lewo \(q\), czyli\(P(A_1) =p, \; \; \; \; \;\;\;\;\; P(A_2) = q, \; \;\;\;\;\; \; \;p+q=1\),
clear all N=32; p=0.5; T=1.0; x(1)=0; t(1)=0; for i=2:N t(i)=t(i-1)+T; x(i)=x(i-1)+((rand(1)>p)*2-1); endfor t=reshape([t;t],[2*length(t),1])(2:2*length(t)); x=reshape([x;x],[2*length(x),1])(1:(2*length(x)-1)); plot(t,x,"-")