WnioskowanieStatystyczne/Zmienne losowe i generatory liczb pseudolosowych: Różnice pomiędzy wersjami

Z Brain-wiki
 
(Nie pokazano 92 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<x}{P(X = x_i) } =  \sum_{x_i<x}p_i</math>
+
: <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 &mdash; 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 &mdash; 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]) \geqslant p
+
P_X((-\infty, x_p]) \leqslant p
 
</math>
 
</math>
 
   
 
   
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 #8191
+
     m=8765132137
     a=42347#101
+
     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):
Linia 241: Linia 248:
  
  
N=50000
+
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 metody ===
+
=== 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));
  
Dla przykładu zobaczmy działanie tych metod dla zmiennych losowych z rozkładu normalnego:
+
====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 dystrybucji ===
+
=== Pomoc do poszczególnych rozkładów ===
Pełną dokumnetację każdej dystrybucji mamy dostępną w postaci docstringów. Dla przykładu:
+
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 docstring o niecentralnym rozkładzie ''t''.
+
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: ilość prób
+
n: liczba prób
 
p: prawdopodobieństwo sukcesu
 
p: prawdopodobieństwo sukcesu
N: ilość liczb do generowania
+
N: liczba liczb do generowania
 
funkcja zwraca:
 
funkcja zwraca:
N liczb będącymi ilościami sukcesów
+
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
for i in range(N):
+
r=np.random.rand(N,n)
r = np.random.random(size=n)
+
k=np.sum(r<=p,axis=1)
k[i]=np.sum(r<=p)
 
 
return k
 
return k
 +
  
 
# kod testujący
 
# kod testujący
Linia 416: Linia 586:
 
  #  for k in {0,1,...,n}
 
  #  for k in {0,1,...,n}
 
Odp: ''p'' = 0,127
 
Odp: ''p'' = 0,127
 +
  
 
<!--
 
<!--
 
=====Rozwiązanie:=====
 
=====Rozwiązanie:=====
 
<source lang = python>  
 
<source lang = python>  
# -*- coding: utf-8 -*-
 
 
import numpy as np
 
import numpy as np
 
import pylab as py
 
import pylab as py
Linia 426: Linia 596:
 
import numpy.random as rnd
 
import numpy.random as rnd
  
p=2./3 #prawdopodobieństwo sukcesu w jednej próbie
+
p=2/3 #prawdopodobieństwo sukcesu w jednej próbie
n=12  # ilość prób w jednej akcji reklamowej
+
n=12  # liczba prób w jednej akcji reklamowej
  
 
# ze wzoru  
 
# ze wzoru  
Linia 442: Linia 612:
 
print(p_wz)
 
print(p_wz)
 
# z symulacji
 
# z symulacji
N_powt =100000 # ilość powtórzeń symulacji
+
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 454: Linia 624:
 
====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ą &minus;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ą &minus;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 ilość punktów wymaganych do zaliczenia zgodnie z równaniem:
+
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:
n-0.5*(30-n) = 15
+
k-0.5*(30-k) = 15
co daje n = 20
+
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 468: Linia 637:
  
 
p = 1./4 #prawdopodobieństwo sukcesu w jednej próbie
 
p = 1./4 #prawdopodobieństwo sukcesu w jednej próbie
n = 30  # ilość prób w jednej akcji reklamowej
+
n = 30  # liczba prób w jednej akcji reklamowej
 
k_kryt = 20
 
k_kryt = 20
 
# ze wzoru  
 
# ze wzoru  
Linia 477: 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) # ilość sukcesów
+
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>
 
-->
 
-->
  
Linia 489: 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 ilość zdarzeń o niskim prawdopodobieństwie zajścia. Np. ilość 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]]
+
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 513: 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 ilość zdarzeń
+
mi: średnia liczba zdarzeń
NT: ilość odcinków czasu
+
NT: liczba odcinków czasu
N: ilość liczb do generowania
+
N: liczba liczb do generowania
 
funkcja zwraca:
 
funkcja zwraca:
 
N liczb będących zliczeniami zdarzeń
 
N liczb będących zliczeniami zdarzeń
 
'''
 
'''
p = mi/NT
+
p = float(mi)/NT
 
k = np.zeros(N)
 
k = np.zeros(N)
 
# w pętli wytwarzam pojedyncze zmienne losowe
 
# w pętli wytwarzam pojedyncze zmienne losowe
 
# 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
for i in range(N):
+
r=np.random.rand(N,NT)
r = np.random.random(size=NT)
+
k=np.sum(r<=p,axis=1)
k[i]=np.sum(r<=p)
 
 
return k
 
return k
+
 
# kod testujący
 
# kod testujący
 
mi=4
 
mi=4
N=10000
+
N=100000
 
NT = 10000
 
NT = 10000
 
x = gen_poiss(mi,NT,N)
 
x = gen_poiss(mi,NT,N)
bins =range(31)
+
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 547: 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}} e^{ -\frac{(x-\mu)^2}{2 \sigma^2}}</math>
+
<math>f(x)=\frac{1}{\sqrt{2 \pi} \sigma} e^{ -\frac{(x-\mu)^2}{2 \sigma^2}}</math>
 
gdzie: <math>\mu</math> &mdash; średnia, <math>\sigma</math> &mdash; odchylenie standardowe (<math>\sigma^2</math> &mdash; wariancja).
 
gdzie: <math>\mu</math> &mdash; średnia, <math>\sigma</math> &mdash; odchylenie standardowe (<math>\sigma^2</math> &mdash; wariancja).
  
Linia 613: Linia 791:
 
-->
 
-->
  
====Problem: transformacja rozkładu normalnego ====
+
====Zadanie: transformacja rozkładu normalnego ====
  
  
Linia 624: Linia 802:
 
* 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 =====
 
===== Rozwiązanie =====
  
Linia 657: Linia 835:
 
py.show()
 
py.show()
 
</source>
 
</source>
 +
-->
  
 
====Transformacja Boxa-Mullera====
 
====Transformacja Boxa-Mullera====
Linia 703: Linia 882:
 
Teraz kilka najprostszych zadań dotyczących rozkładu normalnego:
 
Teraz kilka najprostszych zadań dotyczących rozkładu normalnego:
 
* Znajdźmy prawdopodobieństwo, że ''Z'' < &minus;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'' < &minus;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]
+
[Odp: ''p'' = 0,0068]
 
+
<!--
 
Rozwiązanie:
 
Rozwiązanie:
 
<source lang =python>
 
<source lang =python>
Linia 718: Linia 897:
 
</source>
 
</source>
  
* Znaleźć prawdopodobieństwo ''P''(|''Z''| < 2) [Odp:  ''p'' = 0,9545]
+
* Znaleźć prawdopodobieństwo ''P''(|''Z''| < 2)
 +
[Odp:  ''p'' = 0,9545]
 +
 
  
 
Rozwiązanie:
 
Rozwiązanie:
Linia 735: Linia 916:
  
 
* 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. [Opd: ''p'' = 0,852]
+
Prawdopodobieństwo obliczyć korzystając z dystrybuanty rozkładu normalnego oraz z symulacji.
 +
[Odp: ''p'' = 0,852]
  
 
Rozwiązanie:
 
Rozwiązanie:
Linia 759: Linia 941:
 
py.show()
 
py.show()
 
</source>
 
</source>
 +
-->
  
 
=== Rozkład t ===
 
=== Rozkład t ===
Linia 821: Linia 1004:
 
py.show()
 
py.show()
 
</source>
 
</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()
 +
 +
</source>
 +
-->
  
 
=== Inne rozkłady uzyskiwane przez transformacje ===
 
=== Inne rozkłady uzyskiwane przez transformacje ===
Linia 855: Linia 1081:
  
 
====Przykład====  
 
====Przykład====  
Wygenerować liczby pseudolosowe z rozkładu wykładniczego o dystrybuancie:
+
Wygenerować liczby pseudolosowe z rozkładu wykładniczego o gęstości prawdopodobieństwa:
  
<math>F(X) = 1-e^{-\lambda x} </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:
  
Linia 912: Linia 1142:
 
Zauważmy, że nie zawsze jest łatwo wyznaczyć efektywnie  <math>F^{-1}</math> &mdash; 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> &mdash; 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

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:

Rozklad moneta 1.gif

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]

Funkcja gestosci prawdopodobienstwa.png

  • 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ą:

Rozklad moneta.svg

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:

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:

  1. Wygeneruj 10000 liczb z rozkładu t i narysuj ich histogram.
  2. Wykreśl gęstość prawdopodobieństwa tego rozkładu.
  3. Wykreśl dystrybuantę.
  4. 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:

  1. Rozkład normalny [math]\mu[/math] = 100, [math]\sigma[/math] = 2.
  2. Rozkład dwumianowy p = 0.4, n =10.
  3. Rozkład poissona [math]\lambda[/math] = 4.
  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:

  1. Wartość oczekiwana.
  2. Odchylenie standardowe.
  3. Wariancja.
  4. Mediana.
  5. 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].