MKZR:Symulacje procesów losowych dyskretnych

Z Skrypty dla studentów Ekonofizyki UPGOW

(Różnice między wersjami)
(Próby i schemat Bernoulliego)
(Próby i schemat Bernoulliego)
Linia 77: Linia 77:
Na współczesnym komputerze, taka operacja nie powinna zająć wiecej niż jedna sekunda.
Na współczesnym komputerze, taka operacja nie powinna zająć wiecej niż jedna sekunda.
-
[[Plik:poisson_exper.png|thumb|360px|Eksperyment z 61 próbami Bernouliego]]
+
[[Plik:poisson_exper.png|thumb|360px|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. ]]
<source lang="matlab">
<source lang="matlab">
clear all;
clear all;

Wersja z 16:01, 14 kwi 2010

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} \)

Przeprowadzmy symulację komputerową takiego procesu dla schematu Bernoulliego składającego się z szesciu 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

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

Szum dychotomiczny

Proces Poissona

Ruch Browna