pokračování fitování polynomy různého stupně¶
pravé hodnoty odpovídají polynomu 4. stupně (dominantně kvadratický)
testujeme modely pro polynomy 1. - 10. stupně
In [11]:
%matplotlib inline
from numpy import *
from matplotlib import pyplot as pl
x=r_[-3:3:20j]
sigy=3.
tres=[0.5,0.2,7,-0.5,0] #skutecne parametry
ytrue=polyval(tres,x)
pl.plot(x,ytrue,'k')
y=ytrue+random.normal(0,sigy,size=x.shape)
pl.plot(x,y,'*')
ords=arange(1,10)
In [12]:
res=[polyfit(x,y,i,cov=True) for i in ords]
#[[round(p,3) for p in r[0][::-1]] for r in res]
[" ".join(["%6.3f"%p for p in r[0][::-1]]) for r in res]
Out[12]:
In [3]:
errs=[[round(p,5) for p in sqrt(r[1].diagonal()[::-1])/sigy] for r in res]
print('směrodatná odchylka koeficientů:')
[" ".join(["%6.3f"%p for p in e]) for e in errs]
Out[3]:
test sumy čtverců reziduí¶
In [4]:
smanu=[]
for i in range(1000):
y=ytrue+random.normal(0,sigy,size=x.shape)
res=[polyfit(x,y,i,cov=False) for i in ords]
smanu.append([((y-polyval(r,x))**2).sum() for r in res] )
smanu=array(smanu)
#siges=sqrt(array(s0)/(20-1-ords))
In [7]:
ok=pl.hist(smanu[:,5]/sigy**2,20,alpha=0.6)
ok=pl.hist(smanu[:,8]/sigy**2,20,alpha=0.6)
[stats.chi2.fit(sm/sigy**2, floc=0, fscale=1)[0] for sm in smanu.T]
Out[7]:
In [8]:
# redukovany chi2
msig=smanu.mean(0)/sigy**2/(20-1-ords)
msig
Out[8]:
In [9]:
# rozdeleni rozdilu rezidualni sumy
dif=(smanu[:,5]-smanu[:,8])/sigy**2
ok=pl.hist(dif,20)
ps0=stats.chi2.fit(dif, floc=0, fscale=1)[0]
voo=stats.chi2(ps0).pdf(ok[1])
voo*=sum(ok[0])/sum(voo)
pl.plot(ok[1],voo)
ps0
Out[9]:
In [10]:
from scipy import stats
print("5% kvantil Stud. rozdělení (s kles. poctem stupnu volnosti): \n",stats.t(19-ords).isf(0.05))
statistické významnosti¶
(t-hodnoty pro nulovou hypotézu) jednotlivých parametrů
pro dané měření mohou vycházet významné hodnoty i u vyšších momentů
In [5]:
[" ".join(["%6.3f"%p for p in abs(res[i][0][::-1]/array(errs[i]))]) for i in range(len(res))]
Out[5]:
In [7]:
strue=[((ytrue-polyval(r[0],x))**2).sum() for r in res]
#sigcom=sqrt(array(s0)/(20-1-ords)) # redukovany chi2
pl.semilogy(ords,siges)
pl.semilogy(ords,sqrt(array(strue)/(20-1-ords)))
pl.legend(['red. chi2','true chi2'])
pl.grid()
- (i redukované) chi2 polynomu soustavně klesá - vyšší polynomy "kopírují šum" v modelových datech
- pravé chi2 je rozdíl polynomu daného stupně od nezašuměné křivky (střední hodnota velkého množství opakování) -- obvykle je menší než rezidua zahrnující šum
správným kritériem je validace - chi2 na testovacím vzorku...
opakovaná simulace = validace¶
numerický experiment můžeme snadno opakovat
In [12]:
from numpy.random import normal
def test_poly(x,ytrue,sigy,ords=r_[1:10],rep=1):
'''podle rep:
1: rozdil skutecne funkce a modelu
2: rozdil nove sady dat a modelu
'''
y=ytrue+normal(0,sigy,size=x.shape)
res=[polyfit(x,y,i,cov=True) for i in ords]
if rep==1:
scom=[((ytrue-polyval(r[0],x))**2).sum() for r in res]
return array(scom)/(len(x)-1-ords)
if rep==2:
ycom=ytrue+normal(0,sigy,size=x.shape)
scom=[((ycom-polyval(r[0],x))**2).sum() for r in res]
return array(scom)/(len(x)-1-ords)
return res
In [11]:
# provedeme 100 opakování testu
fullmat=r_[[test_poly(x,ytrue,2,rep=1) for k in range(100)]]
fullmat.mean(0)
Out[11]:
In [13]:
fullmat2=r_[[test_poly(x,ytrue,2,rep=2) for k in range(100)]]
In [16]:
pl.errorbar(ords[1:],fullmat.mean(0)[1:],fullmat.std(0)[1:]/10.)
pl.errorbar(ords[1:],fullmat2.mean(0)[1:],fullmat2.std(0)[1:]/10.)
pl.xlabel("stupen polynomu")
pl.legend(["gen. chyba","validace"])
Out[16]:
záznam hodnot nafitovaných parametrů¶
In [17]:
apars=[test_poly(x,ytrue,2,rep=0) for k in range(100)]
In [18]:
def rep_pars(iord=4):
# vraci prumerne hodnoty fitovanych parametru, teoret. nejistoty a rozptyl parametru
zmat=r_[[a[iord][0] for a in apars]]
cmat=r_[[a[iord][1].diagonal() for a in apars]]
zrep=array([zmat.mean(0),zmat.std(0),sqrt(cmat.mean(0)),sqrt(cmat).std(0)])
return zrep[:,::-1]
zot6=rep_pars(6)
zot6[0]/zot6[1] # podil stredni hodnoty a nejistoty
Out[18]:
kovarianční matice dobře odpovídá rozptylu parametrů
In [19]:
zot6[1]/zot6[2]
Out[19]:
tabulka "významnosti" parametrů¶
In [20]:
zotx=[rep_pars(i) for i in range(7)]
[' '.join(['%6.3f'%ww for ww in zz[0]/zz[1]]) for zz in zotx]
Out[20]: