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 [1]:
%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)
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[1]:
['33.194    0.351',
 '-3.246    0.351   10.990',
 '-3.246   -2.211   10.990    0.431',
 ' 1.206   -2.211    6.460    0.431    0.536',
 ' 1.206   -2.383    6.460    0.513    0.536   -0.008',
 ' 1.281   -2.383    6.298    0.513    0.586   -0.008   -0.004',
 ' 1.281   -4.186    6.298    2.214    0.586   -0.398   -0.004    0.025',
 ' 2.704   -4.186    0.834    2.214    3.778   -0.398   -0.589    0.025    0.033',
 ' 2.704   -5.205    0.834    3.833    3.778   -1.079   -0.589    0.129    0.033   -0.005']
In [2]:
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]

směrodatná odchylka koeficientů:
Out[2]:
[' 2.735    1.502',
 ' 0.622    0.227    0.140',
 ' 0.584    0.538    0.132    0.083',
 ' 0.338    0.249    0.217    0.038    0.025',
 ' 0.352    0.459    0.225    0.186    0.026    0.017',
 ' 0.431    0.479    0.545    0.194    0.155    0.017    0.011',
 ' 0.427    0.727    0.540    0.554    0.153    0.121    0.011    0.008',
 ' 0.386    0.579    0.820    0.441    0.426    0.096    0.075    0.006    0.004',
 ' 0.401    0.851    0.853    1.060    0.443    0.414    0.078    0.062    0.004    0.003']
In [3]:
s0=[((y-polyval(r[0],x))**2).sum() for r in res]
siges=sqrt(array(s0)/(20-1-ords))
print("redukov. chi2: \n",siges)

('redukov. chi2: \n', array([ 34.59892989,   5.21216215,   4.8776368 ,   2.24094107,
         2.31738384,   2.40354435,   2.36482779,   1.86413792,   1.91733631]))
In [4]:
from scipy import stats
print("5% kvantil Stud. rozdělení (s kles. počtem stupňů volnosti): \n",stats.t(19-ords).isf(0.05))

('5% kvantil Stud. rozd\xc4\x9blen\xc3\xad (s kles. po\xc4\x8dtem stup\xc5\x88\xc5\xaf volnosti): \n', array([ 1.73406361,  1.73960673,  1.74588368,  1.75305036,  1.76131014,
        1.7709334 ,  1.78228756,  1.79588482,  1.81246112]))

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]:
['12.135    0.234',
 ' 5.221    1.546   78.509',
 ' 5.556    4.106   83.540    5.183',
 ' 3.564    8.895   29.817   11.227   21.776',
 ' 3.428    5.193   28.675    2.764   20.942    0.454',
 ' 2.971    4.975   11.563    2.648    3.789    0.435    0.330',
 ' 2.996    5.756   11.662    3.997    3.822    3.303    0.332    3.273',
 ' 7.012    7.236    1.017    5.024    8.871    4.152    7.817    4.112    7.819',
 ' 6.742    6.115    0.977    3.615    8.529    2.605    7.515    2.090    7.516    1.691']
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,sigcom)
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í) - proto je menší než rezidua zahrnující šum

správným kritériem je validace - chi2 na testovacím vzorku:

opakovaná simulace = validace

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]:
array([ 1238.85841607,    17.42292297,    17.68858138,     1.41782944,
           1.80194635,     2.20549544,     2.62870594,     3.26606789,
           3.93339854])
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]:
<matplotlib.legend.Legend at 0xb92ae10>
In [17]:
apars=[test_poly(x,ytrue,2,rep=0) for k in range(100)] #zaznam hodnot nafitovanych parametru
In [18]:
def rep_pars(iord=4):
    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]:
array([ 0.08969589, -0.34124228,  5.28565051,  0.24595817,  1.54955129,
       -0.07131472, -0.13749907,  0.06867428])

kovarianční matice dobře odpovídá rozptylu parametrů

In [19]:
zot6[1]/zot6[2]
Out[19]:
array([ 0.98661703,  0.95119805,  0.96080738,  0.87037112,  0.95952248,
        0.85984801,  0.95468107,  0.86476922])

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]:
['73.504  3.083',
 '-5.852  3.083 73.102',
 '-5.852 -0.720 73.102  2.004',
 ' 0.024 -0.720 12.865  2.004  8.473',
 ' 0.024 -0.416 12.865  0.473  8.473 -0.019',
 ' 0.090 -0.416  5.286  0.473  1.550 -0.019 -0.137',
 ' 0.090 -0.341  5.286  0.246  1.550 -0.071 -0.137  0.069']