Fitování 2 parametrů polynomiálního modelu¶

nejistoty z kovarianční matice
sestrojení mapy hodnot $\chi^2$

In [1]:

%pylab inline
x=r_[0:2:.05]
pars=r_[2,3]
y=polyval(pars,x)+normal(0,0.1,size=x.shape)
res1,cov1=polyfit(x,y,1,cov=True)
#cov1=100
err1=sqrt(cov1.diagonal())

Populating the interactive namespace from numpy and matplotlib

In [2]:

plot(x,y-3-2*x,'.')
plot(x,y-res1[1]-res1[0]*x,'v')
grid()
legend(["skut. rez.","fit. rez."])
res1

Out[2]:

array([ 1.96981566,  3.04829809])

In [55]:

# odhad chyby 
std(y-3-2*x),std(y-res1[1]-res1[0]*x)

Out[55]:

(0.094728002242672735, 0.090203226761920607)

namapování v prostoru parametrů

In [5]:

grx,gry=mgrid[-1:1:50j,-1:1:50j]
pos2=res1.reshape(1,1,2)+array([2*err1[0]*grx,2*err1[1]*gry]).swapaxes(0,2) #mozne dvojice parametru
chi2min=sum((y-polyval(res1,x))**2*100)
print('nejistoty parametru: %s'%list(err1)+'\t chi2 v minimu: %.4f'%chi2min)

nejistoty parametru: [0.026047569764186713, 0.029512597726795382]	 chi2 v minimu: 32.5465

In [9]:

chi2=array([[sum((y-polyval(p,x))**2)*100 for p in q] for q in pos2])
contour(chi2,extent=[pos2[0][0][0],pos2[0][-1][0],pos2[0][0][1],pos2[-1][0][1]],levels=r_[chi2min:chi2min+5:1],)
axhline(res1[1]+err1[1])
axvline(res1[0]+err1[0])
axhline(res1[1]-err1[1])
axvline(res1[0]-err1[0])
colorbar()
plot(pars[0],pars[1],'d')
text(pars[0]+0.01,pars[1],'prava hodnota')
xlabel('lin.clen')
ylabel('konst.clen')

Out[9]:

<matplotlib.text.Text at 0xb0be99cc>

hodnota sumy reziduí

$$S - S_0= (P-P_0)^T H (P-P_0)$$

v rozsahu $\pm 2 \sigma$

In [34]:

hess=inv(cov1)
wee=[[dot(dot(p,hess),p) for p in q] for q in array([2*err1[0]*grx,2*err1[1]*gry]).swapaxes(0,2)]
contour(wee,extent=[pos2[0][0][0],pos2[0][-1][0],pos2[0][0][1],pos2[-1][0][1]])
colorbar()
hess

Out[34]:

array([[ 5679.87542352,  4313.82943559],
       [ 4313.82943559,  4424.44044675]])

In [110]:

print('Hessian:')
hess

Out[110]:

array([[ 5259.94549698,  3994.89531417],
       [ 3994.89531417,  4097.32852735]])

In [57]:

xax=r_[pos2[0][0][0]:pos2[0][-1][0]:1j*len(wee[0])]
plot(xax,wee[25])
legend(['func. from Hess.'],loc=3)
pyplot.twinx(gca())
plot(xax,chi2[25],'r')
ylim(chi2min,None)
legend(['chi^2'],loc=4)
#pyplot.twinx?

Out[57]:

<matplotlib.legend.Legend at 0xaf407e4c>

In [51]:

q1,q0=array([polyfit(wee[i],chi2[i],1) for i in range(50)]).transpose()

In [56]:

q0.mean(),q1.mean()

Out[56]:

(0.38002970696801269, 0.011814447074971667)

In [41]:

#wee[chi2==chi2min]=1
imshow(chi2/(wee+chi2min),origin="lower")#,vmin=0.0117,vmax=0.0119)
title('rozdil chi2/Hess')
colorbar()

Out[41]:

<matplotlib.colorbar.Colorbar instance at 0xaf69a3ac>

výsledky při 2000 opakování

In [43]:

#rozdeleni nafitovanych parametru
y0=polyval(pars,x)
fit1=[polyfit(x,y0+normal(0,0.1,size=x.shape),1,cov=True) for i in range(2000)]
set1=array([s[0] for s in fit1]).transpose()
plot(set1[0],set1[1],'.')
grid()

In [46]:

#rozdeleni nejistot (z kovariancni matice)
eset1=array([sqrt(f[1].diagonal()) for f in fit1]).transpose()
q1=hist(eset1[0])

In [58]:

# korelace (vychazi vzdy stejne)
cov1=array([f[1][0,1] for f in fit1])
cor1=cov1/eset1[0]/eset1[1]
cor1.max()-cor1.min()

Out[58]:

5.5511151231257827e-16

In [59]:

set1[0].mean(),set1[0].std(),err1[0],eset1[0].mean()

Out[59]:

(2.0009311454456427,
 0.027880267203110549,
 0.026047569764186713,
 0.028063379961028374)

In [29]:

set1[1].mean(),set1[1].std(),err1[1]

Out[29]:

(2.9991860899454021, 0.029682261213690064, 0.029512597726795382)

In [32]:

dot(set1[0]-set1[0].mean(),set1[1]-set1[1].mean())/set1[0].std()/set1[1].std()/2000.

Out[32]:

-0.85498035205646727

In [33]:

cov1[0][1]/err1[0]/err1[1]

Out[33]:

-0.86052677419934109

In [14]:

dit1=set1-pars[:,newaxis]
repx=(dit1*dot(hess,dit1)).sum(0)

repx je vzdálenost od pravé hodnoty měřená metrikou Hessianů
rozdělení odpovídá $\chi^2_{N-2}$

In [15]:

rechi2=hist(sqrt(repx),20)

In [24]:

from scipy import stats
print("podil do vzdalenosti 1sigma: %.4f (teor. %.4f)"%(sum(repx<1)/2000.,stats.chi2(2).cdf(1)))
print("podil do 1sigma v 1 promenne: %.4f (teor. %.4f)"%(sum(abs(dit1[0])/err1[0]<1)/2000.,(stats.norm().cdf(1)-0.5)*2))
print("podil do 1sigma ve 2 promennych: %.4f (teor. %.4f)"%(sum((abs(dit1[0])/err1[0]<1)*(abs(dit1[1])/err1[1]<1))/2000.,(stats.norm().cdf(1)-0.5)**2*4))

podil do vzdalenosti 1sigma: 0.3875 (teor. 0.3935)
podil do 1sigma v 1 promenne: 0.6755 (teor. 0.6827)
podil do 1sigma ve 2 promennych: 0.5800 (teor. 0.4661)