4.6. Statistika

Sage ima implementirane osnovne statističke funkcije, ali nama će biti zgodnije koristiti paket scipy.stats čije mogućnosti su trenutno veće.

sage: import scipy
sage: import scipy.stats
sage: import matplotlib.pyplot as plt
sage: import numpy as np

Kreirajmo prvo NumPy listu brojeva. NumPy liste imaju metode za izračun osnovnih statističkih veličina poput srednje vrijednosti, varijance, standardne devijacije itd.

sage: xs=np.linspace(1,9,21); xs
array([ 1. ,  1.4,  1.8,  2.2,  2.6,  3. ,  3.4,  3.8,  4.2,  4.6,  5. ,
        5.4,  5.8,  6.2,  6.6,  7. ,  7.4,  7.8,  8.2,  8.6,  9. ])
sage: print "%.2f +- %.2f" % (xs.mean(), xs.std())
5.00 +- 2.42

Obične liste nemaju gornje statističke metode pa ukoliko nas zanima npr. srednja vrijednost brojeva u običnoj listi treba je prvo konvertirati u numpy listu funkcijom np.array() ili možemo primijeniti Sage funkciju mean() (dakle ne metodu).

sage: mean([3, 5, 7, 9, 11])
7

scipy.stats paket ima implementirane brojne statističke raspodjele:

  • norm - normalna (Gaussova) raspodjela
  • poisson - Poissonova raspodjela
  • gamma - gamma raspodjela
  • chi2 - \(\chi^2\) raspodjela
  • t - studentova t raspodjela
  • ...

Svaka od tih raspodjela ima svoje parametre (npr. za Gaussovu raspodjelu su to srednja vrijednost i standardna devijacija). Oni se najčešće specificiraju putem opcionalnih argumenata, a točnu sintaksu saznajemo iz standardne dokumentacije.

Najvažnije metode distribucija su:

  • pdf() - sama distribucija (probability distribution function) \(f(x)\)
  • cdf() - integral distribucije (cumulative distribution function) \(\Phi(x)=\int_{-\infty}^{x} f(y) dy\)
  • rvs() - slučajni brojevi (random variates) koji slijede distribuciju

Tako je npr. vrijednost u \(x=0\) normalne distribucije sa srednjom vrijednosti \(\mu = 5\) i standardnom devijacijom \(\sigma = 1.5\), \(f(x=0, \mu=5, \sigma=1.5)\) dan s:

sage: scipy.stats.norm.pdf(0, loc=5., scale=1.5)
0.0010281859975274036

Skica distribucije i njenog integrala:

sage: xs = np.linspace(-2,12)
sage: fig, ax = plt.subplots(figsize=[5,3])
sage: ax.plot(xs, scipy.stats.norm.pdf(xs, loc=5., scale=1.5), label='pdf')
sage: ax.plot(xs, scipy.stats.norm.cdf(xs, loc=5., scale=1.5), 'r--',
....:                                                          label='cdf')
sage: ax.legend(loc='upper left').draw_frame(0)
sage: fig.savefig('fig')
../_images/gauss.png

Izračunajmo vjerojatnost da se slučajna varijabla nađe unutar jedne standardne devijacije \(\sigma\) od srednje vrijednosti \(\mu\). Ona je dana integralom distribucije od \(\mu-\sigma\) do \(\mu+\sigma\). Kako na raspolaganju imamo funkciju cdf() koja daje integral od \(-\infty\) do \(x\), treba samo oduzeti dva takva integrala (koristimo i to da su defaultne vrijednosti loc=0 i scale=1):

sage: scipy.stats.norm.cdf(1) - scipy.stats.norm.cdf(-1)
0.68268949213708585

(To je čuvenih 68% vjerojatnosti.)

Zadatak 1

Slučajna varijabla \(X\) slijedi normalnu raspodjelu sa srednjom vrijednosti \(\mu=50\) i standardnom devijacijom \(\sigma=2\). Kolika je vjerojatnost da se \(X\) nađe između 47 i 54?

Sad ćemo metodom rvs() izgenerirati uzorak od 10000 brojeva raspodjeljenih po normalnoj raspodjeli s \(\mu=5\) i \(\sigma=1.5\) i testirati da su srednja vrijednost i standardna devijacija prema očekivanjima:

sage: pts=scipy.stats.norm.rvs(loc=5., scale=1.5, size=1e4)
sage: print "broj točaka = %i" % len(pts)
broj točaka = 10000
sage: print "srednja vrijednost = %.3f" % pts.mean()  # rel tol 0.01
srednja vrijednost = 4.998
sage: print "standardna devijacija = %.3f" % pts.std()  # rel tol 0.03
standardna devijacija = 1.493

Za crtanje histograma ovog uzorka koristimo Matplotlibovu hist() funkciju, čiji opcionalni argument bins kontrolira broj “binova” dakle broj intervala apscise u kojima se prebrojavaju točke:

sage: fig, ax = plt.subplots(figsize=[3.5,2])
sage: ax.hist(pts, bins=20, normed=True)
sage: fig.savefig('fig')
../_images/histogram.png

Ovo je bilo za jedan uzorak. Listu npr. standardnih devijacija za 5 nezavisnih uzoraka možemo konstruirati ovako:

sage: [scipy.stats.norm.rvs(loc=5, scale=1.5, size=1e4).std() for k in
....:                                            range(5)] # rel tol 0.05
[1.503340173922008, 1.5040004762159596, 1.4785845284471362,
1.4906248759178866, 1.5136584591557063]

Zadatak 2

Kolika je vjerojatnost da mjerenje slučajne varijable koja slijedi normalnu razdiobu sa srednjom vrijednošću \(\mu\) i standardnom devijacijom \(\sigma\) odstupi od srednje vrijednosti za više od \(3\sigma\)?

Važna raspodjela je tzv. \(\chi^2\)-raspodjela

\[P(\chi^2, \nu) = \frac{(\chi^2)^{\nu/2-1}e^{-\chi^2/2}}{2^{\nu/2} \Gamma(\nu/2)}\]

Vidimo da \(\chi^2\) raspodjela, pored argumenta x, ima samo jedan parametar: \(\nu\) = df = broj stupnjeva slobode (degrees of freedom) [1].

Srednja vrijednost \(\chi^2\) raspodjele jednaka je broju stupnjeva slobode:

sage: scipy.stats.chi2.rvs(3, size=1e5).mean()  # rel tol 0.01
3.0033578764955071

Jedno od svojstava \(\chi^2\) raspodjele je da integral repa te raspodjele (s jednim stupnjem) od 1, 4, 9 do \(\infty\) daje vjerojatnost da mjerenje slučajne varijable raspodjeljene po Gaussovoj raspodjeli sa standardnom devijacijom \(\sigma\) odstupi od srednje vrijednosti za \(\sigma, 2\sigma, 3\sigma, \ldots\).

sage: [1-scipy.stats.chi2.cdf(eps^2,1) for eps in [1,2,3]]
[0.31731050786291404, 0.045500263896358528, 0.0026997960632602069]

Zadatak 3

Središnji granični teorem kaže da su srednje vrijednosti dovoljno velikih uzoraka neke raspodjele raspodjeljene prema normalnoj (Gaussovoj) raspodjeli bez obzira na to kakva je originalna raspodjela iz koje te uzorke uzimamo. U to ćemo se uvjeriti na primjeru \(\chi^2\) raspodjele.

  1. Nacrtajte tu raspodjelu za df=3 stupnja slobode i uvjerite se da je asimetrična oko srednje vrijednosti.
  2. Uzmite jedan uzorak od milijun točaka te raspodjele i uvjerite se da je standardna devijacija \(\sigma=\sqrt{2\cdot df}\)
  3. Uzmite 500 uzoraka od po 400 točaka svaki i napravite listu srednjih vrijednosti tih uzoraka. Uvjerite se da je standardna devijacija vrijednosti u listi \(\sigma/\sqrt{400}\), kako traži teorem.
  4. Nacrtajte histogram srednjih vrijednosti u listi i uvjerite se (superponiranjem krivulje normalne razdiobe) da su one zaista raspodjeljene prema normalnoj razdiobi srednje vrijednosti \(\mu=df=3\) i standardne devijacije \(\sigma/\sqrt{400}\)

Za ove iste primjere, ali izvedene pomoću Sage funkcija, vidi odjeljak Statistika - alternative.

Bilješke

[1]Naziv tog parametra postaje jasan u kontekstima određivanja dobrote prilagodbe i statističkog testiranja hipoteze.

Prijašnja tema

4.5. Diferencijalne jednadžbe

Slijedeća tema

4.7. Prilagodba funkcije podacima