Z Skrypty dla studentów Ekonofizyki UPGOW
(→Próby i schemat Bernoulliego) |
(→Próby i schemat Bernoulliego) |
||
Linia 78: | Linia 78: | ||
[[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. ]] | [[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. ]] | ||
+ | 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ę: | ||
<source lang="matlab"> | <source lang="matlab"> | ||
clear all; | clear all; |
Wersja z 16:04, 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.
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")