MKZR:Symulacje procesów losowych dyskretnych

Z Skrypty dla studentów Ekonofizyki UPGOW

(Różnice między wersjami)

Aktualna wersja na dzień 12:51, 17 mar 2011

Spis treści

1 Próby i schemat Bernoulliego
2 Proces Poissona
- 2.1 Generowanie procesu Poissona
- 2.2 Funkcja korelacyjna procesu Poissona
3 Proces urodzin i śmierci
4 Błądzenie przypadkowe

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.

Eksperyment z 61 próbami Bernouliego, punkty połączone linią ciągła sa wynikiem dokładnym a niebieskie punkty są średnią z 1000 eksperymentów.

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

(1)

$Pr\{N(t_2) - N(t_1) =k\} = e^{-\mu (t_2 - t_1)} \; \frac{[\mu (t_2 - t_1)]^k}{k!}$

Generowanie procesu Poissona

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

Realizacja procesu 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$

Realizacja procesu urodzin i śmierci.

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$ ,

Realizacja procesu błądzenia przypadkowego.

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,"-")