WnioskowanieStatystyczne/Zmienne losowe i generatory liczb pseudolosowych: Różnice pomiędzy wersjami
m (→Zadanie: Poprawka dla Python3) |
|||
(Nie pokazano 111 wersji utworzonych przez 6 użytkowników) | |||
Linia 57: | Linia 57: | ||
czyli wartość dystrybuanty w punkcie <math>x</math> jest prawdopodobieństwem, że zmienna losowa przyjmie wartość nie większą niż <math>x</math>. | czyli wartość dystrybuanty w punkcie <math>x</math> jest prawdopodobieństwem, że zmienna losowa przyjmie wartość nie większą niż <math>x</math>. | ||
* Dla zmiennej dyskretnej: | * Dla zmiennej dyskretnej: | ||
− | : <math>F(x)=P(X \le x) = \sum_{x_i | + | : <math>F(x)=P(X \le x) = \sum_{x_i \le x}{P(X = x_i) } = \sum_{x_i \le x}p_i</math> |
* Dla zmiennej losowej ciągłej z funkcją gęstości prawdopodobieństwa <math>f</math>: | * Dla zmiennej losowej ciągłej z funkcją gęstości prawdopodobieństwa <math>f</math>: | ||
: <math>F(x) = P(X \le x)= \int_{-\infty}^x f(s) ds </math> | : <math>F(x) = P(X \le x)= \int_{-\infty}^x f(s) ds </math> | ||
Linia 98: | Linia 98: | ||
* Trzeci moment centralny przydaje się do badania symetrii rozkładu. Przyjmuje on wartość 0 dla zmiennych o rozkładzie symetrycznym, wartości ujemne dla zmiennych o rozkładzie wydłużonym w lewą stronę i wartości dodatnie dla zmiennych o rozkładzie wydłużonym w prawą stronę. | * Trzeci moment centralny przydaje się do badania symetrii rozkładu. Przyjmuje on wartość 0 dla zmiennych o rozkładzie symetrycznym, wartości ujemne dla zmiennych o rozkładzie wydłużonym w lewą stronę i wartości dodatnie dla zmiennych o rozkładzie wydłużonym w prawą stronę. | ||
* Czwarty moment centralny jest przydatny do konstrukcji miary spłaszczenia rozkładu zmiennej losowej '''kurtozy'''. Definiuje się ją następująco: | * Czwarty moment centralny jest przydatny do konstrukcji miary spłaszczenia rozkładu zmiennej losowej '''kurtozy'''. Definiuje się ją następująco: | ||
− | : <math> \mathrm{Kurt} = \frac{\mu_4}{\mu_2^4} - 3 </math> | + | : <math> \mathrm{Kurt} = \frac{\mu_4}{\mu_2^2} - 3 = \frac{\mu_4}{\sigma^4} - 3</math> |
Rozkłady prawdopodobieństwa można podzielić ze względu na wartość kurtozy na rozkłady: | Rozkłady prawdopodobieństwa można podzielić ze względu na wartość kurtozy na rozkłady: | ||
* mezokurtyczne — wartość kurtozy wynosi 0, spłaszczenie rozkładu jest podobne do spłaszczenia rozkładu normalnego (dla którego kurtoza wynosi dokładnie 0) | * mezokurtyczne — wartość kurtozy wynosi 0, spłaszczenie rozkładu jest podobne do spłaszczenia rozkładu normalnego (dla którego kurtoza wynosi dokładnie 0) | ||
Linia 108: | Linia 108: | ||
: <math> | : <math> | ||
− | P_X((-\infty, x_p]) \ | + | P_X((-\infty, x_p]) \leqslant p |
</math> | </math> | ||
Linia 166: | Linia 166: | ||
for i in range(N): | for i in range(N): | ||
x=(a*x+c) % m | x=(a*x+c) % m | ||
− | y[i]= | + | y[i]= x/m |
return y | return y | ||
x=gen(23,100000) | x=gen(23,100000) | ||
Linia 207: | Linia 207: | ||
===Testy jakości=== | ===Testy jakości=== | ||
+ | ==== Zadanie ==== | ||
Dla wyżej wymienionych generatorów, proszę wykonać następujące testy: | Dla wyżej wymienionych generatorów, proszę wykonać następujące testy: | ||
* zobrazować rozkład prawdopodobieństwa (narysować histogram); | * zobrazować rozkład prawdopodobieństwa (narysować histogram); | ||
Linia 223: | Linia 224: | ||
print 'mu= %(m).4f, var= %(v).4f, kurtoza= %(k).4f'% {'m':0.5,'v':1.0/12,'k':-6.0/5} | print 'mu= %(m).4f, var= %(v).4f, kurtoza= %(k).4f'% {'m':0.5,'v':1.0/12,'k':-6.0/5} | ||
print u'wartości z próby:' | print u'wartości z próby:' | ||
− | print 'mu= %(m).4f, var= %(v).4f, kurtoza= %(k).4f'% {'m':mu,'v':var,'k':kurtoza} | + | print 'mu= %(m).4f, var= %(v).4f, kurtoza= %(k).4f'% {'m':mu,'v':var,'k':kurtoza} |
<source lang =python> | <source lang =python> | ||
Linia 231: | Linia 232: | ||
import pylab as py | import pylab as py | ||
def gen(x,N): | def gen(x,N): | ||
− | m=8765132137 | + | m=8765132137 |
− | a=42347 | + | a=42347 |
c=1731 | c=1731 | ||
+ | |||
+ | ## warto porównać wynik dla poniższych parametrów | ||
+ | #m=8191 | ||
+ | #a=101 | ||
+ | #a=107 | ||
+ | |||
y=np.zeros(shape=(N,)) | y=np.zeros(shape=(N,)) | ||
for i in range(N): | for i in range(N): | ||
x=(a*x+c) % m | x=(a*x+c) % m | ||
− | y[i]= | + | y[i] = x/m |
return y | return y | ||
− | N= | + | N=10000 |
x=gen(23,N) | x=gen(23,N) | ||
Linia 253: | Linia 260: | ||
</source> | </source> | ||
− | + | --> | |
Bardziej wyczerpujące baterie testów jakości generatorów można znaleźć na przykład w: | Bardziej wyczerpujące baterie testów jakości generatorów można znaleźć na przykład w: | ||
* http://csrc.nist.gov/groups/ST/toolkit/rng/index.html | * http://csrc.nist.gov/groups/ST/toolkit/rng/index.html | ||
Linia 268: | Linia 275: | ||
Przykładowo: zmienne losowe o rozkładzie normalnym są reprezentowane przez klasę <tt>norm</tt>. | Przykładowo: zmienne losowe o rozkładzie normalnym są reprezentowane przez klasę <tt>norm</tt>. | ||
− | === Funkcje przydatne do pracy ze zmiennymi losowymi dostępne są jako | + | === Funkcje przydatne do pracy ze zmiennymi losowymi dostępne są jako metody === |
Obiekty reprezentujące zmienne losowe podlegające konkretnym rozkładom posiadają następujące metody: | Obiekty reprezentujące zmienne losowe podlegające konkretnym rozkładom posiadają następujące metody: | ||
Linia 281: | Linia 288: | ||
* <tt>isf</tt>: funkcja odwrotna do funkcji przetrwania (ang. Inverse Survival Function (Inverse of SF)); | * <tt>isf</tt>: funkcja odwrotna do funkcji przetrwania (ang. Inverse Survival Function (Inverse of SF)); | ||
− | + | ====Przykład==== | |
+ | Zobaczmy działanie tych metod dla zmiennych losowych z rozkładu normalnego: | ||
<source lang =python> | <source lang =python> | ||
import scipy.stats as st | import scipy.stats as st | ||
Linia 290: | Linia 298: | ||
x = st.norm.rvs(size=10000) | x = st.norm.rvs(size=10000) | ||
py.subplot(2,2,4) | py.subplot(2,2,4) | ||
− | py.hist(x) | + | bins=np.linspace(-3,3,50) |
+ | py.hist(x,bins=bins) | ||
py.title('Histogram') | py.title('Histogram') | ||
py.xlabel(u'wartość zmiennej') | py.xlabel(u'wartość zmiennej') | ||
Linia 335: | Linia 344: | ||
-2.93456088085 | -2.93456088085 | ||
− | === Pomoc do poszczególnych | + | === Pomoc do poszczególnych rozkładów === |
− | Pełną | + | Pełną dokumentację każdego rozkładu mamy dostępną w postaci docstringów. Dla przykładu: |
<source lang = python> | <source lang = python> | ||
help(st.nct) | help(st.nct) | ||
</source> | </source> | ||
− | wypisuje pełną informację z | + | wypisuje pełną informację z docstring o niecentralnym rozkładzie ''t''. |
Podstawowe informacje można uzyskać wypisując treść pola <tt>extradoc</tt> | Podstawowe informacje można uzyskać wypisując treść pola <tt>extradoc</tt> | ||
Linia 349: | Linia 358: | ||
2**df*exp(nc**2/2)*(df+x**2)**(df/2) * gamma(df/2) | 2**df*exp(nc**2/2)*(df+x**2)**(df/2) * gamma(df/2) | ||
for df > 0, nc > 0. | for df > 0, nc > 0. | ||
+ | |||
+ | ==== Zadanie ==== | ||
+ | Wypisz pomoc dla (centralnego) rozkładu ''t'' i sprawdź jakie parametry on przyjmuje. Wybierz dowolnie ich wartości, a następnie posiłkując się [[WnioskowanieStatystyczne/Zmienne_losowe_i_generatory_liczb_pseudolosowych#Przykład|przykładem]]: | ||
+ | # Wygeneruj 10000 liczb z rozkładu ''t'' i narysuj ich histogram. | ||
+ | # Wykreśl gęstość prawdopodobieństwa tego rozkładu. | ||
+ | # Wykreśl dystrybuantę. | ||
+ | # Wykreśl funkcję odwrotną do dystrybuanty. | ||
+ | |||
+ | <!-- | ||
+ | ===== Rozwiązanie ===== | ||
+ | <source lang=python> | ||
+ | import scipy.stats as st | ||
+ | import pylab as py | ||
+ | import numpy as np | ||
+ | |||
+ | N=10000 # ile liczb generujemy | ||
+ | |||
+ | # przykładowe parametry rozkładu | ||
+ | df = 3 # liczba stopni swobody | ||
+ | |||
+ | # generowanie 10000 liczb z rozkładu normalnego | ||
+ | py.subplot(2,2,4) | ||
+ | bins=np.linspace(-3,3,50) | ||
+ | x = st.t.rvs(size=N, df=df) | ||
+ | py.hist(x,bins=bins, histtype="step") #histtype="step" daje nam przezroczysty histogram | ||
+ | x = st.norm.rvs(size=N) | ||
+ | py.hist(x,bins=bins, histtype="step") | ||
+ | |||
+ | py.title('Histogram') | ||
+ | py.xlabel(u'wartość zmiennej') | ||
+ | py.ylabel(u'ilość zliczeń') | ||
+ | |||
+ | # wykreślenie gęstości prawdopodobieństwa | ||
+ | os_x = np.arange(-3,3,0.01) | ||
+ | y=st.t.pdf(os_x, df=df) | ||
+ | py.subplot(2,2,3) | ||
+ | py.plot(os_x,y,'g') | ||
+ | y=st.norm.pdf(os_x) | ||
+ | py.plot(os_x,y) | ||
+ | py.title(u'rozkład gęstości prawdopodobieństwa: pdf') | ||
+ | py.xlabel(u'wartość zmiennej') | ||
+ | py.ylabel(u'gęstość prawdopodobieństwa') | ||
+ | |||
+ | # wykreślenie dystrybuanty | ||
+ | z= st.t.cdf(os_x, df=df) | ||
+ | py.subplot(2,2,1) | ||
+ | py.plot(os_x,z,'g') | ||
+ | z= st.norm.cdf(os_x) | ||
+ | py.plot(os_x,z) | ||
+ | py.title(u'Dystrybuanta: cdf') | ||
+ | |||
+ | py.xlabel(u'wartość zmiennej') | ||
+ | py.ylabel(u'prawdopodobieństwo') | ||
+ | |||
+ | # wykreślenie funkcji odwrotnej do dystrybuanty | ||
+ | os_p = np.arange(0,1,0.01) | ||
+ | f=st.t.ppf(os_p, df=df) | ||
+ | py.subplot(2,2,2) | ||
+ | py.plot(os_p,f,'g') | ||
+ | f=st.norm.ppf(os_p) | ||
+ | py.plot(os_p,f) | ||
+ | py.title(u'funkcja odwrotna do dystrybuanty: ppf') | ||
+ | py.ylabel(u'wartość zmiennej') | ||
+ | py.xlabel(u'prawdopodobieństwo') | ||
+ | |||
+ | py.legend(["Rozkład t-Studenta","Rozkład normalny"]) | ||
+ | py.show() | ||
+ | |||
+ | </source> | ||
+ | --> | ||
+ | |||
+ | === Statystyki i estymatory === | ||
+ | Estymator jest statystyką służącą do szacowania pewnej wielkości (np. parametru) na podstawie próby losowej (np. danych eksperymentalnych). Na wykładzie poznaliśmy m.in. takie statystki jak [[WnioskowanieStatystyczne/Momenty|wartość oczekiwana, wariancja i mediana]], z których ostatnia jest szczególnym przypadkiem [[WnioskowanieStatystyczne/Zmienne_losowe_i_generatory_liczb_pseudolosowych#Kwantyle|kwantylu]]. Moduły numpy oraz scipy.stats umożliwiają łatwe estymowanie niektórych statystyk, a także obliczanie ich rzeczywistych wartości, jeśli znamy parametry rozkładu. | ||
+ | |||
+ | ====Przykład==== | ||
+ | Jako przykład niech posłuży poniższy kod, w którym szukamy wartości estymatorów wartości oczekiwanej oraz wariancji rozkładu normalnego na podstawie próby losowej, a wynik porównujemy z wartościami rzeczywistymi: | ||
+ | <source lang =python> | ||
+ | import numpy as np | ||
+ | import scipy.stats as st | ||
+ | |||
+ | # definiujemy parametry rozkładu, a następnie sam rozkład | ||
+ | m=2 | ||
+ | s=0.5 | ||
+ | rozklad = st.norm(m,s) #definiujemy rozkład normalny | ||
+ | |||
+ | #określamy wielkość próby losowej | ||
+ | N = 10000 | ||
+ | |||
+ | #obliczamy statystyki korzystając z metod obiektu 'rozklad' | ||
+ | mu = rozklad.mean() # rzeczywista wartość oczekiwana naszego rozkładu | ||
+ | var = rozklad.var() # rzeczywiste odchylenie standardow | ||
+ | |||
+ | #losujemy próbę o rozmiarze N | ||
+ | x = rozklad.rvs(size=N) | ||
+ | |||
+ | #estymujemy statystyki | ||
+ | mu_est = np.mean(x) # estymowana wartość oczekiwana | ||
+ | var_est = np.var(x) # estymowana wariancja | ||
+ | |||
+ | # wypisujemy wyniki z dokładnością do 5 cyfr znaczących | ||
+ | print("Wartość oczekiwana rozkładu wynosi {:.5g}, wartość estymowana z próby wynosi {:.5g}".format(mu,mu_est)) | ||
+ | print("Wariancja rozkładu wynosi {:.5g}, wartość estymowana z próby wynosi {:.5g}".format(var,var_est)) | ||
+ | </source> | ||
+ | |||
+ | ==== Zadanie ==== | ||
+ | Posiłkując się powyższym przykładem wygeneruj próby losowe dla następujących rozkładów: | ||
+ | # Rozkład normalny <math>\mu</math> = 100, <math>\sigma</math> = 2. | ||
+ | # Rozkład dwumianowy p = 0.4, n =10. | ||
+ | # Rozkład poissona <math>\lambda</math> = 4. | ||
+ | # Rozkład jednostajny na przedziale [2, 6). | ||
+ | Następnie estymuj na ich podstawie wartości poniższych statystyk i porównaj je z wartościami obliczonymi na podstawie znanych parametrów rozkładu: | ||
+ | # Wartość oczekiwana. | ||
+ | # Odchylenie standardowe. | ||
+ | # Wariancja. | ||
+ | # Mediana. | ||
+ | # Kurtoza. | ||
+ | |||
+ | <!-- | ||
+ | ===== Rozwiązanie ===== | ||
+ | Poniżej znajduje się przykładowe rozwiązanie dla rozkładu dwumianowego. Aby uzyskać rozwiązanie dla pozostałych rozkładów wystarczy dostosować pierwszą część kodu, gdzie definiowany jest rozkład. | ||
+ | <source lang=python> | ||
+ | import numpy as np | ||
+ | import scipy.stats as st | ||
+ | |||
+ | # definiujemy parametry rozkładu, a następnie sam rozkład | ||
+ | n=10 | ||
+ | p=0.4 | ||
+ | rozklad = st.binom(n,p) #definiujemy nasz rozkład | ||
+ | |||
+ | #określamy wielkość próby losowej | ||
+ | N = 10000 | ||
+ | |||
+ | # dalsza część jest jest niezależna od rodzaju rozkładu, który zdefiniowaliśmy powyżej | ||
+ | ########################################################################################### | ||
+ | |||
+ | #obliczamy statystyki korzystając z metod obiektu 'rozklad' | ||
+ | mu = rozklad.mean() # rzeczywista wartość oczekiwana naszego rozkładu | ||
+ | std = rozklad.std() # rzeczywiste odchylenie standardowe | ||
+ | var = rozklad.var() # rzeczywiste odchylenie standardowe | ||
+ | med = rozklad.median() # rzeczywista mediana | ||
+ | kurt = float(rozklad.stats('k'))# rzecz. kurtoza (funkcja zwraca jednoelementową tablicę | ||
+ | # numpy.ndarray, więc trzeba ją ręcznie zrzutować na float) | ||
+ | |||
+ | #losujemy próbę o rozmiarze N | ||
+ | x = rozklad.rvs(size=N) | ||
+ | |||
+ | #estymujemy statystyki | ||
+ | mu_est = np.mean(x) # estymowana wartość oczekiwana | ||
+ | std_est = np.std(x) # estymowane odchylenie standardowe | ||
+ | var_est = np.var(x) # estymowana wariancja | ||
+ | med_est = np.median(x) # estymowana mediana | ||
+ | kurt_est = st.kurtosis(x) # estymowana kurtoza | ||
+ | |||
+ | # wypisujemy wyniki z dokładnością do 5 cyfr znaczących | ||
+ | print("Wartość oczekiwana rozkładu wynosi {:.5g}, wartość estymowana z próby wynosi {:.5g}".format(mu,mu_est)) | ||
+ | print("Odchylenie standardowe rozkładu wynosi {:.5g}, wartość estymowana z próby wynosi {:.5g}".format(std,std_est)) | ||
+ | print("Wariancja rozkładu wynosi {:.5g}, wartość estymowana z próby wynosi {:.5g}".format(var,var_est)) | ||
+ | print("Mediana rozkładu wynosi {:.5g}, wartość estymowana z próby wynosi {:.5g}".format(med,med_est)) | ||
+ | print("Kurtoza rozkładu wynosi {:.5g}, wartość estymowana z próby wynosi {:.5g}".format(kurt,kurt_est)) | ||
+ | </source> | ||
+ | --> | ||
== Generowanie zmiennych losowych podlegającym różnym rozkładom prawdopodobieństwa == | == Generowanie zmiennych losowych podlegającym różnym rozkładom prawdopodobieństwa == | ||
Linia 376: | Linia 546: | ||
ilość sukcesów w n próbach przy prawdopodobieństwie sukcesu p | ilość sukcesów w n próbach przy prawdopodobieństwie sukcesu p | ||
parametry: | parametry: | ||
− | n: | + | n: liczba prób |
p: prawdopodobieństwo sukcesu | p: prawdopodobieństwo sukcesu | ||
− | N: | + | N: liczba liczb do generowania |
funkcja zwraca: | funkcja zwraca: | ||
− | N liczb będącymi | + | N liczb będącymi liczbami sukcesów |
''' | ''' | ||
k = np.zeros(N) | k = np.zeros(N) | ||
Linia 386: | Linia 556: | ||
# losuję n liczb z rozkładu jednostajnego [0,1) i jako sukcesy | # losuję n liczb z rozkładu jednostajnego [0,1) i jako sukcesy | ||
#traktuję wylosowanie liczby nie większej niż p | #traktuję wylosowanie liczby nie większej niż p | ||
− | + | r=np.random.rand(N,n) | |
− | + | k=np.sum(r<=p,axis=1) | |
− | |||
return k | return k | ||
+ | |||
# kod testujący | # kod testujący | ||
Linia 417: | Linia 587: | ||
Odp: ''p'' = 0,127 | Odp: ''p'' = 0,127 | ||
+ | |||
+ | <!-- | ||
=====Rozwiązanie:===== | =====Rozwiązanie:===== | ||
<source lang = python> | <source lang = python> | ||
− | |||
import numpy as np | import numpy as np | ||
import pylab as py | import pylab as py | ||
Linia 425: | Linia 596: | ||
import numpy.random as rnd | import numpy.random as rnd | ||
− | p=2 | + | p=2/3 #prawdopodobieństwo sukcesu w jednej próbie |
− | n=12 # | + | n=12 # liczba prób w jednej akcji reklamowej |
# ze wzoru | # ze wzoru | ||
Linia 441: | Linia 612: | ||
print(p_wz) | print(p_wz) | ||
# z symulacji | # z symulacji | ||
− | N_powt =100000 # | + | N_powt =100000 # liczba powtórzeń symulacji |
k = np.zeros(N_powt) | k = np.zeros(N_powt) | ||
for i in range(N_powt): | for i in range(N_powt): | ||
Linia 449: | Linia 620: | ||
print(p_sym) | print(p_sym) | ||
</source> | </source> | ||
+ | --> | ||
====Zadanie==== | ====Zadanie==== | ||
Na egzaminie testowym jest 30 pytań. Na każde pytanie są podane cztery możliwe odpowiedzi, z czego tylko jedna jest prawdziwa. Za prawidłową odpowiedź student otrzymuje 1 punkt a za fałszywą −0,5 punktu. Próg zaliczenia wynosi 15 punktów. Oblicz prawdopodobieństwo, że udzielając czysto losowych odpowiedzi student zaliczy egzamin. | Na egzaminie testowym jest 30 pytań. Na każde pytanie są podane cztery możliwe odpowiedzi, z czego tylko jedna jest prawdziwa. Za prawidłową odpowiedź student otrzymuje 1 punkt a za fałszywą −0,5 punktu. Próg zaliczenia wynosi 15 punktów. Oblicz prawdopodobieństwo, że udzielając czysto losowych odpowiedzi student zaliczy egzamin. | ||
− | + | <!-- | |
=====Rozwiązanie:===== | =====Rozwiązanie:===== | ||
− | Liczba prawidłowych odpowiedzi w teście może być potraktowana jako liczba sukcesów k w n = 30 próbach. Prawdopodobieństwo sukcesu pojedynczego zdarzenia (prawidłowa odpowiedź na jedno pytanie) wynosi p = 1/4. Ze względu na obecność ujemnych punktów za fałszywe odpowiedzi zmienia się efektywna | + | Liczba prawidłowych odpowiedzi w teście może być potraktowana jako liczba sukcesów k w n = 30 próbach. Prawdopodobieństwo sukcesu pojedynczego zdarzenia (prawidłowa odpowiedź na jedno pytanie) wynosi p = 1/4. Ze względu na obecność ujemnych punktów za fałszywe odpowiedzi zmienia się efektywna liczba punktów wymaganych do zaliczenia zgodnie z równaniem: |
− | + | k-0.5*(30-k) = 15 | |
− | co daje | + | co daje k = 20 |
Zatem problem redukuje się do wyznaczenia prawdopodobieństwa co najmniej 20 sukcesów w 30 próbach, przy czym prawdopodobieństwo sukcesu w 1 próbie wynosi p=1/4: | Zatem problem redukuje się do wyznaczenia prawdopodobieństwa co najmniej 20 sukcesów w 30 próbach, przy czym prawdopodobieństwo sukcesu w 1 próbie wynosi p=1/4: | ||
<source lang=python> | <source lang=python> | ||
Linia 465: | Linia 637: | ||
p = 1./4 #prawdopodobieństwo sukcesu w jednej próbie | p = 1./4 #prawdopodobieństwo sukcesu w jednej próbie | ||
− | n = 30 # | + | n = 30 # liczba prób w jednej akcji reklamowej |
k_kryt = 20 | k_kryt = 20 | ||
# ze wzoru | # ze wzoru | ||
Linia 474: | Linia 646: | ||
N_powt =10000000 # liczba powtórzeń symulacji | N_powt =10000000 # liczba powtórzeń symulacji | ||
r = np.random.random(size=(n,N_powt)) # symulowany test | r = np.random.random(size=(n,N_powt)) # symulowany test | ||
− | k = np.sum(r<=p,0) # | + | k = np.sum(r<=p,0) # liczba sukcesów |
p_sym = sum(k>=20)/N_powt | p_sym = sum(k>=20)/N_powt | ||
print(p_sym) | print(p_sym) | ||
</source> | </source> | ||
+ | --> | ||
=== Rozkład Poissona === | === Rozkład Poissona === | ||
Linia 485: | Linia 658: | ||
P(k)=\frac{\mu^k e^{-\mu}}{k!} | P(k)=\frac{\mu^k e^{-\mu}}{k!} | ||
</math> | </math> | ||
− | Rozkładowi Poissona podlegają zmienne losowe zliczające w jednostce czasu | + | Rozkładowi Poissona podlegają zmienne losowe zliczające w jednostce czasu liczbę zdarzeń o niskim prawdopodobieństwie zajścia. Np. liczba rozpadów promieniotwórczych na jednostkę czasu [[http://brain.fuw.edu.pl/edu/STAT:Przyk%C5%82adowe_rozk%C5%82ady#Rozk.C5.82ad_Poissona więcej o rozkładzie Poissona]] |
====Przykład==== | ====Przykład==== | ||
Linia 509: | Linia 682: | ||
import numpy as np | import numpy as np | ||
import pylab as py | import pylab as py | ||
− | + | import scipy.stats as st | |
+ | |||
def gen_poiss(mi,NT,N): | def gen_poiss(mi,NT,N): | ||
'''funkcja zwracająca zmienne losowe z rozkładu poissona: | '''funkcja zwracająca zmienne losowe z rozkładu poissona: | ||
średnia ilość zdarzeń o niskim prawdopodobieństwie w jednostce czasu | średnia ilość zdarzeń o niskim prawdopodobieństwie w jednostce czasu | ||
parametry: | parametry: | ||
− | mi: średnia | + | mi: średnia liczba zdarzeń |
− | NT: | + | NT: liczba odcinków czasu |
− | N: | + | N: liczba liczb do generowania |
funkcja zwraca: | funkcja zwraca: | ||
N liczb będących zliczeniami zdarzeń | N liczb będących zliczeniami zdarzeń | ||
Linia 525: | Linia 699: | ||
# losuję NT liczb z rozkładu jednostajnego [0,1) i jako zajście zdarzenia | # losuję NT liczb z rozkładu jednostajnego [0,1) i jako zajście zdarzenia | ||
# traktuję wylosowanie liczby nie większej niż p. Ilość zdarzeń to liczba z rozkładu Poissona | # traktuję wylosowanie liczby nie większej niż p. Ilość zdarzeń to liczba z rozkładu Poissona | ||
− | + | r=np.random.rand(N,NT) | |
− | + | k=np.sum(r<=p,axis=1) | |
− | |||
return k | return k | ||
− | + | ||
# kod testujący | # kod testujący | ||
mi=4 | mi=4 | ||
− | N= | + | N=100000 |
NT = 10000 | NT = 10000 | ||
x = gen_poiss(mi,NT,N) | x = gen_poiss(mi,NT,N) | ||
− | bins = | + | bins =np.arange(int(x.max()+1))-0.5 |
− | py.hist(x,bins) | + | py.hist(x,bins,normed=True) |
+ | |||
+ | |||
+ | # porównanie z generatorami scipy.stats | ||
+ | p=mi/NT | ||
+ | k=np.arange(x.max()+1) | ||
+ | B=st.binom.pmf(k,NT,p) | ||
+ | py.plot(k,B,'ro') | ||
+ | P=st.poisson.pmf(k,mi) | ||
+ | py.plot(k,P,'go') | ||
py.show() | py.show() | ||
</source> | </source> | ||
Linia 543: | Linia 725: | ||
Gęstość prawdopodobieństwa [[WnioskowanieStatystyczne/Rozklady-przyklady#Rozk.C5.82ad_Gaussa| rozkładu normalnego]] dana jest wzorem: | Gęstość prawdopodobieństwa [[WnioskowanieStatystyczne/Rozklady-przyklady#Rozk.C5.82ad_Gaussa| rozkładu normalnego]] dana jest wzorem: | ||
− | <math>f(x)=\frac{1}{\sqrt{2 \pi \sigma | + | <math>f(x)=\frac{1}{\sqrt{2 \pi} \sigma} e^{ -\frac{(x-\mu)^2}{2 \sigma^2}}</math> |
gdzie: <math>\mu</math> — średnia, <math>\sigma</math> — odchylenie standardowe (<math>\sigma^2</math> — wariancja). | gdzie: <math>\mu</math> — średnia, <math>\sigma</math> — odchylenie standardowe (<math>\sigma^2</math> — wariancja). | ||
Linia 560: | Linia 742: | ||
* Jakie parametry charakteryzują rozkład do którego zbiegają sumy? Przy każdym z powyższych histogramów wypisz średnią i wariancję z próby (funkcje <tt>np.mean</tt> i <tt>np.var</tt>). | * Jakie parametry charakteryzują rozkład do którego zbiegają sumy? Przy każdym z powyższych histogramów wypisz średnią i wariancję z próby (funkcje <tt>np.mean</tt> i <tt>np.var</tt>). | ||
− | |||
<!-- | <!-- | ||
+ | =====Rozwiązanie:===== | ||
+ | <source lang=python> | ||
# -*- coding: utf-8 -*- | # -*- coding: utf-8 -*- | ||
import numpy as np | import numpy as np | ||
import pylab as py | import pylab as py | ||
+ | import numpy.random as rnd | ||
+ | |||
+ | def gen_CTG_1(N): | ||
+ | '''Funkcja generująca pojedynczą zmienną losową będącą sumą N liczb z rozkładu jednostajnego [0,1) | ||
+ | pomnożoną przez czynnik normujący zapewniający wariancję równą 1. | ||
+ | parametry: | ||
+ | N: ile liczb sumować | ||
+ | funkcja zwraca: | ||
+ | jedną liczbę losową''' | ||
+ | X = np.sqrt(12/N) * (np.sum(rnd.random(N)) - N/2) | ||
+ | return X | ||
− | + | def gen_CTG(N,N_liczb): | |
− | def gen_CTG(N, N_liczb): | + | '''Funkcja zwraca wektor (tablicę) zmiennych losowych będących sumami N liczb z rozkładu jednostajnego [0,1). |
− | ''' | ||
parametry: | parametry: | ||
N: ile liczb sumować | N: ile liczb sumować | ||
N_liczb: ilość liczb do generowania | N_liczb: ilość liczb do generowania | ||
funkcja zwraca: | funkcja zwraca: | ||
− | N_liczb | + | N_liczb losowych''' |
− | + | ||
− | + | X=np.zeros(N_liczb) | |
− | + | ||
− | + | for i in range(N_liczb): | |
− | + | X[i] = gen_CTG_1(N) | |
− | + | ||
− | for i in | + | return X |
− | |||
− | |||
− | return | ||
− | |||
− | |||
# kod testujący | # kod testujący | ||
N_liczb=10000 | N_liczb=10000 | ||
− | + | Nmax=12 | |
− | for | + | |
− | + | for N in range(1,Nmax+1): | |
− | py.subplot(4, | + | X=gen_CTG(N,N_liczb) |
+ | py.subplot(3,4,N) | ||
bins = np.arange(-4,4,0.1) | bins = np.arange(-4,4,0.1) | ||
− | py.hist( | + | py.hist(X,bins=bins,normed=True) |
− | py.title( | + | py.title("N=%i, µ=%.2g, σ²=%.2g"%(N, X.mean(), X.var(ddof=1))) |
+ | |||
py.show() | py.show() | ||
+ | </source> | ||
--> | --> | ||
− | ==== | + | ====Zadanie: transformacja rozkładu normalnego ==== |
Linia 610: | Linia 801: | ||
* narysować histogramy zmiennych przed i po transformacji. | * narysować histogramy zmiennych przed i po transformacji. | ||
* narysować rozkład gęstości prawdopodobieństwa dla rozkładu <math>N(2,9)</math> (skorzystać z funkcji <tt>st.norm.pdf(...)</tt>) . | * narysować rozkład gęstości prawdopodobieństwa dla rozkładu <math>N(2,9)</math> (skorzystać z funkcji <tt>st.norm.pdf(...)</tt>) . | ||
− | |||
<!-- | <!-- | ||
+ | ===== Rozwiązanie ===== | ||
+ | |||
+ | <source lang =python> | ||
+ | # -*- coding: utf-8 -*- | ||
+ | import numpy as np | ||
+ | import pylab as py | ||
import scipy.stats as st | import scipy.stats as st | ||
− | |||
− | |||
+ | Z=st.norm.rvs(loc=0,scale=1,size=100000) #losowanie liczb z rozkładu N(0,1) (czyli Z) | ||
+ | |||
+ | std=9 | ||
+ | mu=2 | ||
+ | X=Z*std+mu #transformacja do σ=9, µ=2 | ||
+ | |||
+ | szerokosc_binu=0.5 | ||
+ | bins=np.arange(-40,40.5,szerokosc_binu) #te same granice binów dla obu histogramów | ||
+ | |||
+ | py.subplot(211) | ||
+ | py.hist(Z,bins=bins,normed=True,label="histogram przed transformacją") | ||
+ | py.title("µ = %.3g σ = %.3g"%(Z.mean(),Z.std())) | ||
+ | py.legend() | ||
− | + | py.subplot(212) | |
+ | py.hist(X,bins=bins,normed=True,label="histogram po transformacji") | ||
− | std = | + | PDF=st.norm.pdf(bins,loc=mu,scale=std) #modelowa gęstość prawdopodobieństwa dla σ=9, µ=2 |
− | + | py.plot(bins+szerokosc_binu/2,PDF,"ro",label="gęstość prawd.") | |
− | X | + | py.title("µ = %.3g σ = %.3g"%(X.mean(),X.std())) |
+ | py.legend() | ||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
py.show() | py.show() | ||
+ | </source> | ||
--> | --> | ||
− | + | ====Transformacja Boxa-Mullera==== | |
+ | To inna popularna metoda metoda generowania liczb losowych o rozkładzie normalnym, na podstawie dwóch wartości zmiennej o rozkładzie jednostajnym na przedziale [0,1) http://en.wikipedia.org/wiki/Box-Muller_transform, ale nie będziemy się nią tu szerzej zajmować. | ||
====Przykład==== | ====Przykład==== | ||
Linia 669: | Linia 873: | ||
py.plot([m_kryt, m_kryt],[0, max(h[0])],'r') | py.plot([m_kryt, m_kryt],[0, max(h[0])],'r') | ||
− | p1= | + | p1=np.sum(srednia<=m_kryt)/N_rep |
print('prawdopodobieństwo uzyskane z symulacji: %(0).4f' %{'0':p1}) | print('prawdopodobieństwo uzyskane z symulacji: %(0).4f' %{'0':p1}) | ||
py.show() | py.show() | ||
Linia 675: | Linia 879: | ||
</source> | </source> | ||
− | ==== | + | ====Zadania==== |
Teraz kilka najprostszych zadań dotyczących rozkładu normalnego: | Teraz kilka najprostszych zadań dotyczących rozkładu normalnego: | ||
* Znajdźmy prawdopodobieństwo, że ''Z'' < −2,47. Proszę zrobić to na dwa sposoby: raz z użyciem wygenerowanych zmiennych z rozkładu normalnego, drugi raz z użyciem dystrybuanty <tt>norm.cdf</tt> | * Znajdźmy prawdopodobieństwo, że ''Z'' < −2,47. Proszę zrobić to na dwa sposoby: raz z użyciem wygenerowanych zmiennych z rozkładu normalnego, drugi raz z użyciem dystrybuanty <tt>norm.cdf</tt> | ||
− | [ | + | [Odp: ''p'' = 0,0068] |
− | |||
<!-- | <!-- | ||
+ | Rozwiązanie: | ||
+ | <source lang =python> | ||
import scipy.stats as st | import scipy.stats as st | ||
Linia 686: | Linia 891: | ||
z_crit = -2.47 | z_crit = -2.47 | ||
Z=st.norm.rvs(loc=0, scale=1, size=N) | Z=st.norm.rvs(loc=0, scale=1, size=N) | ||
− | p= | + | p= sum(Z < z_crit )/N |
p_cdf = st.norm.cdf(z_crit, loc=0, scale=1) | p_cdf = st.norm.cdf(z_crit, loc=0, scale=1) | ||
print('symulowane p: %.4f' % p ) | print('symulowane p: %.4f' % p ) | ||
print('p z dystrybuanty: %.4f' % p_cdf) | print('p z dystrybuanty: %.4f' % p_cdf) | ||
− | + | </source> | |
+ | |||
+ | * Znaleźć prawdopodobieństwo ''P''(|''Z''| < 2) | ||
+ | [Odp: ''p'' = 0,9545] | ||
+ | |||
− | + | Rozwiązanie: | |
− | + | <source lang =python> | |
− | < | ||
import scipy.stats as st | import scipy.stats as st | ||
+ | import numpy as np | ||
+ | |||
+ | N=100000 | ||
+ | z_crit = 2 | ||
+ | Z=st.norm.rvs(loc=0, scale=1, size=N) | ||
+ | p= sum(np.abs(Z) < z_crit )/N | ||
+ | p_cdf = st.norm.cdf( z_crit, loc=0, scale=1) - st.norm.cdf( -z_crit, loc=0, scale=1) | ||
+ | print('symulowane p: %.4f' % p ) | ||
+ | print('p z dystrybuanty: %.4f' % p_cdf) | ||
+ | </source> | ||
− | |||
− | |||
− | |||
* Koncentracja zanieczyszczeń w półprzewodniku używanym do produkcji procesorów podlega rozkładowi normalnemu o średniej 127 cząsteczek na milion i odchyleniu standardowemu 22. Półprzewodnik może zostać użyty jedynie gdy koncentracja zanieczyszczeń spada poniżej 150 cząstek na milion. Jaka proporcja półprzewodników nadaje się do użycia? | * Koncentracja zanieczyszczeń w półprzewodniku używanym do produkcji procesorów podlega rozkładowi normalnemu o średniej 127 cząsteczek na milion i odchyleniu standardowemu 22. Półprzewodnik może zostać użyty jedynie gdy koncentracja zanieczyszczeń spada poniżej 150 cząstek na milion. Jaka proporcja półprzewodników nadaje się do użycia? | ||
− | Prawdopodobieństwo obliczyć korzystając z dystrybuanty rozkładu normalnego oraz z symulacji. [ | + | Prawdopodobieństwo obliczyć korzystając z dystrybuanty rozkładu normalnego oraz z symulacji. |
− | + | [Odp: ''p'' = 0,852] | |
− | < | + | |
+ | Rozwiązanie: | ||
+ | <source lang =python> | ||
import scipy.stats as st | import scipy.stats as st | ||
import pylab as py | import pylab as py | ||
import numpy as np | import numpy as np | ||
− | + | ||
− | |||
− | |||
− | |||
− | |||
− | |||
mu=127 | mu=127 | ||
sig=22 | sig=22 | ||
m_kryt=150 | m_kryt=150 | ||
p=st.norm.cdf(m_kryt, loc = mu, scale = sig ) | p=st.norm.cdf(m_kryt, loc = mu, scale = sig ) | ||
− | print('proporcja odczytana z dystrybuanty rozkładu normalnego: % | + | print('proporcja odczytana z dystrybuanty rozkładu normalnego: %4g' %p) |
− | |||
# symulacja | # symulacja | ||
N_prob=10000 | N_prob=10000 | ||
seria=st.norm.rvs(loc=mu, scale=sig, size=N_prob) | seria=st.norm.rvs(loc=mu, scale=sig, size=N_prob) | ||
− | p = | + | p = sum(seria<=m_kryt)/N_prob |
bins = np.arange(60,200,5) | bins = np.arange(60,200,5) | ||
h=py.hist(seria,bins); | h=py.hist(seria,bins); | ||
− | |||
py.plot([m_kryt, m_kryt],[0, max(h[0])],'r') | py.plot([m_kryt, m_kryt],[0, max(h[0])],'r') | ||
− | p1= | + | p1=np.sum(seria<=m_kryt)/N_prob |
− | print('prawdopodobieństwo uzyskane z symulacji: % | + | print('prawdopodobieństwo uzyskane z symulacji: %.4g' %p1) |
py.show() | py.show() | ||
+ | </source> | ||
--> | --> | ||
Linia 747: | Linia 958: | ||
import numpy as np | import numpy as np | ||
import pylab as py | import pylab as py | ||
− | + | import scipy.stats as st | |
+ | |||
N=3 | N=3 | ||
N_liczb=10000 | N_liczb=10000 | ||
− | + | mu=2 | |
− | + | sigma=1 | |
− | + | ||
− | + | y_wyidealizowane=np.zeros(N_liczb) | |
− | + | m_populacji=np.zeros(N_liczb) | |
− | + | s_populacji=np.zeros(N_liczb) | |
+ | y_estymowane=np.zeros(N_liczb) | ||
+ | m_estymowane=np.zeros(N_liczb) | ||
+ | s_estymowane=np.zeros(N_liczb) | ||
for i in range(N_liczb): | for i in range(N_liczb): | ||
− | x= | + | x=st.norm.rvs(loc=mu,scale=sigma,size=N) |
− | + | m_populacji[i]=mu | |
− | + | s_populacji[i]=sigma/np.sqrt(N) | |
− | |||
− | |||
− | |||
+ | m_estymowane[i]=np.mean(x) | ||
+ | s_estymowane[i]=np.std(x,ddof=1)/np.sqrt(N) | ||
+ | |||
+ | y_wyidealizowane[i]=(m_estymowane[i]-mu)/(s_populacji[i]) | ||
+ | y_estymowane[i]=(m_estymowane[i]-mu)/(s_estymowane[i]) | ||
+ | |||
py.subplot(311) | py.subplot(311) | ||
− | py.plot( | + | py.title(u'średnia (wartość oczekiwana)') |
− | py. | + | py.plot(m_estymowane,label="estymowana") |
− | py. | + | py.plot(m_populacji,'r',label="populacji") |
− | py. | + | py.legend() |
− | py.title('odchylenie standardowe | + | |
− | py. | + | py.subplot(312) |
− | py.plot( | + | py.title('odchylenie standardowe') |
− | py. | + | py.plot(s_estymowane,label="estymowane") |
+ | py.plot(s_populacji,'r',label='populacji') | ||
+ | py.legend() | ||
+ | |||
py.subplot(325) | py.subplot(325) | ||
− | py.hist( | + | py.title(u'rozkład ilorazu Y dla wersji wyidealizowanej') |
− | + | py.hist(y_wyidealizowane,30,(-4,4)) | |
+ | |||
py.subplot(326) | py.subplot(326) | ||
− | py.hist( | + | py.title(u'rozkład ilorazu Y dla wersji estymowanej') |
− | + | py.hist(y_estymowane,30,(-4,4)) | |
+ | py.show() | ||
+ | </source> | ||
+ | ====Zadanie==== | ||
+ | Dział kontroli jakości producenta silników testuje 10 egzemplarzy silnika nowego typu. Uzyskano wartość średnią 220 KM oraz odchylenie standardowe równe 15 KM. Jakie jest prawdopodobieństwo, że klient, który zamówi 100 silników, otrzyma partię, w której średnia będzie mniejsza niż 217 KM? Wykonaj obliczenia oraz symulację. Porównaj wynik z wcześniejszym [[WnioskowanieStatystyczne/Zmienne_losowe_i_generatory_liczb_pseudolosowych#Przykład_5|przykładem]], gdzie znaliśmy odchylenie standardowe z innego źródła (nie było ono estymowane z próby). | ||
+ | |||
+ | <!-- =====Rozwiązanie===== | ||
+ | |||
+ | Tym razem spoglądamy na problem z punktu widzenia działu kontroli jakości, dla którego zarówno średnia, jak i odchylenie standardowe są wartościami estymowanymi z próby. Zatem szukamy: | ||
+ | |||
+ | <math>P(\bar x <217) =P\left( t<\frac{217-\mu}{\frac{s}{ \sqrt{n_{klient}}}}, df=n_{kontrola}-1 \right) = P\left( t<\frac{217-220}{\frac{15}{ \sqrt{100}}}, df=10-1 \right)=</math> <math>P(t<-2, df=9) = 0,0383</math> | ||
+ | |||
+ | <source lang= python> | ||
+ | import scipy.stats as st | ||
+ | import pylab as py | ||
+ | import numpy as np | ||
+ | |||
+ | # ze wzoru | ||
+ | p=st.t.cdf(-2,df=9) | ||
+ | print('prawdopodobieństwo odczytane z dystrybuanty rozkładu t-Studenta: %(0).4f' %{'0':p}) | ||
+ | |||
+ | # symulacja | ||
+ | |||
+ | mu=220 | ||
+ | s=15 | ||
+ | N_prob_kontrola= 10 #liczba silników skontrolowana przez dział kontroli jakości | ||
+ | N_prob_klient =100 #liczba silników kupowanych przez klienta | ||
+ | |||
+ | m_kryt=217 | ||
+ | |||
+ | N_rep=int(1e5) | ||
+ | srednia=np.zeros(N_rep) | ||
+ | for i in range(N_rep): | ||
+ | seria=st.t.rvs(df=N_prob_kontrola-1,loc=mu, scale=s, size=N_prob_klient) | ||
+ | srednia[i]=np.mean(seria) | ||
+ | |||
+ | h=py.hist(srednia,30); | ||
+ | |||
+ | py.plot([m_kryt, m_kryt],[0, max(h[0])],'r') | ||
+ | p1=np.sum(srednia<=m_kryt)/N_rep | ||
+ | print('prawdopodobieństwo uzyskane z symulacji: %(0).4f' %{'0':p1}) | ||
py.show() | py.show() | ||
</source> | </source> | ||
+ | --> | ||
=== Inne rozkłady uzyskiwane przez transformacje === | === Inne rozkłady uzyskiwane przez transformacje === | ||
Linia 789: | Linia 1052: | ||
<!--[http://wazniak.mimuw.edu.pl/index.php?title=Rachunek_prawdopodobieństwa_i_statystyka/Wykład_14:_Komputerowe_metody_statystyki] | <!--[http://wazniak.mimuw.edu.pl/index.php?title=Rachunek_prawdopodobieństwa_i_statystyka/Wykład_14:_Komputerowe_metody_statystyki] | ||
--> | --> | ||
− | Niech <math>\xi</math> będzie zmienną losową o rozkładzie jednostajnym na przedziale <math> | + | Niech <math>\xi</math> będzie zmienną losową o rozkładzie jednostajnym na przedziale <math>[0,1)</math>. Jej dystrybuanta <math>J</math> dana jest więc worem: |
<math> | <math> | ||
Linia 818: | Linia 1081: | ||
====Przykład==== | ====Przykład==== | ||
− | Wygenerować liczby pseudolosowe z rozkładu wykładniczego o | + | Wygenerować liczby pseudolosowe z rozkładu wykładniczego o gęstości prawdopodobieństwa: |
− | <math> | + | <math>f(x) = \lambda e^{-\lambda x} </math> |
Rozwiązanie: | Rozwiązanie: | ||
+ | * Znajdujemy dystrybuantę tego rozkładu: | ||
+ | |||
+ | : <math>F(x) = 1-e^{-\lambda x} </math> | ||
+ | |||
* Znajdujemy funkcję odwrotną do dystrybuanty: | * Znajdujemy funkcję odwrotną do dystrybuanty: | ||
: <math>G(u)=F^{-1}(u) = -\frac{1}{\lambda}\ln(1-u) </math> | : <math>G(u)=F^{-1}(u) = -\frac{1}{\lambda}\ln(1-u) </math> | ||
− | * Generujemy liczby losowe <math>u</math> z przedziału <math> | + | * Generujemy liczby losowe <math>u</math> z przedziału <math>[0,1)</math> i stosujemy transformację <math>G(u)</math> |
* Test generowanego rozkładu robimy przez porównanie dystrybuanty teoretycznej i empirycznej. Proszę zwrócić uwagę na funkcję estymującą dystrybuantę empiryczną | * Test generowanego rozkładu robimy przez porównanie dystrybuanty teoretycznej i empirycznej. Proszę zwrócić uwagę na funkcję estymującą dystrybuantę empiryczną | ||
<source lang =python> | <source lang =python> | ||
+ | # -*- coding: utf-8 -*- | ||
import pylab as py | import pylab as py | ||
import numpy as np | import numpy as np | ||
Linia 852: | Linia 1120: | ||
N_liczb = len(X) | N_liczb = len(X) | ||
D=np.zeros(nx) #wyzerowany wektor na naszą dystrybuantę | D=np.zeros(nx) #wyzerowany wektor na naszą dystrybuantę | ||
− | for i in | + | for i in range(nx): |
− | D[i] = | + | D[i] = sum(X <= os_x[i])/N_liczb |
return (os_x, D) | return (os_x, D) | ||
Linia 866: | Linia 1134: | ||
py.plot(x,F,'r',label = 'dystrybuanta teoretyczna') | py.plot(x,F,'r',label = 'dystrybuanta teoretyczna') | ||
py.plot(os_x, F_emp,'.b', label = 'dystrybuanta empiryczna') | py.plot(os_x, F_emp,'.b', label = 'dystrybuanta empiryczna') | ||
− | py.legend( ) | + | py.legend() |
py.show() | py.show() | ||
</source> | </source> | ||
Linia 874: | Linia 1142: | ||
Zauważmy, że nie zawsze jest łatwo wyznaczyć efektywnie <math>F^{-1}</math> — jest tak na przykład w przypadku rozkładu normalnego. W takich sytuacjach można szukać przybliżonej wartości <math>F^{-1}</math>, rozwiązując numerycznie równanie: | Zauważmy, że nie zawsze jest łatwo wyznaczyć efektywnie <math>F^{-1}</math> — jest tak na przykład w przypadku rozkładu normalnego. W takich sytuacjach można szukać przybliżonej wartości <math>F^{-1}</math>, rozwiązując numerycznie równanie: | ||
<math>F(x)=u</math>, co sprowadza się do znalezienia miejsc zerowych funkcji <math>H(x)=F(x)-u</math>. | <math>F(x)=u</math>, co sprowadza się do znalezienia miejsc zerowych funkcji <math>H(x)=F(x)-u</math>. | ||
+ | |||
+ | ==== Zadanie ==== | ||
+ | Posiłkując się powyższym przykładem napisz funkcję ''MyGen(N)'' dla rozkładu, którego funkcja gęstości prawdopodobieństwa jest dana wzorem: | ||
+ | <math> f(x) = \frac{1}{\pi (1+x^{2})} </math> | ||
+ | Następnie narysuj dla niego dystrybuantę i dystrybuantę empiryczną na przedziale [-10, 10]. | ||
+ | |||
+ | <!-- | ||
+ | <source lang=python> | ||
+ | import numpy as np | ||
+ | import pylab as py | ||
+ | |||
+ | #f(x) = 1/(π(1+x²)) -- f. gęstości prawdopodobieństwa | ||
+ | #zatem F(x) = arctg(x)/π + 1/2 -- dystrybuanta | ||
+ | #zatem F⁻¹(x) = tg(π(x-1/2)) -- f. odwrotna do dystrybuanty | ||
+ | |||
+ | def MyGen(size=1): | ||
+ | x=np.random.random(size) # liczby z generatora jednostajnego na przedziale [0,1) | ||
+ | y=np.tan(np.pi*(x-1/2)) # korzystamy z metody odwracania dystrybuanty | ||
+ | return y | ||
+ | |||
+ | x=MyGen(int(1e6)) | ||
+ | bins=np.linspace(-10,10,101) | ||
+ | py.subplot(211) | ||
+ | py.hist(x,bins,normed=True,histtype='step') | ||
+ | py.plot(bins,1./(np.pi*(1+bins**2)),'r') | ||
+ | py.subplot(212) | ||
+ | py.hist(x,bins,normed=True,histtype='step',cumulative=True) | ||
+ | py.plot(bins,np.arctan(bins)/np.pi+0.5,'r') | ||
+ | py.show() | ||
+ | </source> | ||
+ | --> |
Aktualna wersja na dzień 10:15, 23 cze 2023
Spis treści
- 1 Pojęcia wstępne
- 2 Pseudolosowość
- 3 Zmienne losowe w module scipy.stats
- 4 Generowanie zmiennych losowych podlegającym różnym rozkładom prawdopodobieństwa
Pojęcia wstępne
Zmienna losowa
- Intuicyjnie
- Zmienna losowa to taka zmienna, która przyjmuje pewną konkretną wartość w wyniku przeprowadzenia pomiaru lub eksperymentu. Na przykład zmienną losową jest liczba oczek, która wypada na kostce do gry. Dopóki trzymamy kostkę w ręce to nie wiemy jaką wartość ta zmienna przyjmie, konkretną wartość zmienna ta przyjmuje po rzucie. Jednak w kolejnych rzutach nie wiemy jakie będą jej wartości.
- Bardziej formalnie
- Zmienna losowa to funkcja przypisująca zdarzeniom elementarnym liczby.
- Zupełnie formalnie
Do formalnej definicji potrzebne jest nam pojęcie przestrzeni probabilistycznej: Przestrzeń probabilistyczna to układ trzech elementów [math](\Omega, F,P)[/math], gdzie:
- [math]\Omega[/math] jest pewnym zbiorem, zwanym przestrzenią zdarzeń elementarnych,
- [math]F[/math] jest [math]\sigma[/math]-ciałem podzbiorów zbioru [math]\Omega[/math]. Elementy tego [math]\sigma[/math]-ciała nazywane są zdarzeniami losowymi,
- [math]P: F\rightarrow [0,1][/math] jest miarą probabilistyczną, tzn.
- [math]P(A) \ge 0 [/math] dla każdego [math]A \in F[/math],
- [math]P(\Omega) = 1[/math],
- jeżeli zbiory [math]A_1, A_2,\dots,A_n \in F[/math] są parami rozłączne, to [math]P\left( \bigcup\limits_{i=1}^{\infty} A_i \right ) = \sum_{i=1}^{\infty} P(A_i) [/math].
Niech [math](\Omega, F,P)[/math] będzie przestrzenią probabilistyczną. Zmienną losową nazywamy dowolna funkcję [math]X[/math], określoną na przestrzeni zdarzeń elementarnych [math]\Omega[/math], o wartościach ze zbioru liczb rzeczywistych i mierzalną względem przestrzeni [math](\Omega, F,P)[/math]. Zmienna losowa jest ciągła jeśli jej zbiór wartości jest ciągły, lub dyskretna, jeśli jej zbiór wartości jest dyskretny.
Rozkład prawdopodobieństwa
Rozkład prawdopodobieństwa — (rozkład zmiennej losowej) miara probabilistyczna określona na [math]\sigma[/math]-ciele podzbiorów zbioru wartości zmiennej losowej, pozwalająca przypisywać prawdopodobieństwa zbiorom wartości tej zmiennej, odpowiadającym zdarzeniom losowym.
Dla zmiennej losowej dyskretnej możemy wprowadzić pojęcie funkcji prawdopodobieństwa: [math]P(X=x_i) = p_i[/math]
Przykład: Dla rzutów monetą funkcję prawdopodobieństwa można zobrazować jako tabelkę:
Zdarzenie: | orzeł | reszka | |
---|---|---|---|
zmienna losowa | [math]x_i[/math] | 0 | 1 |
prawdopodobieństwo | [math]p_i[/math] | 1/2 | 1/2 |
Można też przedstawić ją na wykresie:
Dla zmiennej losowej ciągłej zamiast funkcji prawdopodobieństwa wprowadzamy funkcję gęstości prawdopodobieństwa. Jest to funkcja określona na zbiorze liczb rzeczywistych posiadająca następujące własności:
- jest nieujemna: [math]f(x) \ge 0 [/math]
- prawdopodobieństwo tego, że zmienna losowa przyjmie wartość z przedziału [math](a,b][/math]: [math] \int_a^b f(x) dx = P(a\lt X \le b) [/math]
- Prawdopodobieństwo tego, że zmienna losowa przyjmie dowolną wartość: [math] \int_{-\infty}^{+\infty} f(x) dx = P(-\infty \lt X \le +\infty) = 1 [/math]
Dystrybuanta
Dystrybuantą zmiennej losowej [math]X[/math] nazywamy funkcję [math]F(x)[/math] określoną na zbiorze liczb rzeczywistych jako: [math]F(x) = P(X \le x)[/math]
czyli wartość dystrybuanty w punkcie [math]x[/math] jest prawdopodobieństwem, że zmienna losowa przyjmie wartość nie większą niż [math]x[/math].
- Dla zmiennej dyskretnej:
- [math]F(x)=P(X \le x) = \sum_{x_i \le x}{P(X = x_i) } = \sum_{x_i \le x}p_i[/math]
- Dla zmiennej losowej ciągłej z funkcją gęstości prawdopodobieństwa [math]f[/math]:
- [math]F(x) = P(X \le x)= \int_{-\infty}^x f(s) ds [/math]
Zauważmy ścisły związek dystrybuanty z funkcją prawdopodobieństwa w przypadku dyskretnym i funkcją gęstości prawdopodobieństwa w przypadku ciągłym.
Przykład z monetą:
Momenty
Dla zmiennej losowej można obliczyć momenty zwykłe lub centralne.
- Moment zwykły
- rzędu [math]k[/math] zmiennej losowej [math]X[/math] to wartość oczekiwana [math]k[/math]-tej potęgi tej zmiennej. Oblicza się go tak:
- dla dyskretnej zmiennej losowej
- [math]m_k = E(X^k) = \sum_{i} {x_i^k p_i} [/math]
- tu [math]p_i[/math] oznacza prawdopodobieństwo, że [math]X = x_i[/math]
- dla ciągłej zmiennej losowej
- [math] m_k = E(X^k) = \int\limits_{-\infty}^{\infty} {x^k f(x)dx}[/math]
- tu [math]f(x)[/math] jest funkcją gęstości prawdopodobieństwa.
Zauważmy, że pierwszy moment zwykły to średnia.
- Moment centralny
- rzędu [math]k[/math] zmiennej losowej [math]X[/math] to wartość oczekiwana [math]k[/math]-tej potęgi funkcji [math][x_i - E(x_i)][/math]. Zatem możemy go obliczyć tak:
- dla dyskretnej zmiennej losowej
- [math] \mu_k = E\left( \left[ X-E(X) \right] ^k \right) = \sum_{i} {\left[x_i-E(X)\right]^k p_i} [/math]
- dla ciągłej zmiennej losowej
- [math] \mu_k = E\left(\left[X-E(X)\right]^k\right) = \int\limits_{-\infty}^{\infty} { \left[X-E(X)\right]^k f(x)dx} [/math]
Uwaga:
- Drugi moment centralny ma swoją nazwę. Jest to wariancja, często oznaczana przez [math]\sigma^2[/math].
- Trzeci moment centralny przydaje się do badania symetrii rozkładu. Przyjmuje on wartość 0 dla zmiennych o rozkładzie symetrycznym, wartości ujemne dla zmiennych o rozkładzie wydłużonym w lewą stronę i wartości dodatnie dla zmiennych o rozkładzie wydłużonym w prawą stronę.
- Czwarty moment centralny jest przydatny do konstrukcji miary spłaszczenia rozkładu zmiennej losowej kurtozy. Definiuje się ją następująco:
- [math] \mathrm{Kurt} = \frac{\mu_4}{\mu_2^2} - 3 = \frac{\mu_4}{\sigma^4} - 3[/math]
Rozkłady prawdopodobieństwa można podzielić ze względu na wartość kurtozy na rozkłady:
- mezokurtyczne — wartość kurtozy wynosi 0, spłaszczenie rozkładu jest podobne do spłaszczenia rozkładu normalnego (dla którego kurtoza wynosi dokładnie 0)
- leptokurtyczne — kurtoza jest dodatnia, wartości cechy bardziej skoncentrowane niż przy rozkładzie normalnym
- platokurtyczne — kurtoza jest ujemna, wartości cechy mniej skoncentrowane niż przy rozkładzie normalnym
Kwantyle
Kwantylem rzędu p, gdzie [math]0 \leqslant p \leqslant 1[/math], w rozkładzie empirycznym [math] P_{X} [/math] zmiennej losowej [math]X[/math] nazywamy każdą liczbę [math]x_{p}[/math], dla której spełnione są nierówności
- [math] P_X((-\infty, x_p]) \leqslant p [/math]
oraz
- [math] P_X([x_p,\infty)) \geqslant 1-p. [/math]
W szczególności, kwantylem rzędu p jest taka wartość [math]x_{p}[/math] zmiennej losowej, że wartości mniejsze lub równe od [math]x_{p}[/math] są przyjmowane z prawdopodobieństwem co najmniej p, zaś wartości większe lub równe od [math]x_{p}[/math] są przyjmowane z prawdopodobieństwem co najmniej 1−p.
Możemy zatem myśleć o obliczaniu kwantyla tak:
- zaobserwowaliśmy n wartości zmiennej losowej,
- zapiszemy te wartości na liście,
- następnie listę tę sortujemy w porządku rosnącym,
- kwantylem rzędu p jest albo element znajdujący się na liście w pozycji np jeśli np jest całkowite albo średnią z elementów położonych najbliżej np w przeciwnym wypadku.
Nazwy poszczególnych kwantyli
Kwantyl rzędu 1/2 to inaczej mediana (ściślej zależy to od definicji mediany, przy jej obliczaniu z próbki o parzystej liczbie elementów często stosuje się średnią arytmetyczną dwóch środkowych elementów, szczegóły są w artykule).
Kwantyle rzędu 1/4, 2/4, 3/4 są inaczej nazywane kwartylami.
Kwantyle rzędu 1/5, 2/5, 3/5, 4/5 to inaczej kwintyle.
Kwantyle rzędu 1/10, 2/10, ..., 9/10 to inaczej decyle.
Kwantyle rzędu 1/100, 2/100, ..., 99/100 to inaczej percentyle.
Pseudolosowość
W obliczeniach numerycznych często korzystamy z liczb „losowych”. Tak naprawdę sekwencje liczb, których używamy są jedynie pseudolosowe, tzn. są wytwarzane w deterministyczny (algorytmiczny) sposób, ale sekwencja generowanych wartości ma pewne cechy losowości. W idealnym przypadku cechą tą jest nieprzewidywalność: na podstawie obserwacji dotychczasowych wartości sekwencji niemożliwe jest podanie kolejnych wartości. W praktyce oznacza to nie możność ustalenie ziarna lub stanu wewnętrznego generatora na podstawie obserwacji dowolnie długiego ciągu wygenerowanych bitów. Stan generatora przechowywany jest na zmiennych o skończonej precyzji. Jeśli stan jest przechowywany na [math]n[/math] bitach to górną granicą na długość unikalnego ciągu liczb jest [math]2^n[/math]. Zazwyczaj unikalna sekwencja jest jednak krótsza. Wynika stąd, że generatory liczb pseudolosowych mają okres, po którym sekwencja liczb powtarza się. Podstawowe generatory liczb pseudolosowych wytwarzają liczby podlegające rozkładowi jednostajnemu (płaskiemu).
Istnieją generatory „prawdziwych” liczb losowych, które są osobnymi urządzeniami. Liczy przez nie dostarczane są pewnymi parametrami jakichś procesów fizycznych, przebiegających w tych urządzeniach (zależnie od producenta). Do naszych zastosowań w zupełności wystarczą liczby generowane przez algorytmy komputerowe.
Proste generatory liczb pseudolosowych można otrzymać jako funkcje używające operacji dzielenia modulo (generatory kongruencyjne):
[math]x_{n+1}=f(x_n,x_{n-1}, ..., x_{n-k})\mod M[/math].
Zakładamy, że argumentami funkcji [math]f[/math] są liczby całkowite ze zbioru [math]{0, 1, ...,M-1}[/math]. Dla ustalenia uwagi mogą to być generatory liniowe typu:
[math] x_{n+1}= (ax_{n} + c )\mod M[/math]
Przykład generatora liniowego
Dla przykładu zaimplementujmy jeden taki generator:
# -*- coding: utf-8 -*-
import numpy as np
import pylab as py
def gen(x,N=1):
'''Przykład generatora liniowego.
funkcja zwraca liczby pseudolosowe z przedziału [0,1)
x - ziarno,
N - ile liczb zwrócić'''
m=8191
a=101
c=1731
y=np.zeros(shape=(N,))
for i in range(N):
x=(a*x+c) % m
y[i]= x/m
return y
x=gen(23,100000)
py.hist(x)
py.show()
Liczby m, a, c są parametrami generatora. Liczba x, od której startuje generator nazywana jest ziarnem generatora (ang. seed). Funkcja py.hist() zlicza i wyrysowuje histogram (czy wiemy co to jest histogram?).
- Wywołaj funkcję gen dla tego samego ziarna kilka razy zaobserwuj powtarzalność wartości.
W module numpy mamy do dyspozycji generator liczb pseudolosowych oparty na algorytmie Mersenne Twister
# -*- coding: utf-8 -*-
"""
kod demonstrujacy pseudolosowosc
"""
import numpy as np
seed=3
np.random.seed(seed)
x=np.random.random(5)
np.random.seed(seed)
y=np.random.random(5)
print('x:', x)
print('y:', y)
- Zaobserwuj powtarzalność wartości.
- Wygeneruj i wykreśl 200 kolejnych liczb pseudolosowych.
# -*- coding: utf-8 -*-
import numpy as np
import pylab as py
seed=3
np.random.seed(seed)
x=np.random.random(200)
py.plot(x,'.')
py.show()
Testy jakości
Zadanie
Dla wyżej wymienionych generatorów, proszę wykonać następujące testy:
- zobrazować rozkład prawdopodobieństwa (narysować histogram);
- test korelacji na rysunku: sporządzić wykres gdzie na jednej osi są wartości [math]x_n[/math] na drugiej zaś wartości [math]x_{n+1}[/math].
Bardziej wyczerpujące baterie testów jakości generatorów można znaleźć na przykład w:
- http://csrc.nist.gov/groups/ST/toolkit/rng/index.html
- http://www.phy.duke.edu/~rgb/General/rand_rate.php
Zmienne losowe w module scipy.stats
Wstęp
Zmienne losowe w module scipy.stats są zaimplementowane jako dwie klasy: stats.rv_continuous (bazuje na niej ponad 80 typów rozkładów ciągłych zmiennych losowych) i stats.rv_discrete (bazuje na niej 10 typów rozkładów dyskretnych zmiennych losowych). Aktualną listę dostępnych rozkładów można uzyskać np. przy pomocy następującego skryptu:
from scipy import stats
dc= [d for d in dir(stats) if isinstance(getattr(stats,d), stats.rv_continuous)]
print(dc)
Przykładowo: zmienne losowe o rozkładzie normalnym są reprezentowane przez klasę norm.
Funkcje przydatne do pracy ze zmiennymi losowymi dostępne są jako metody
Obiekty reprezentujące zmienne losowe podlegające konkretnym rozkładom posiadają następujące metody:
- rvs: generator zmiennych losowych danego typu;
- cdf: dystrybuanta (ang. Cumulative Distribution Function);
- ppf: funkcja odwrotna do dystrybuanty (ang. Percent Point Function (Inverse of CDF));
- stats: zwraca momenty centralne rozkładu (średnią, wariancję, skośność i kurtozę);
- moment: zwraca niecentralne momenty rozkładu.
- pdf: funkcja gęstości prawdopodobieństwa (ang. Probability Density Function);
- pmf dla rozkładów dyskretnych pdf jest zamienione na funkcję prawdopodobieństwa (ang. probability mass function; masa przez analogię z gęstością masy i masą punktową);
- sf: funkcja przeżycia (ang. Survival Function) (1−CDF);
- isf: funkcja odwrotna do funkcji przetrwania (ang. Inverse Survival Function (Inverse of SF));
Przykład
Zobaczmy działanie tych metod dla zmiennych losowych z rozkładu normalnego:
import scipy.stats as st
import pylab as py
import numpy as np
# generowanie 10000 liczb z rozkładu normalnego
x = st.norm.rvs(size=10000)
py.subplot(2,2,4)
bins=np.linspace(-3,3,50)
py.hist(x,bins=bins)
py.title('Histogram')
py.xlabel(u'wartość zmiennej')
py.ylabel(u'ilość zliczeń')
# wykreślenie gęstości prawdopodobieństwa
os_x = np.arange(-3,3,0.01)
y=st.norm.pdf(os_x)
py.subplot(2,2,3)
py.plot(os_x,y)
py.title(u'rozkład gęstości prawdopodobieństwa: pdf')
py.xlabel(u'wartość zmiennej')
py.ylabel(u'gęstość prawdopodobieństwa')
# wykreślenie dystrybuanty
z= st.norm.cdf(os_x)
py.subplot(2,2,1)
py.plot(os_x,z)
py.title(u'Dystrybuanta: cdf')
py.xlabel(u'wartość zmiennej')
py.ylabel(u'prawdopodobieństwo')
# wykreślenie funkcji odwrotnej do dystrybuanty
os_p = np.arange(0,1,0.01)
f=st.norm.ppf(os_p)
py.subplot(2,2,2)
py.plot(os_p,f)
py.title(u'funkcja odwrotna do dystrybuanty: ppf')
py.ylabel(u'wartość zmiennej')
py.xlabel(u'prawdopodobieństwo')
py.show()
Zmienne losowe mogą być używane na jeden z dwóch sposobów: można podawać wszystkie parametry opisujące rozkład w każdym wywołaniu metody, albo wytworzyć obiekt reprezentujący rozkład o konkretnych parametrach (w dokumentacji scipy jest to nazywane zamrażaniem parametrów rozkładu, ang. freezing) Jako przykład zobaczmy funkcję odwrotną do dystrybuanty rozkładu normalnego N(2, 9):
>>> import scipy.stats as st >>> print(st.norm.ppf(0.05, loc=2, scale=3)) -2.93456088085 >>> my_norm = st.norm(loc=2,scale=3) >>> print(my_norm.ppf(0.05)) -2.93456088085
Pomoc do poszczególnych rozkładów
Pełną dokumentację każdego rozkładu mamy dostępną w postaci docstringów. Dla przykładu:
help(st.nct)
wypisuje pełną informację z docstring o niecentralnym rozkładzie t. Podstawowe informacje można uzyskać wypisując treść pola extradoc
>>> print(st.nct.extradoc) Non-central Student T distribution df**(df/2) * gamma(df+1) nct.pdf(x,df,nc) = -------------------------------------------------- 2**df*exp(nc**2/2)*(df+x**2)**(df/2) * gamma(df/2) for df > 0, nc > 0.
Zadanie
Wypisz pomoc dla (centralnego) rozkładu t i sprawdź jakie parametry on przyjmuje. Wybierz dowolnie ich wartości, a następnie posiłkując się przykładem:
- Wygeneruj 10000 liczb z rozkładu t i narysuj ich histogram.
- Wykreśl gęstość prawdopodobieństwa tego rozkładu.
- Wykreśl dystrybuantę.
- Wykreśl funkcję odwrotną do dystrybuanty.
Statystyki i estymatory
Estymator jest statystyką służącą do szacowania pewnej wielkości (np. parametru) na podstawie próby losowej (np. danych eksperymentalnych). Na wykładzie poznaliśmy m.in. takie statystki jak wartość oczekiwana, wariancja i mediana, z których ostatnia jest szczególnym przypadkiem kwantylu. Moduły numpy oraz scipy.stats umożliwiają łatwe estymowanie niektórych statystyk, a także obliczanie ich rzeczywistych wartości, jeśli znamy parametry rozkładu.
Przykład
Jako przykład niech posłuży poniższy kod, w którym szukamy wartości estymatorów wartości oczekiwanej oraz wariancji rozkładu normalnego na podstawie próby losowej, a wynik porównujemy z wartościami rzeczywistymi:
import numpy as np
import scipy.stats as st
# definiujemy parametry rozkładu, a następnie sam rozkład
m=2
s=0.5
rozklad = st.norm(m,s) #definiujemy rozkład normalny
#określamy wielkość próby losowej
N = 10000
#obliczamy statystyki korzystając z metod obiektu 'rozklad'
mu = rozklad.mean() # rzeczywista wartość oczekiwana naszego rozkładu
var = rozklad.var() # rzeczywiste odchylenie standardow
#losujemy próbę o rozmiarze N
x = rozklad.rvs(size=N)
#estymujemy statystyki
mu_est = np.mean(x) # estymowana wartość oczekiwana
var_est = np.var(x) # estymowana wariancja
# wypisujemy wyniki z dokładnością do 5 cyfr znaczących
print("Wartość oczekiwana rozkładu wynosi {:.5g}, wartość estymowana z próby wynosi {:.5g}".format(mu,mu_est))
print("Wariancja rozkładu wynosi {:.5g}, wartość estymowana z próby wynosi {:.5g}".format(var,var_est))
Zadanie
Posiłkując się powyższym przykładem wygeneruj próby losowe dla następujących rozkładów:
- Rozkład normalny [math]\mu[/math] = 100, [math]\sigma[/math] = 2.
- Rozkład dwumianowy p = 0.4, n =10.
- Rozkład poissona [math]\lambda[/math] = 4.
- Rozkład jednostajny na przedziale [2, 6).
Następnie estymuj na ich podstawie wartości poniższych statystyk i porównaj je z wartościami obliczonymi na podstawie znanych parametrów rozkładu:
- Wartość oczekiwana.
- Odchylenie standardowe.
- Wariancja.
- Mediana.
- Kurtoza.
Generowanie zmiennych losowych podlegającym różnym rozkładom prawdopodobieństwa
Podstawowe typy generatorów liczb losowych wytwarzają liczby losowe podlegające rozkładowi jednostajnemu. Zastanowimy sie teraz jak mając zmienne losowe z rozkładu jednostajnego wytworzyć zmienne losowe o innych rozkładach.
Rozkład dwumianowy
Zmienna losowa, która zlicza liczbę sukcesów [math]k[/math] w [math]n[/math] próbach, gdzie [math]p[/math] jest prawdopodobieństwem sukcesu w pojedynczej próbie, podlega rozkładowi dwumianowemu [więcej o rozkładzie dwumianowym]:
[math]P_n(k)=\left( \begin{array}{c} n \\ k \end{array} \right) p^kq^{n-k}=\frac{n!}{k!(n-k)!}p^kq^{n-k} [/math]
Jak z rozkładu jednostajnego wytworzyć zmienne losowe o rozkładzie dwumianowym?
- Wykonać symulację!
# -*- coding: utf-8 -*-
import numpy as np
import pylab as py
def gen_binom(n,p,N):
'''funkcja zwracająca zmienną losową z rozkładu dwumianowego:
ilość sukcesów w n próbach przy prawdopodobieństwie sukcesu p
parametry:
n: liczba prób
p: prawdopodobieństwo sukcesu
N: liczba liczb do generowania
funkcja zwraca:
N liczb będącymi liczbami sukcesów
'''
k = np.zeros(N)
# w pętli wytwarzam pojedyncze zmienne losowe
# losuję n liczb z rozkładu jednostajnego [0,1) i jako sukcesy
#traktuję wylosowanie liczby nie większej niż p
r=np.random.rand(N,n)
k=np.sum(r<=p,axis=1)
return k
# kod testujący
n=20
N=100000
p=0.8
x = gen_binom(n,p,N)
bins =range(22) #granice binów
py.hist(x,bins)
py.show()
Zadanie
Z oszacowań agencji wynika, że średnio 2 z 3 reklam spotyka się z pozytywnym odzewem. Akcja marketingowa obejmuje 12 reklam. Niech [math]X[/math] oznacza liczbę reklam skutecznych. Czy [math]X[/math] podlega rozkładowi dwumianowemu? Jakie jest prawdopodobieństwo, że 10 reklam będzie skutecznych? Prawdopodobieństwo obliczyć ze wzoru oraz korzystając z symulacji.
Wskazówka informacje o zmiennych z rozkładu dwumianowego (binom) można uzyskać tak:
# >>> print(st.binom.extradoc) # Binomial distribution # # Counts the number of successes in *n* independent # trials when the probability of success each time is *pr*. # # binom.pmf(k,n,p) = choose(n,k)*p**k*(1-p)**(n-k) # for k in {0,1,...,n}
Odp: p = 0,127
Zadanie
Na egzaminie testowym jest 30 pytań. Na każde pytanie są podane cztery możliwe odpowiedzi, z czego tylko jedna jest prawdziwa. Za prawidłową odpowiedź student otrzymuje 1 punkt a za fałszywą −0,5 punktu. Próg zaliczenia wynosi 15 punktów. Oblicz prawdopodobieństwo, że udzielając czysto losowych odpowiedzi student zaliczy egzamin.
Rozkład Poissona
[math] P(k)=\frac{\mu^k e^{-\mu}}{k!} [/math] Rozkładowi Poissona podlegają zmienne losowe zliczające w jednostce czasu liczbę zdarzeń o niskim prawdopodobieństwie zajścia. Np. liczba rozpadów promieniotwórczych na jednostkę czasu [więcej o rozkładzie Poissona]
Przykład
Lekarz pełniący dyżur w szpitalu jest wzywany do pacjentów średnio 3 razy w ciągu nocy. Załóżmy, że liczba wezwań na noc podlega rozkładowi Poissona. Jakie jest prawdopodobieństwo, że noc upłynie lekarzowi spokojnie? U nas [math]\mu=3[/math], [math]x=0[/math] więc
[math] P(0) = \frac{3^0e^{-3}}{0!} = e^{-3} = 0,0498 [/math]
Przykład
Jak ze zmiennych podlegających rozkładowi jednostajnemu uzyskać zmienne podlegające rozkładowi Poissona?
Wykonać symulację wynikającą z następujących spostrzeżeń:
- Weźmy jednostkę czasu. Ma w niej zachodzić średnio [math]\mu[/math] zdarzeń.
- Podzielmy jednostkę czasu na NT odcinków o jednakowej długości.
- Prawdopodobieństwo, że zdarzenie zajdzie w konkretnym odcinku jest [math]p=\frac{\mu}{NT}[/math] Takie prawdopodobieństwo ma też wylosowanie przy pomocy generatora liczb z rozkładu jednostajnego [0,1) liczby nie większej niż p.
- Zatem, aby wytworzyć jedną liczbę z rozkładu Poissona można policzyć ile spośród NT liczb z rozkładu jednostajnego [0,1) jest nie większych niż p.
# -*- coding: utf-8 -*-
import numpy as np
import pylab as py
import scipy.stats as st
def gen_poiss(mi,NT,N):
'''funkcja zwracająca zmienne losowe z rozkładu poissona:
średnia ilość zdarzeń o niskim prawdopodobieństwie w jednostce czasu
parametry:
mi: średnia liczba zdarzeń
NT: liczba odcinków czasu
N: liczba liczb do generowania
funkcja zwraca:
N liczb będących zliczeniami zdarzeń
'''
p = float(mi)/NT
k = np.zeros(N)
# w pętli wytwarzam pojedyncze zmienne losowe
# losuję NT liczb z rozkładu jednostajnego [0,1) i jako zajście zdarzenia
# traktuję wylosowanie liczby nie większej niż p. Ilość zdarzeń to liczba z rozkładu Poissona
r=np.random.rand(N,NT)
k=np.sum(r<=p,axis=1)
return k
# kod testujący
mi=4
N=100000
NT = 10000
x = gen_poiss(mi,NT,N)
bins =np.arange(int(x.max()+1))-0.5
py.hist(x,bins,normed=True)
# porównanie z generatorami scipy.stats
p=mi/NT
k=np.arange(x.max()+1)
B=st.binom.pmf(k,NT,p)
py.plot(k,B,'ro')
P=st.poisson.pmf(k,mi)
py.plot(k,P,'go')
py.show()
Rozkład normalny i Centralne Twierdzenie Graniczne
Gęstość prawdopodobieństwa rozkładu normalnego dana jest wzorem:
[math]f(x)=\frac{1}{\sqrt{2 \pi} \sigma} e^{ -\frac{(x-\mu)^2}{2 \sigma^2}}[/math] gdzie: [math]\mu[/math] — średnia, [math]\sigma[/math] — odchylenie standardowe ([math]\sigma^2[/math] — wariancja).
Jeden ze sposóbw wytwarzania zmiennych losowych o rozkładzie normalnym polega na zastosowaniu Centralnego Twierdzenia Granicznego.
[math]X=\frac{\sqrt{12}}{\sqrt{N}}\left(\sum_{n=1}^N Y_n -\frac{N}{2}\right)[/math]
gdzie [math]Y[/math] zmienna losowa z rozkładu jednostajnego [math][0,1)[/math] W granicy dużych [math]N[/math] rozkład zmiennej [math]X[/math] dąży do rozkładu normalnego [math]N(0,1)[/math].
Zadanie
Proszę:
- napisać funkcję gen_CTG_1(N), która zwraca liczbę losową X otrzymaną zgodnie powyższym wzorem przez sumowanie N liczb z rozkładu jednostajnego.
- napisać funkcję gen_CTG(N, N_liczb) korzystającą z gen_CTG_1(N) i zwracającą zadaną ilość liczb N_liczb
- narysować histogramy 10 000 liczb [math]X[/math] uzyskanych kolejno dla [math]N={1,2,...,12}[/math]
- Jakie parametry charakteryzują rozkład do którego zbiegają sumy? Przy każdym z powyższych histogramów wypisz średnią i wariancję z próby (funkcje np.mean i np.var).
Zadanie: transformacja rozkładu normalnego
Rozkład o średniej 0 i wariancji 1 (notacja [math]N(0,1)[/math]) jest nazywany rozkładem standardowym i często jest oznaczany literą [math]Z[/math]. Dokonując odpowiedniej transformacji można z rozkładu [math]Z[/math] uzyskać dowolny inny rozkład normalny.
Proszę:
- wygenerować 100000 zmiennych losowych z rozkładu [math]N(0,1)[/math] (skorzystać z funkcji st.norm.rvs(loc=0, scale=1, size=100000))
- przetransformować je do zmiennych losowych [math]N(2,9)[/math].
- narysować histogramy zmiennych przed i po transformacji.
- narysować rozkład gęstości prawdopodobieństwa dla rozkładu [math]N(2,9)[/math] (skorzystać z funkcji st.norm.pdf(...)) .
Transformacja Boxa-Mullera
To inna popularna metoda metoda generowania liczb losowych o rozkładzie normalnym, na podstawie dwóch wartości zmiennej o rozkładzie jednostajnym na przedziale [0,1) http://en.wikipedia.org/wiki/Box-Muller_transform, ale nie będziemy się nią tu szerzej zajmować.
Przykład
Producent silników twierdzi, że jego silniki mają średnią moc 220 KM, a odchylenie standardowe wynosi 15 KM. Potencjalny klient testuje 100 silników. Jakie jest prawdopodobieństwo, że średnia z próby będzie mniejsza niż 217 KM?
Przypomnijmy, że z CTG dla dużych liczebności próby [math]n[/math] [math]\bar x \sim N(\mu,\sigma^2/n)[/math] (dowód). Zatem szukamy
[math]P(\bar x \lt 217) =P\left( Z\lt \frac{217-\mu}{\frac{\sigma}{ \sqrt{n}}} \right) = P\left( Z\lt \frac{217-220}{\frac{15}{ \sqrt{100}}} \right)=P(Z\lt -2)=0,0228[/math]
import scipy.stats as st
import pylab as py
import numpy as np
# ze wzoru
p=st.norm.cdf(-2)
print('prawdopodobieństwo odczytane z dystrybuanty rozkładu normalnego: %(0).4f' %{'0':p})
# symulacja
mu=220
sig=15
N_prob=100
m_kryt=217
N_rep=int(1e4)
srednia=np.zeros(N_rep)
for i in range(N_rep):
seria=st.norm.rvs(loc=mu, scale=sig, size=N_prob)
srednia[i]=np.mean(seria)
h=py.hist(srednia,30);
py.plot([m_kryt, m_kryt],[0, max(h[0])],'r')
p1=np.sum(srednia<=m_kryt)/N_rep
print('prawdopodobieństwo uzyskane z symulacji: %(0).4f' %{'0':p1})
py.show()
Zadania
Teraz kilka najprostszych zadań dotyczących rozkładu normalnego:
- Znajdźmy prawdopodobieństwo, że Z < −2,47. Proszę zrobić to na dwa sposoby: raz z użyciem wygenerowanych zmiennych z rozkładu normalnego, drugi raz z użyciem dystrybuanty norm.cdf
[Odp: p = 0,0068]
Rozkład t
Gdy nie znamy wariancji populacji [math]\sigma^2[/math] używamy w jego miejsce estymatora wariancji próby [math]s^2[/math] danego wzorem:
[math]s^2 =\frac{\sum(\bar x - x_i)^2}{n-1} [/math] gdzie n — rozmiar próby
Wtedy rozkład [math]Y=\frac{\bar X -\mu}{\frac{s}{\sqrt{n}}}[/math] nie podlega rozkładowi normalnemu. Jeśli populacja [math]X[/math] podlega rozkładowi normalnemu to [math]Y[/math] podlega rozkładowi t z n−1 stopniami swobody.
Zobaczmy na symulacji co to zmienia. Proszę obejrzeć wyniki dla kilku wartości N
# -*- coding: utf-8 -*-
import numpy as np
import pylab as py
import scipy.stats as st
N=3
N_liczb=10000
mu=2
sigma=1
y_wyidealizowane=np.zeros(N_liczb)
m_populacji=np.zeros(N_liczb)
s_populacji=np.zeros(N_liczb)
y_estymowane=np.zeros(N_liczb)
m_estymowane=np.zeros(N_liczb)
s_estymowane=np.zeros(N_liczb)
for i in range(N_liczb):
x=st.norm.rvs(loc=mu,scale=sigma,size=N)
m_populacji[i]=mu
s_populacji[i]=sigma/np.sqrt(N)
m_estymowane[i]=np.mean(x)
s_estymowane[i]=np.std(x,ddof=1)/np.sqrt(N)
y_wyidealizowane[i]=(m_estymowane[i]-mu)/(s_populacji[i])
y_estymowane[i]=(m_estymowane[i]-mu)/(s_estymowane[i])
py.subplot(311)
py.title(u'średnia (wartość oczekiwana)')
py.plot(m_estymowane,label="estymowana")
py.plot(m_populacji,'r',label="populacji")
py.legend()
py.subplot(312)
py.title('odchylenie standardowe')
py.plot(s_estymowane,label="estymowane")
py.plot(s_populacji,'r',label='populacji')
py.legend()
py.subplot(325)
py.title(u'rozkład ilorazu Y dla wersji wyidealizowanej')
py.hist(y_wyidealizowane,30,(-4,4))
py.subplot(326)
py.title(u'rozkład ilorazu Y dla wersji estymowanej')
py.hist(y_estymowane,30,(-4,4))
py.show()
Zadanie
Dział kontroli jakości producenta silników testuje 10 egzemplarzy silnika nowego typu. Uzyskano wartość średnią 220 KM oraz odchylenie standardowe równe 15 KM. Jakie jest prawdopodobieństwo, że klient, który zamówi 100 silników, otrzyma partię, w której średnia będzie mniejsza niż 217 KM? Wykonaj obliczenia oraz symulację. Porównaj wynik z wcześniejszym przykładem, gdzie znaliśmy odchylenie standardowe z innego źródła (nie było ono estymowane z próby).
Inne rozkłady uzyskiwane przez transformacje
Metoda odwracania dystrybuanty
Niech [math]\xi[/math] będzie zmienną losową o rozkładzie jednostajnym na przedziale [math][0,1)[/math]. Jej dystrybuanta [math]J[/math] dana jest więc worem:
[math] J(x) = P(\xi \le x) = \left\{ \begin{array}{lrl} 0 & \mathrm{dla} & x\lt 0\\ x & \mathrm{dla} & 0 \le x \lt 1 \\ 1 & \mathrm{dla} & 1 \le x \end{array} \right. [/math]
oraz niech [math]F[/math] będzie dystrybuantą interesującego nas rozkładu. Jeśli [math]F[/math] jest funkcją odwracalną to zmienna losowa zdefiniowana jako:
[math]\eta = F^{-1}(\xi) [/math]
jest zmienną losową podlegającą rozkładowi o dystrybuancie [math] F[/math]. Własność tę łatwo sprawdzić. Ponieważ dla każdego [math]x \in R [/math] mamy:
[math] P(\eta \le x) = P(F^{-1} (\xi) \le x) = P(\xi \le F(x)) = F(x) [/math] więc [math]F[/math] jest dystrybuantą zmiennej losowej [math]\eta[/math].
Czyli jeśli znamy dystrybuantę interesującego nas rozkładu prawdopodobieństwa i dystrybuanta ta jest odwracalna, to możemy generować liczby z rozkładu o dystrybuancie [math]F[/math].
Przykład
Wygenerować liczby pseudolosowe z rozkładu wykładniczego o gęstości prawdopodobieństwa:
[math]f(x) = \lambda e^{-\lambda x} [/math]
Rozwiązanie:
- Znajdujemy dystrybuantę tego rozkładu:
- [math]F(x) = 1-e^{-\lambda x} [/math]
- Znajdujemy funkcję odwrotną do dystrybuanty:
- [math]G(u)=F^{-1}(u) = -\frac{1}{\lambda}\ln(1-u) [/math]
- Generujemy liczby losowe [math]u[/math] z przedziału [math][0,1)[/math] i stosujemy transformację [math]G(u)[/math]
- Test generowanego rozkładu robimy przez porównanie dystrybuanty teoretycznej i empirycznej. Proszę zwrócić uwagę na funkcję estymującą dystrybuantę empiryczną
# -*- coding: utf-8 -*-
import pylab as py
import numpy as np
def gen_wykladniczy(lam, N):
'''funkcja zwraca N liczb z rozkładu wykładniczego o parametrze lam
korzystając z metody odwracania dystrybuanty
Dystrybuanta rozkładu wykładniczego: F(x) = 1- exp(-lam*x)
Funkcja do niej odwrotna: G(u) = -1/lam ln(1-u)'''
u = np.random.random(N)
return -1.0/lam* np.log(1.-u)
def dystrybuanta(X,od=0, do=4, krok=0.05):
'''funkcja zwraca dystrybuantę empiryczną
zmiennej X na przedziale (od, do) z krokiem krok'''
os_x = np.arange(od,do,krok)
nx = len(os_x)
N_liczb = len(X)
D=np.zeros(nx) #wyzerowany wektor na naszą dystrybuantę
for i in range(nx):
D[i] = sum(X <= os_x[i])/N_liczb
return (os_x, D)
# testujemy: narysujemy dystrybuantę oczekiwaną i empiryczną
lam=3
N=10000
x = np.arange(0,4,0.01)
F = 1-np.exp(-lam*x)
X = gen_wykladniczy(lam, N)
(os_x, F_emp) = dystrybuanta(X,od=0, do=4, krok=0.1)
py.plot(x,F,'r',label = 'dystrybuanta teoretyczna')
py.plot(os_x, F_emp,'.b', label = 'dystrybuanta empiryczna')
py.legend()
py.show()
Zauważmy, że nie zawsze jest łatwo wyznaczyć efektywnie [math]F^{-1}[/math] — jest tak na przykład w przypadku rozkładu normalnego. W takich sytuacjach można szukać przybliżonej wartości [math]F^{-1}[/math], rozwiązując numerycznie równanie: [math]F(x)=u[/math], co sprowadza się do znalezienia miejsc zerowych funkcji [math]H(x)=F(x)-u[/math].
Zadanie
Posiłkując się powyższym przykładem napisz funkcję MyGen(N) dla rozkładu, którego funkcja gęstości prawdopodobieństwa jest dana wzorem: [math] f(x) = \frac{1}{\pi (1+x^{2})} [/math] Następnie narysuj dla niego dystrybuantę i dystrybuantę empiryczną na przedziale [-10, 10].