Sage ima implementirane osnovne statističke funkcije (vidi niže odjeljak "čisti Sage"), ali nama će biti zgodnije koristiti paket scipy čije mogućnosti su trenutno veće.
{{{id=269| import scipy import scipy.stats import pylab /// }}}Kreirajmo prvo (numpy) listu brojeva. Numpy liste imaju metode za izračun osnovnih statističkih veličina poput srednje vrijednosti, varijance, standardne devijacije itd.
{{{id=386| xs=pylab.linspace(1,9,21); xs /// }}} {{{id=387| print "sr. vrijednost +- std. devijacija = %.2f +- %.2f" % (xs.mean(), xs.std()) print "varijanca = %.3f" % xs.var() /// }}}(Ovdje smo za formatiranje ispisa primjenili operator %, te standardne stringove za formatiranje (poznate iz sprintf() C funkcije): Npr. %5.3f formatira realni broj na 5 mjesta ukupno i 3 mjesta iza decimalne točke.)
Obične liste nemaju gornje statističke metode ...
{{{id=388| t1 = [3, 5, 7, 9, 11] t1.mean() /// }}}... pa ukoliko nas zanima npr. srednja vrijednost brojeva u običnoj listi treba je prvo konvertirati u numpy listu funkcijom pylab.array():
{{{id=389| t1num = pylab.array(t1) t1num.mean() /// }}}... ili možemo primijeniti istoimenu Sage funkciju (dakle ne metodu):
{{{id=427| mean(t1) /// }}}scipy.stats paket ima implementirane brojne statističke raspodjele:
Svaka od tih raspodjela ima svoje parametre (npr. srednja vrijednost i standardna devijacija), koji se najčešće specificiraju putem opcionalnih argumenata loc i scale, čije točno značenje saznamo iz standardne dokumentacije ili isprintavanjem specijalnog konciznog atributa 'extradoc':
{{{id=407| print scipy.stats.norm.extradoc /// }}}Najvažnije metode distribucija su:
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:
{{{id=390| scipy.stats.norm.pdf(0, loc=5., scale=1.5) /// }}}Skica distribucije i njenog integrala:
{{{id=273| pylab.figure(figsize=[5,3]) xs = pylab.linspace(-2,12) pylab.plot(xs, scipy.stats.norm.pdf(xs, loc=5., scale=1.5), label='pdf') pylab.plot(xs, scipy.stats.norm.cdf(xs, loc=5., scale=1.5), 'r--', label='cdf') pylab.legend(loc='upper left').draw_frame(0) pylab.savefig('tmp') /// }}}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:
{{{id=413| scipy.stats.norm.cdf(1) - scipy.stats.norm.cdf(-1) # defaulti su mu=0, sigma=1 /// }}}(To je čuvenih 68% vjerojatnosti.)
♦ Zadatak 2-7.1: ("Primjer 2" is SOM skripti M. Požeka) 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 radiobi i testirati da su srednja vrijednost i standardna devijacija prema očekivanjima:
{{{id=274| pts=scipy.stats.norm.rvs(loc=5., scale=1.5, size=1e4) print "broj točaka = %i" % len(pts) print "srednja vrijednost = %.3f" % pts.mean() print "standardna devijacija = %.3f" % pts.std() /// }}}Za crtanje histograma ovog uzorka koristimo pylab.hist() funkciju, čiji opcionalni argument bins kontrolira broj "binova" dakle broj intervala apscise u kojima se prebrojavaju točke:
{{{id=369| pylab.figure(figsize=[3.5, 2]) pylab.hist(pts, bins=20, normed=True) pylab.savefig('tmp') /// }}}Ovo je bilo za jedan uzorak. Listu npr. standardnih devijacija za 5 nezavisnih uzoraka možemo konstruirati obuhvaćanjem liste uz korištenje funkcije range() koja kreira listu cijelih brojeva:
{{{id=370| print range(5) [k^2 for k in range(5)] /// }}} {{{id=367| [scipy.stats.norm.rvs(loc=5, scale=1.5, size=1e4).std() for k in range(5)] /// }}}♦ Zadatak 2-7.2: Koja 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
{{{id=418| print scipy.stats.chi2.extradoc /// }}}Usporedi sa SOM (Požek p. 80):
$$ 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). Naziv tog parametra će biti jasniji niže kod analize prilagodbi funkcija točkama.
Srednja vrijednost $\chi^2$ raspodjele jednaka je broju stupnjeva slobode:
{{{id=420| scipy.stats.chi2.rvs(3, size=1e5).mean() /// }}}Vjerojatnost da mjerenje slučajne varijable raspodjeljene po normalnoj raspodjeli sa standardnom devijacijom $\sigma$ odstupi od srednje vrijednosti za $\sigma, 2\sigma, 3\sigma, \ldots$ je dano integralom repa $\chi^2$-raspodjele s jednim stupnjem slobode od 1, 4, 9 do $\infty$:
{{{id=159| [1-scipy.stats.chi2.cdf(eps^2,1) for eps in [1,2,3]] /// }}}♦ Zadatak 2-7.3: Središnji granični teorem (cf. Požek p. 94) između ostalog 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.
U ovom odjeljku ćemo, bez puno komentara, napraviti identične primjere kao gore, ali koristeći isključivo "native" Sage funkcije, a ne scipy i pylab.
{{{id=423| xs =[ 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. ] /// }}} {{{id=424| print "srednja vrijednost = %.1f" % mean(xs) print "varijanca = %.3f" % variance(xs, bias=True) # default bias=False dijeli s 1/(len(xs)-1) print "standardna devijacija = %.3f" % std(xs, bias=True) /// }}} {{{id=435| G = RealDistribution('gaussian', 1.5) G.distribution_function(-5) /// }}} {{{id=372| G.plot((-7, 7), figsize=[3.5,2]) /// }}} {{{id=437| G.cum_distribution_function(1.5)-G.cum_distribution_function(-1.5) /// }}} {{{id=436| pts = [5 + G.get_random_element() for k in range(1e4)] print "broj točaka = %i" % len(pts) print "srednja vrijednost = %.1f" % mean(pts) print "standardna devijacija = %.3f" % std(pts) # znatno sporije nego gore /// }}} {{{id=440| G.generate_histogram_plot('fig', num_samples=1e4, bins=20) /// }}} {{{id=443| [mean([5 + G.get_random_element() for k in range(10000)]) for k in range(5)] /// }}} {{{id=375| CHISQ = RealDistribution('chisquared', 3) # argument je broj stupnjeva slobode mean([CHISQ.get_random_element() for k in range(1000)]) /// }}} {{{id=377| CHISQ = RealDistribution('chisquared', 1) [1-CHISQ.cum_distribution_function(eps^2) for eps in [1,2,3]] /// }}}Često je potrebno neke podatke dobivene npr. mjerenjima grafički prikazati i približno opisati nekom funkcijom (tzv. prilagodba ili "fitanje").
{{{id=8| data = [[0.2, 0.1], [0.35, 0.2], [1, 0.6], [1.8, 0.9], [2.5, 1.3], [4, 2.06], [5, 2.6]] list_plot(data, figsize=[3,3]) /// }}}Klasična situacija je potreba za prilagodbom pravca $f(x) = a + b x$ metodom najmanjih kvadrata (linearna regresija). Pogledajmo prvo kako bismo sami implementirali formule za metodu najmanjih kvadrata:
{{{id=303| print "Prilagodba funkcije f(x) = a + b*x" xs = [x for x,y in data] ys = [y for x,y in data] xs2 = [x^2 for x,y in data] ys2 = [y^2 for x,y in data] xys = [x*y for x,y in data] n = len(data) delc = n*sum(xs2) - sum(xs)^2 b = (n*sum(xys)-sum(xs)*sum(ys))/delc a = (sum(xs2)*sum(ys)-sum(xs)*sum(xys))/delc sigsq = sum([(y-b*x-a)^2 for x,y in data])/(n-2) erra = sqrt(sum(xs2)*sigsq/delc) errb = sqrt(n*sigsq/delc) print "a = %.3f +- %.3f" % (a, erra) print "b = %.3f +- %.3f" % (b, errb) /// }}}Sage sadrži funkciju find_fit() koja radi otprilike to isto:
{{{id=444| var('a b') model(x) = a + b*x /// }}} {{{id=171| fit = find_fit(data, model, solution_dict=True); fit /// }}}Tu smo, kao i ranije kod rješavanja jednadžbi zatražili rješenje u obliku riječnika kako bismo ga lako uvrstili u model putem metode subs():
{{{id=517| list_plot(data, figsize=[3,3]) + plot(model(x).subs(fit), (x, 0, 5), color='red') /// }}}♦ Zadatak 2-7.4: Donjim podacima iz liste data2 prilagodite polinom trećeg reda, te funkciju $f(x)=a \sin(x)$ te nacrtajte sve na jednom grafu. Radi lakšeg snalaženja neka krivulje budu različite boje.
{{{id=487| data2 = [[0, 0.1], [0.3, 0.4], [1, 0.7], [1.8, 1], [2.5, 0.5], [4, -0.6], [5, -0.9]] /// }}}Nažalost, find_fit() ne daje greške parametara pravca $a$ i $b$. Za to trebamo koristiti naprednije funkcije za prilagodbu. Zgodna takva funkcija je leastsq() iz paketa scipy.optimize. Ona kao svoj prvi argument traži funkciju koju minimiziramo (kod metode najmanjih kvadrata to je razlika ordinate točke i vrijednosti funkcije u točki: $y_i - f(x_i)$. leastsq() sam radi kvadriranje te razlike!). Tu funkciju treba definirati kao python funkciju (vidi slijedeće poglavlje za detalje) što znači da mora imati sintaksu:
def ime_funkcije(argumenti):
...
return rezultat
gdje je tijelo funkcije uvučeno (To je python sintaktičko pravilo: ono što bi u C-u bilo unutar vitičastih zagrada ovdje se uvlači (tijelo funkcija, if-then naredbi, petlji). Obično se uvlači za širinu od 4 znaka.)
Za leastsq(), parametri koje prilagođavamo trebaju biti u listi (lista p[] dolje) , a drugi argument leastsq() su početne vrijednosti tih parametera.
Funkcija koju minimiziramo (dolje je zovemo distances()) treba vraćati listu odstupanja pravca od točaka.
Treći argument leastsq() je tuple s eventualnim dodatnim argumentima funkcije koju minimiziramo (prvi argument te funkcije je uvijek lista parametara koje prilagođavamo). Ovdje ćemo taj dodatni argument iskoristiti da proslijedimo toj funkciji same podatke (data) na koje prilagođujemo funkciju.
{{{id=454| def distances(p, data): return [(y - p[0] - p[1]*x) for x,y in data] /// }}}Po defaultu, leastsq() isto vraća samo tupl s listom prilagođenih vrijednosti parametara (bez grešaka) i zastavicom (flag) koja označava tip pronađenog rješenja ili eventualne probleme:
{{{id=481| scipy.optimize.leastsq(distances, [1.,1.], (data,)) /// }}}Ukoliko želimo i greške koristimo opcionalni argument full_output=True, koji nam daje i tzv. matricu kovarijanci (cov_x dolje) iz koje odredimo greške parametara na slijedeći način:
{{{id=456| p_final, cov_x, infodict, msg, ier = scipy.optimize.leastsq(distances, [1.,1.], (data,), full_output=True) vr = scipy.dot(distances(p_final, data), distances(p_final, data))/(len(data)-len(p_final)) p_errs = sqrt(scipy.diagonal(vr*cov_x)) print "p[0] =%8.3f +- %.4f" % (p_final[0], p_errs[0]) print "p[1] =%8.3f +- %.4f" % (p_final[1], p_errs[1]) /// }}}Gornji način određivanja greške parametara se može prihvatiti kao recept, no za one koji žele znati evo objašnjenja: (*) Kad se minimizira ispravno normalizirana funkcija (poput $\chi^2$ niže) onda dijagonalni elementi matrice kovarijanci izravno daju greške parametara. No, za običnu prilagodbu metodom najmanjih kvadrata, minimizira se kvadrat apsolutne greške pa je potrebno renormalizirati matricu kovarijanci varijancom mjerenja, a procjenitelj te varijance je upravo var=(suma kvadrata odstupanja)/(broj stupnjeva slobode). Inače, kad se tako vrijednost kvadrata odstupanja iskoristi za procjenu greške mjerenja ista se više ne može jednostavno iskoristiti za procjenu dobrote fita. (Vidi diskusiju u Bevingtonu ispod Eq. (6.15).)
leastsq() također omogućuje i prilagodbu složenijih funkcija. Npr. možemo se pitati trebamo li gornjoj funkciji dodati i kvadratni član:
{{{id=458| def distances(p, data): return [(y - p[0] - p[1]*x - p[2]*x^2) for x,y in data] /// }}} {{{id=459| p_final, cov_x, infodict, msg, ier = scipy.optimize.leastsq(distances, [1.,1.,1.], (data,), full_output=True) vr = scipy.dot(distances(p_final, data), distances(p_final, data))/(len(data)-len(p_final)) p_errs = sqrt(scipy.diagonal(vr*cov_x)) for k in range(len(p_final)): print "p[%i] =%8.3f +- %.4f" % (k, p_final[k], p_errs[k]) /// }}}Činjenica da je greška paramera p[2] znatno veća od vrijednosti samog parametra, sugerira da je kvadratni član suvišan. (Štoviše, vidimo da je i slobodni član p[0] problematične vrijednosti.)
Radi informacije spomenimo dva alternativna paketa za linearnu regresiju (dakle prilagođava se isključivo pravac). Prvi je OLS:
{{{id=261| load('http://www.phy.hr/~kkumer/sp/ols.py') /// }}}OLS traži svoje argumente u obliku numpy listi:
OLS nam daje i statistički test dobrote prilagodbe putem tzv F-testa, gdje je prilagodba dobra ukoliko je atribut Fpv (F-test p-value) manji od 0.01 ili barem 0.1.
{{{id=210| ft.Fpv # F-test probability /// }}}Druga mogućnost je poznati specijalizirani jezik za statističku analizu "R", koji je sadržan u Sage-u (premda je sučelje još uvijek relativno primitivno). Linearna regresija pomoću R-sučelja ide ovako:
{{{id=464| rx = r([x for x,y in data]) ry = r([y for x,y in data]) r.summary(r.lm(ry.tilde(rx))) /// }}}Ovdje gore imamo prvo "residuals", što su greške pojedinih točaka, zatim "coefficients", što su parametri pravca gdje imamo i statistički test signifikantnosti parametara (zadnji stupac Pr(>|t|) gdje je manja brojka bolja). Na kraju, zadnji broj (p-value) označava dobrotu prilagodbe određenu F-testom jednaku atributu Fpv OLS paketa gore.
Više o sučelju s R programskim jezikom vidi ovdje.
Često je potrebno kvantificirati dobrotu neke prilagodbe te odrediti da li predloženi teorijski model dobro opisuje mjerenja. Uzmimo slijedeći niz mjerenja s greškama $(x_i, y_i, \sigma_i\equiv\Delta y_i$):
{{{id=494| errdata = [[0.55, 0.74, 0.20], [0.9, 1.07, 0.12], [1.2, 1.28, 0.12], [1.5, 1.02, 0.12], [1.9, 1.63, 0.20]] xs = [x for x,y,err in errdata] ys = [y for x,y,err in errdata] errs = [err for x,y,err in errdata] /// }}}Za crtanje točaka s greškama (errorbars), koristimo pylabovu funkciju pylab.errorbar() gdje greške idu u opcionalni argument yerr:
{{{id=341| pylab.figure(figsize=[5,3]) pylab.errorbar(xs, ys, yerr=errs, linestyle='None') pylab.savefig('fig') /// }}}Promotrimo tri moguća teorijska opisa ovih mjerenja: potenciju $x^a$, logaritam $\log(a x)$ i obični pravac $a + b x$. (Zasad zanemarimo različite greške pojedinih mjerenja.)
{{{id=216| powmodel(x) = x^a logmodel(x) = log(a*x) linmodel(x) = a*x+b /// }}} {{{id=304| powfit=find_fit([[x, y] for x,y,err in errdata], powmodel, solution_dict=True); print powfit logfit=find_fit([[x, y] for x,y,err in errdata], logmodel, solution_dict=True); print logfit linfit=find_fit([[x, y] for x,y,err in errdata], linmodel, solution_dict=True); linfit /// }}} {{{id=348| xvals = pylab.linspace(0.4, 2) pylab.plot(xvals, [powmodel(x).subs(powfit).n() for x in xvals],'r') pylab.plot(xvals, [logmodel(x).subs(logfit).n() for x in xvals],'b-.') pylab.plot(xvals, [linmodel(x).subs(linfit).n() for x in xvals], color='black', linestyle=':') pylab.savefig('fig') /// }}}Teško je "od oka" procijeniti koji je opis najbolji, a još teže koji opisi su prihvatljivi a koji nisu. Standardna mjera dobrote fita (M. Požek, SOM, p. 128) je
$$ \chi^2 = \sum_i \frac{(y_i - f(x_i))^2}{\sigma_{i}^2} $$
gdje su $\sigma_i$ greške mjerenja. (U SOM skriptama $\sigma_i=\sqrt{y_i}$ jer je tamo riječ o mjerenju frekvencija.). Dobre su prilagodbe s $\chi^2 \approx df$, gdje je $df$ broj stupnjeva slobode (broj mjerenja minus broj slobodnih parametara funkcije koju prilagođavamo). Za točnije kvantificiranje dobrote prilagodbe gleda se integral statističke $\chi^2$ raspodjele s $df$ stupnjeva slobode od izračunatog $\chi^2$ do beskonačnosti koji se često zove $p$-vrijednost i koji ne bi trebao biti manji od 0.1 ili barem 0.01 ili hipotezu treba odbaciti. Taj integral je implementiran kao funkcija scipy.stats.chisqprob(chisq, df).
Definirajmo tri funkcije koje prilagođavamo ...
{{{id=126| def powfun(p, x): return x^p[0] def logfun(p, x): return log(p[0]*x) def linfun(p, x): return p[1]*x+p[0] /// }}}... i odgovarajuću funkciju udaljenosti (čije kvadrate će minimizirati leastsq()). Tu trebamo staviti i težinski faktor $1/\sigma_i$ za svaku točku kako bi točke s manjom greškom više utjecale na prilagodbu. Ovdje ćemo dodati još jedan argument putem kojeg ćemo prosljeđivati jednu od gornje tri funkcije:
{{{id=231| def dist(p, fun, data): return pylab.array([(y - fun(p, x))/err for (x,y,err) in data]) /// }}}Definirajmo i odgovarajući $\chi^2$:
{{{id=142| def chisq(p, fun, data): return sum( dist(p, fun, data)^2 ) /// }}}Nakon minimizacije ćemo testirati i zastavicu ier tako da ukoliko leastsq() nije uspješna dobijemo odgovarajuću poruku o grešci:
{{{id=232| fitfun = powfun ppow_final, cov_x, infodict, msg, ier = scipy.optimize.minpack.leastsq(dist, [1.], (fitfun, errdata), full_output=True) if ier not in (1,2,3,4): print " -----> No fit! <--------" print msg else: print "a =%8.4f +- %.4f" % (ppow_final, sqrt(cov_x[0,0])) chi2 = chisq([ppow_final], fitfun, errdata) df = len(errdata)-1 print "chi^2 = %.3f" % chi2 print "df = %i" % df print "p-vrijednost = %.4f" % scipy.stats.chisqprob(chi2, df) /// }}} {{{id=504| fitfun = logfun plog_final, cov_x, infodict, msg, ier = scipy.optimize.minpack.leastsq(dist, [1.], (fitfun, errdata), full_output=True) if ier not in (1,2,3,4): print " -----> No fit! <--------" print msg else: print "a =%8.4f +- %.4f" % (plog_final, sqrt(cov_x[0,0])) chi2 = chisq([plog_final], fitfun, errdata) df = len(errdata)-1 print "chi^2 = %.3f" % chi2 print "df = %i" % df print "p-vrijednost = %.4f" % scipy.stats.chisqprob(chi2, df) /// }}} {{{id=237| fitfun = linfun plin_final, cov_x, infodict, msg, ier = scipy.optimize.minpack.leastsq(dist, [1., 1.], (fitfun, errdata), full_output=True) if ier not in (1,2,3,4): print " -----> No fit! <--------" print msg else: parerrs = sqrt(scipy.diagonal(cov_x)) print "p[0] =%8.3f +- %.4f" % (plin_final[0], parerrs[0]) print "p[1] =%8.3f +- %.4f" % (plin_final[1], parerrs[1]) chi2 = chisq(plin_final, fitfun, errdata) df = len(errdata)-2 # dva parametra print "chi^2 = %.3f" % chi2 print "df = %i" % df print "p-vrijednost = %.4f" % scipy.stats.chisqprob(chi2, df) /// }}} {{{id=507| pylab.figure(figsize=[5,3]) pylab.errorbar(xs, ys, yerr=errs, linestyle='None') xvals = pylab.linspace(0.4, 2) pylab.plot(xvals, [powfun([ppow_final], x) for x in xvals],'r') pylab.plot(xvals, [logfun([plog_final], x) for x in xvals],'b-.') pylab.plot(xvals, [linfun(plin_final, x) for x in xvals], color='black', linestyle=':') pylab.savefig('fig') /// }}}Ova analiza sugerira da možemo odbaciti logaritamski model, dok su druga dva modela, u nedostaku boljih, prihvatljiva.
♦ Zadatak 2-7.5 Prilagodite funkciju $f(x) = a x$ donjim podacima. Odredite parametar $a$ i njegovu grešku te ocijenite dobrotu fita.
{{{id=253| xs2, ys2, errs2 = [range(2,25,2), [5.3,14.4,20.7,30.1,35.0,41.3,52.7,55.7,63.0,72.1,80.5,87.9], [1.5 for k in range(2,25,2)]] errdata2 = zip(xs2, ys2, errs2) /// }}}