2-7-Statistika_i_obrada_podataka-S
system:sage


<h2>Statistika (scipy.stats)</h2>

<p>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.</p>

{{{id=269|
import scipy
import scipy.stats
import pylab
///
}}}

<p>Kreirajmo prvo (numpy) listu brojeva. Numpy liste imaju metode za izračun osnovnih statističkih veličina poput srednje vrijednosti, varijance, standardne devijacije itd.</p>

{{{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()
///
}}}

<p>(Ovdje smo za formatiranje ispisa primjenili operator %, te standardne <a href="http://docs.python.org/release/2.5.2/lib/typesseq-strings.html">stringove za formatiranje</a> (poznate iz sprintf() C funkcije): Npr. %5.3f formatira realni broj na 5 mjesta ukupno i 3 mjesta iza decimalne točke.)</p>
<p>Obične liste nemaju gornje statističke metode ...</p>

{{{id=388|
t1 = [3, 5, 7, 9, 11]
t1.mean()
///
}}}

<p>... pa ukoliko nas zanima npr. srednja vrijednost brojeva u običnoj listi treba je prvo konvertirati u numpy listu funkcijom pylab.array():</p>

{{{id=389|
t1num = pylab.array(t1)
t1num.mean()
///
}}}

<p>... ili možemo primijeniti istoimenu Sage funkciju (dakle ne metodu):</p>

{{{id=427|
mean(t1)
///
}}}

<p>scipy.stats paket ima implementirane brojne statističke raspodjele:</p>
<ul>
<li><strong><span style="font-family: courier new,courier;">norm</span></strong> - normalna (Gaussova) raspodjela</li>
<li><strong><span style="font-family: courier new,courier;">poisson</span></strong> - Poissonova raspodjela</li>
<li><strong><span style="font-family: courier new,courier;">gamma</span></strong> - gamma raspodjela</li>
<li><strong><span style="font-family: courier new,courier;">chi2</span></strong> - $\chi^2$ raspodjela</li>
<li><strong><span style="font-family: courier new,courier;">t</span></strong> - studentova t raspodjela</li>
<li>...</li>
</ul>
<p>Svaka od tih raspodjela ima svoje parametre (npr. srednja vrijednost i standardna devijacija), koji se najče&scaron;će specificiraju putem opcionalnih argumenata loc i scale, čije točno značenje saznamo iz standardne dokumentacije ili isprintavanjem specijalnog konciznog atributa 'extradoc':</p>

{{{id=407|
print scipy.stats.norm.extradoc
///
}}}

<p>Najvažnije metode distribucija su:</p>
<ul>
<li><strong><span style="font-family: courier new,courier;">pdf()</span></strong>&nbsp; - sama distribucija (<em>probability distribution function</em>) $f(x)$</li>
<li><strong><span style="font-family: courier new,courier;">cdf()</span></strong>&nbsp; - integral distribucije&nbsp; (<em>cumulative distribution function</em>) $\Phi(x)=\int_{-\infty}^{x} f(y) dy$</li>
<li><strong><span style="font-family: courier new,courier;">rvs()</span></strong>&nbsp; - slučajni brojevi (<em>random variates</em>) koji slijede distribuciju</li>
</ul>
<p>Tako je npr. vrijednost u $x=0$ normalne distribucije sa srednjom vrijednosti $\mu = 5$ i&nbsp; standardnom devijacijom $\sigma = 1.5$, $f(x=0, \mu=5, \sigma=1.5)$ dan s:</p>

{{{id=390|
scipy.stats.norm.pdf(0, loc=5., scale=1.5)
///
}}}

<p>Skica distribucije&nbsp; i njenog integrala:</p>

{{{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')
///
}}}

<p>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:</p>

{{{id=413|
scipy.stats.norm.cdf(1) - scipy.stats.norm.cdf(-1)   # defaulti su mu=0, sigma=1
///
}}}

<p>(To je čuvenih 68% vjerojatnosti.)</p>
<p><strong>&diams;&nbsp;Zadatak 2-7.1</strong>: ("Primjer 2" is SOM skripti M. Požeka)&nbsp; 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?</p>

<p>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:</p>

{{{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()
///
}}}

<p>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:</p>

{{{id=369|
pylab.figure(figsize=[3.5, 2])
pylab.hist(pts, bins=20, normed=True)
pylab.savefig('tmp')
///
}}}

<p>Ovo je bilo za jedan uzorak. Listu npr. standardnih devijacija za 5 nezavisnih uzoraka možemo konstruirati obuhvaćanjem liste uz kori&scaron;tenje funkcije range() koja kreira listu cijelih brojeva:</p>

{{{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)]
///
}}}

<p><strong>&diams;&nbsp;Zadatak 2-7.2</strong>:&nbsp;Koja je vjerojatnost da mjerenje slučajne varijable koja slijedi normalnu razdiobu sa srednjom vrijedno&scaron;ću $\mu$ i standardnom devijacijom $\sigma$ odstupi od srednje vrijednosti za vi&scaron;e od $3\sigma$?</p>

<p>Važna raspodjela je tzv. $\chi^2$<span id="cell_outer_166"><span><span></span></span>-raspodjela</span></p>

{{{id=418|
print scipy.stats.chi2.extradoc
///
}}}

<p>Usporedi sa SOM&nbsp; (Požek p. 80):</p>
<p>$$ P(\chi^2, \nu) = \frac{(\chi^2)^{\nu/2-1}e^{-\chi^2/2}}{2^{\nu/2} \Gamma(\nu/2)}$$</p>
<p>Vidimo da $\chi^2$ raspodjela, pored argumenta x, ima samo jedan parametar: $\nu$ = df = broj stupnjeva slobode (<em>degrees of freedom</em>). Naziv tog parametra će biti jasniji niže kod analize prilagodbi funkcija točkama.</p>
<p>Srednja vrijednost $\chi^2$ raspodjele jednaka je broju stupnjeva slobode:</p>

{{{id=420|
scipy.stats.chi2.rvs(3, size=1e5).mean()
///
}}}

<p>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$:</p>

{{{id=159|
[1-scipy.stats.chi2.cdf(eps^2,1) for eps in [1,2,3]]
///
}}}

<p><strong>&diams;&nbsp;Zadatak 2-7.3</strong>: Sredi&scaron;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.</p>
<ol>
<li>Nacrtajte tu raspodjelu za df=3 stupnja slobode i uvjerite se da je asimetrična oko srednje vrijednosti.</li>
<li>Uzmite jedan uzorak od milijun točaka te raspodjele i uvjerite se da je standardna devijacija $\sigma=\sqrt{2\cdot df}$</li>
<li>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.</li>
<li>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}$</li>
</ol> 

<h2>(*) Statistika &nbsp;("čisti" Sage)</h2>
<p>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.</p>

{{{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]]
///
}}}

<h2>Prilagodba funkcije podacima</h2>

<p>Često je potrebno neke podatke dobivene npr. mjerenjima grafički prikazati i približno opisati nekom funkcijom (tzv. prilagodba ili "fitanje").</p>

{{{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])
///
}}}

<p>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:</p>

{{{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)
///
}}}

<p>Sage sadrži funkciju find_fit() koja radi otprilike to isto:</p>

{{{id=444|
var('a b')
model(x) = a + b*x
///
}}}

{{{id=171|
fit = find_fit(data, model, solution_dict=True); fit
///
}}}

<p>Tu smo, kao i ranije kod rje&scaron;avanja jednadžbi zatražili rje&scaron;enje u obliku riječnika kako bismo ga lako uvrstili u model putem metode subs():</p>

{{{id=517|
list_plot(data, figsize=[3,3]) + plot(model(x).subs(fit), (x, 0, 5), color='red')
///
}}}

<p><span id="cell_outer_11"><span id="cell_outer_10"><span id="cell_outer_5">&diams; </span></span><strong>Zadatak 2-7.4</strong>: </span>Donjim podacima iz liste data2 prilagodite polinom trećeg reda, te funkciju $f(x)=a&nbsp; \sin(x)$ te nacrtajte sve na jednom grafu. Radi lak&scaron;eg snalaženja neka krivulje budu različite boje.</p>

{{{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]]
///
}}}

<p>Nažalost, find_fit() ne daje gre&scaron;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 <strong>prvi argument</strong> traži funkciju koju minimiziramo (kod metode najmanjih kvadrata to je razlika ordinate točke i vrijednosti funkcije u točki:&nbsp; $y_i - f(x_i)$. leastsq() sam radi kvadriranje te razlike!). Tu funkciju treba definirati kao python funkciju (vidi slijedeće poglavlje za detalje) &scaron;to znači da mora imati sintaksu:</p>
<p>def&nbsp; ime_funkcije(argumenti):</p>
<p>&nbsp;&nbsp;&nbsp; ...</p>
<p>&nbsp;&nbsp; return rezultat</p>
<p>gdje je tijelo funkcije uvučeno (To je python sintaktičko pravilo: ono &scaron;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 &scaron;irinu od 4 znaka.)</p>
<p>Za leastsq(), parametri koje prilagođavamo trebaju biti u listi (lista p[] dolje) , a <strong>drugi argument</strong> leastsq() su početne vrijednosti tih parametera.</p>
<p>Funkcija koju minimiziramo (dolje je zovemo distances()) treba vraćati listu odstupanja pravca od točaka.</p>
<p><strong>Treći argument</strong> 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.</p>

{{{id=454|
def distances(p, data):
    return [(y - p[0] - p[1]*x) for x,y in data]
///
}}}

<p>Po defaultu, leastsq() isto vraća samo tupl s listom prilagođenih vrijednosti parametara (bez gre&scaron;aka) i zastavicom (flag) koja označava tip pronađenog rje&scaron;enja ili eventualne probleme:</p>

{{{id=481|
scipy.optimize.leastsq(distances, [1.,1.],  (data,))
///
}}}

<p>Ukoliko želimo i gre&scaron;ke koristimo opcionalni argument full_output=True, koji nam daje i tzv. matricu kovarijanci (cov_x dolje) iz koje odredimo gre&scaron;ke parametara na slijedeći način:</p>

{{{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])
///
}}}

<p>Gornji način određivanja gre&scaron;ke parametara se može prihvatiti kao recept, no za one koji žele znati evo obja&scaron;njenja: (*) Kad se minimizira ispravno normalizirana funkcija (poput $\chi^2$ niže) onda dijagonalni elementi matrice kovarijanci izravno daju gre&scaron;ke parametara. No, za običnu prilagodbu metodom najmanjih kvadrata, minimizira se kvadrat <em>apsolutne</em> gre&scaron;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&scaron;ke mjerenja ista se vi&scaron;e ne može jednostavno iskoristiti za procjenu dobrote fita. (Vidi diskusiju u Bevingtonu ispod Eq. (6.15).)</p>
<p>leastsq() također omogućuje i prilagodbu složenijih funkcija. Npr. možemo se pitati trebamo li gornjoj funkciji dodati i kvadratni član:</p>

{{{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])
///
}}}

<p>Činjenica da je gre&scaron;ka paramera p[2] znatno veća od vrijednosti samog parametra, sugerira da je kvadratni član suvi&scaron;an. (&Scaron;tovi&scaron;e, vidimo da je i slobodni član p[0] problematične vrijednosti.)</p>
<h2>(*) Prilagodba funkcije podacima - alternativne metode</h2>
<p>Radi informacije spomenimo dva alternativna paketa za linearnu regresiju (dakle prilagođava se isključivo pravac). Prvi je OLS:</p>

{{{id=261|
load('http://www.phy.hr/~kkumer/sp/ols.py')
///
}}}

<p><span style="font-size: medium;">OLS traži svoje argumente u obliku numpy listi:<br /></span></p>

{{{id=467|
xs = scipy.array([x for x,y in data])
ys = scipy.array([y for x,y in data])
///
}}}

{{{id=466|
ft = ols(ys, xs, 'y', ['x1'])   # ols() kreira 'fit' objekt ft
///
}}}

{{{id=213|
ft.b # parametri pravca ...
///
}}}

{{{id=214|
ft.se # ... i njihove greške
///
}}}

<p><span style="font-size: small;"><br /></span></p>

<p>OLS nam daje i statistički test dobrote prilagodbe putem tzv F-testa, gdje je prilagodba dobra ukoliko je atribut Fpv (<em>F-test p-value</em>) manji od 0.01 ili barem 0.1.</p>

{{{id=210|
ft.Fpv  # F-test probability
///
}}}

<p><span id="cell_outer_21"><span id="cell_outer_10"><span id="cell_outer_5"><br /></span></span></span></p>

<p>Druga mogućnost je poznati specijalizirani jezik za statističku analizu "R", koji je sadržan u Sage-u (premda je sučelje jo&scaron; uvijek relativno primitivno). Linearna regresija pomoću R-sučelja ide ovako:</p>

{{{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)))
///
}}}

<p>Ovdje gore imamo prvo "residuals", &scaron;to su gre&scaron;ke pojedinih točaka, zatim "coefficients", &scaron;to su parametri pravca gdje imamo i statistički test signifikantnosti parametara (zadnji stupac Pr(&gt;|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.</p>
<p>Vi&scaron;e o sučelju s R programskim jezikom vidi <a href="https://localhost:8000/src/interfaces/r.py">ovdje</a>.</p>
<h3>Testiranje hipoteze</h3>

<p>Č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&scaron;kama $(x_i, y_i, \sigma_i\equiv\Delta y_i$):</p>

{{{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]
///
}}}

<p>Za crtanje točaka s gre&scaron;kama (errorbars), koristimo pylabovu funkciju pylab.errorbar() gdje gre&scaron;ke idu u opcionalni argument yerr:</p>

{{{id=341|
pylab.figure(figsize=[5,3])
pylab.errorbar(xs, ys, yerr=errs, linestyle='None')
pylab.savefig('fig')
///
}}}

<p><span style="font-size: medium;">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&scaron;ke pojedinih mjerenja.)</span></p>

{{{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')
///
}}}

<p>Te&scaron;ko je "od oka" procijeniti koji je opis najbolji, a jo&scaron; teže koji opisi su prihvatljivi a koji nisu. Standardna mjera dobrote fita (M. Požek, SOM, p. 128) je</p>
<p>$$ \chi^2 = \sum_i \frac{(y_i - f(x_i))^2}{\sigma_{i}^2} $$</p>
<p>gdje su $\sigma_i$ gre&scaron;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&nbsp; 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).</p>
<p>Definirajmo tri funkcije koje prilagođavamo ...</p>

{{{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]
///
}}}

<p>... 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&scaron;kom vi&scaron;e utjecale na prilagodbu. Ovdje ćemo dodati jo&scaron; jedan argument putem kojeg ćemo prosljeđivati jednu od gornje tri funkcije:</p>

{{{id=231|
def dist(p, fun, data):
    return pylab.array([(y - fun(p, x))/err for (x,y,err) in data])
///
}}}

<p>Definirajmo i odgovarajući $\chi^2$:</p>

{{{id=142|
def chisq(p, fun, data):
    return sum( dist(p, fun, data)^2 )
///
}}}

<p>Nakon minimizacije ćemo testirati i zastavicu ier tako da ukoliko leastsq() nije uspje&scaron;na dobijemo odgovarajuću poruku o gre&scaron;ci:</p>

{{{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')
///
}}}

<p>Ova analiza sugerira da možemo odbaciti logaritamski model, dok su druga dva modela, u nedostaku boljih, prihvatljiva.</p>

<p><span id="cell_outer_21"><span id="cell_outer_10"><span id="cell_outer_5">&diams; </span></span><strong>Zadatak 2-7</strong>.<strong>5 </strong></span>Prilagodite funkciju $f(x) = a x$ donjim podacima. Odredite parametar $a$ i njegovu gre&scaron;ku te ocijenite dobrotu fita.</p>

{{{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)
///
}}}