data se snažíme modelovat funkcí
$$c\sin(a-a_0)+ b =c \sin(a)\cos(a_0)-c \cos(a)\sin(a_0) + b$$tedy 3 parametrický model bude mít odpovídající vektor parametrů
$$p=[b, c \cos(a_0), c \sin(a_0)]$$a amplitudu resp. nulový polarizační úhel budeme rekonstruovat jako
$c=\sqrt(p_2^2+p_3^2),\ a_0=\arctan(p3/p2)$
%matplotlib inline
from matplotlib import pyplot as pl
import numpy as np
#ve kterych bodech merime
xval=np.r_[0:10:50j]
#modelová matice
mat_M=np.array([np.ones(shape=len(xval)), np.cos(xval), np.sin(xval)]).T
hessian=mat_M.T.dot(mat_M)
print("hessian")
hessian
#kovarianci matici ziskame jako inverzi
covar=np.linalg.inv(hessian)
covar
teď nasimulujeme nějaká data: nejprve předpokládejme konstantní nejistotu měřených hodnot $\sigma$
a0,b,c=30*np.pi/180.,50,20
sigma=5
pars0=np.array([b,c*np.cos(a0),c*np.sin(a0)]) #hodnoty parametru prepocitane podle vzorce 2
yteor=mat_M.dot(pars0)
yval=yteor+np.random.normal(0,sigma,size=len(xval))
pl.plot(xval,yteor,label="teor. prubeh")
pl.plot(xval,yval,"v",label="mereni")
pl.xlabel("xval")
pl.ylabel("yval")
pl.legend(loc=0)
# jednoduchy vzorecek
# p = H^(-1) M^T Y
pars1=covar.dot(mat_M.T.dot(yval))
pars1
ymodel=mat_M.dot(pars1)
pl.plot(xval,ymodel,'r',label="fitovany. prubeh")
pl.plot(xval,yval,"v",label="mereni")
pl.xlabel("xval")
pl.ylabel("yval")
pl.legend(loc=0)
pl.plot(xval,ymodel-yteor)
pl.title("odchylka od praveho modelu")
#suma rezidui
s0=sum((yval-ymodel)**2)
rs0=s0/(len(xval)-3) #3 urcovane parametry
print("reduk. chi^2: %.2f vs. skutecna sigma^2: %.2f"%(rs0,sigma**2))
errs=np.sqrt(covar.diagonal())
corel=covar/errs.reshape(1,3)/errs.reshape(3,1)
print("korelacni matice parametru")
corel
errs*=np.sqrt(rs0)
print("nejistoty parametru:"+str(errs))
na výpočet parametru $c$ aplikujeme šíření nejistot $\sigma_2, \sigma_3$ parametrů $p_2, p_3$
$$\frac{\partial c}{\partial p_2}=\frac{p_2}{\sqrt(p_2^2+p_3^2)}$$$$\frac{\partial c}{\partial p_3}=\frac{p_3}{\sqrt(p_2^2+p_3^2)}$$V případě nekorelovaných parametrů $p_2, p_3$ by platilo (v rozvoji do prvního řádu)
$$\sigma_c^2 = \left( \frac{\partial c}{\partial p_2}\right)^2 \sigma_2^2 + \left( \frac{\partial c}{\partial p_3}\right)^2 \sigma_2^3$$Zde ale máme nezanedbatelnou antikorelaci (-25%), tedy je nutno ještě přidat člen
$$\sigma_c^2 = \left( \frac{\partial c}{\partial p_2}\right)^2 \sigma_2^2 + \left( \frac{\partial c}{\partial p_3}\right)^2 \sigma_3^2 + \frac{\partial c}{\partial p_2} \sigma_2 \frac{\partial c}{\partial p_3} \sigma_3\ \rho_{23} $$# konkretne
c=np.sqrt(pars1[1]**2+pars1[2]**2)
sigma2_c =(pars1[1]/c)**2 * errs[1]**2
sigma2_c+=(pars1[2]/c)**2 * errs[2]**2
sigma2_c+=pars1[1]*pars1[2]/c**2 * errs[1]*errs[2]*corel[1,2]
print("parametr c: %.4f +- %.4f"%(c,np.sqrt(sigma2_c)))
sigma2x_c =(pars1[1]/c)**2 * errs[1]**2
sigma2x_c+=(pars1[2]/c)**2 * errs[2]**2
print("nejistota c bez zapocteni korelace : %.4f"%np.sqrt(sigma2x_c))
$y$ budou data s Poissonovským rozdělením s rozptylem $\sigma^2_y=y_{teor}$
Poněvadž ale $y_{teor}$ při analýze neznáme, budeme dosazovat experimentální hodnoty $y_{val}$
yval2=np.random.poisson(yteor)
mat_W=np.linalg.inv(yval2*np.eye(len(xval))) #vahova matice jako inverze kovarianci merenych hodnot
hessian2=mat_M.T.dot(mat_W.dot(mat_M))
covar2=np.linalg.inv(hessian2)
pars2=covar2.dot(mat_M.T.dot(mat_W.dot(yval2)))
ymodel2=mat_M.dot(pars2)
pl.plot(xval,ymodel2,'r',label="fitovany. prubeh")
pl.plot(xval,yval2,"v",label="mereni")
pl.xlabel("xval")
pl.ylabel("yval")
pl.legend(loc=0)
pars2
s02=sum((yval2-ymodel2)**2)
rs02=s02/(len(xval)-3) #3 urcovane parametry
print("reduk. chi^2: %.2f vs. \"prumerna\" sigma^2: %.2f (%.2f)"%(rs02,yval2.mean(),(np.sqrt(yval2).mean())**2))
minimalizovaná hodnota je ale ve skutečnosti $(Y-Y_m)^T W (Y-Y_m)$, po vydělení počtem stupňů volnosti (očekávaná hodnota $\chi^2$) máme dostat hodnotu blízkou 1
s02m=(yval2-ymodel2).T.dot(mat_W.dot(yval2-ymodel2))
rs02m=s02m/(len(xval)-3)
rs02m
errs2=np.sqrt(covar2.diagonal())
corel2=covar2/errs2.reshape(1,3)/errs2.reshape(3,1)
print("nejistoty parametru:"+str(errs2))
print("korelacni matice parametru")
corel2
c2=np.sqrt(pars2[1]**2+pars2[2]**2)
sigma2_c2 =(pars2[1]/c2)**2 * errs2[1]**2
sigma2_c2+=(pars2[2]/c2)**2 * errs2[2]**2
sigma2_c2+=pars1[1]*pars1[2]/c**2 * errs2[1]*errs2[2]*corel2[1,2]
print("parametr c: %.4f +- %.4f"%(c2,np.sqrt(sigma2_c2)))
pokusíme se o analýzu pomocí momentů: šikmosti (skew) a špičatosti (curtosis)
#mame relativne malo bodu pro histogram
xbins=np.r_[-15:15:2]
ok1=pl.hist(yval-ymodel,xbins,alpha=0.5,label="konst.sigma")
ok2=pl.hist(yval2-ymodel2,xbins,alpha=0.5,label="poisson")
pl.legend(loc=0)
def moments(meas):
cmeas=meas-meas.mean()
return (cmeas**2).mean(),(cmeas**3).mean(),(cmeas**4).mean()
moms1=moments(yval-ymodel)
skew1,curt1=moms1[1]/moms1[0]**1.5,moms1[2]/moms1[0]**2-3
skew1,curt1
moms2=moments(yval2-ymodel2)
skew2,curt2=moms1[1]/moms2[0]**1.5,moms2[2]/moms2[0]**2-3
skew2,curt2
print("nejistoty skew %.2f a curt. %.2f "%(np.sqrt(6./len(xval)),np.sqrt(24./len(xval))))
ani v jednom případě (pokud šikmosti a špičatosti přisoudíme v nulové hypotéze normální rozdělení) není odchylka o 0 větší jako 2 $\sigma$, tedy hypotézu, že původní měřená data mají normální rozdělené, nemůžeme zamítnout (ne s rizikem menším jako $2\sigma$ kvantil (4.5%).
from scipy import stats
stats.norm().sf(2)*2