Pozorování pan Hubblea

snažíme se nalézt parametry fitu

$$\log(v) = b + a m$$

kdy $m$ je magnituda dané galaxie (úměrná záporně vzatému logaritmu zářivého výkonu).

In [3]:
#vtup
import numpy as np
dd='''Virgo 890 7 12:5
Pegasus 3810 5 15:5
Pisces 4630 4 15:4
Cancer 4820 2 16:0
Perseus 5230 4 16:4
Coma 7500 3 17:0
Ursa Major 11;800 1 18:0
Leo 19;600 1 19:0
(No name) 2350 16 13:8
(No name) 630 21 11:6'''.replace(':','.').replace(';','').split('\n')
data=np.array([np.array(d.split()[-3:]).astype(float) for d in dd])
v,N,m=data.T#.astype(float)
In [4]:
lv=np.log10(v)
#lv.var(),((lv-lv.mean())**2).mean()
sigma_n=1.6
w=1/(0.1*lv/N)**2/sigma_n**2 #vahy
znorm=w.sum()
zmean=lambda z:sum(z*w)/znorm
zmean(lv),lv.mean()
Out[4]:
(3.004443978629423, 3.600705245422423)
In [5]:
sum(lv*w),sum(m*w),znorm
Out[5]:
(10532.319335671613, 43697.293040294855, 3505.58020405369)
In [6]:
covar=(lv*m).mean()-lv.mean()*m.mean()
b=covar/m.var()
a=(lv-b*m).mean()
b,a
Out[6]:
(0.19666889208742763, 0.548404040225546)
In [7]:
covar2=zmean(lv*m)-zmean(lv)*zmean(m)
mvar=zmean(m**2)-zmean(m)**2
b2=covar2/mvar
a2=zmean(lv-b2*m)
b2,a2
Out[7]:
(0.22789515450957298, 0.163715548695727)
In [8]:
from matplotlib import pyplot as pl
%matplotlib inline
pl.errorbar(m,lv,np.sqrt(1/w),fmt='s')
x=np.r_[11:20:1]
pl.plot(x,a+b*x)
pl.plot(x,a2+b2*x)
pl.grid()
In [9]:
dder=[[1,sum(w*m)/znorm],[sum(w*m)/znorm,sum(w*m**2)/znorm]]
np.linalg.det(dder)
Out[9]:
1.6555449506902173
In [10]:
cmat=np.linalg.inv(dder)
np.sqrt(cmat.diagonal())
Out[10]:
array([9.73925178, 0.77719413])
In [11]:
sigma2=zmean(lv-a2+b2*m)
np.sqrt(sigma2/(len(lv)-2))
Out[11]:
0.8427230312999782
In [12]:
cmat,dder
Out[12]:
(array([[94.85302514, -7.52928349],
        [-7.52928349,  0.60403071]]),
 [[1, 12.465067263263677], [12.465067263263677, 157.03344682837803]])
In [13]:
sum(w)/znorm
Out[13]:
1.0
In [65]:
wmat=np.eye(len(lv))*w/znorm#len(lv)
amat=np.array([np.ones(len(lv)),m])
cmat@amat@wmat@lv
Out[65]:
array([0.16371555, 0.22789515])
In [70]:
cmat=cov
errs=np.sqrt(cmat.diagonal())
cmat[1,0]/errs[0]/errs[1]
Out[70]:
-0.9947147197430005
In [66]:
b2,a2
Out[66]:
(0.2278951545095726, 0.16371554869573177)

maticový přístup

In [15]:
import numpy as np
amat=np.array([m,np.ones(len(m))]).T
wmat=np.eye(len(m))/(0.1*lv/N)**2/sigma_n**2
hess=amat.T@wmat@amat
hess
Out[15]:
array([[550493.34257588,  43697.29304029],
       [ 43697.29304029,   3505.58020405]])
In [16]:
np.linalg.det(hess)
Out[16]:
20345145.14794945
In [17]:
cov=np.linalg.inv(hess)
np.sqrt(cov.diagonal())
Out[17]:
array([0.01312652, 0.16449233])
In [18]:
pars=cov@amat.T@wmat@lv
pars #indenticke jako vyse
Out[18]:
array([0.22789515, 0.16371555])
In [19]:
model=amat@pars
res=lv-model
S0=res@wmat@res
S0
Out[19]:
7.94228360453256
In [20]:
len(m)
Out[20]:
10
In [21]:
sigma_n2=np.sqrt(S0/(len(m)-2))
sigma_n2
Out[21]:
0.9963861954917731
In [22]:
errs=np.sqrt(cov.diagonal())
cov/errs[:,np.newaxis]/errs[np.newaxis,:]
Out[22]:
array([[ 1.        , -0.99471472],
       [-0.99471472,  1.        ]])

modifikace modelu

přidáme kvadratický člen (před ostatní)

In [23]:
amat2=np.array([m**2,m,np.ones(len(m))]).T
hess2=amat2.T@wmat@amat2
cov2=np.linalg.inv(hess2)
In [25]:
errs2=np.sqrt(cov2.diagonal())
errs2,errs
Out[25]:
(array([0.00837316, 0.22415133, 1.48354789]), array([0.01312652, 0.16449233]))

Nejistoty s dalším parametrem výrazně narostou!
Důvodem jsou značné korelace mezi parametry.

In [26]:
cov2/errs2[:,np.newaxis]/errs2[np.newaxis,:]
Out[26]:
array([[ 1.        , -0.99828383,  0.99383406],
       [-0.99828383,  1.        , -0.99858726],
       [ 0.99383406, -0.99858726,  1.        ]])
In [27]:
pars2=cov2 @ amat2.T @ wmat @ lv
pars2
Out[27]:
array([-0.01794391,  0.70743338, -2.99596601])

Studentův t-test

porovnání (absolutní) hodnoty parametru a směrodatné odchylky - podíl má t-rozdělení s počtem stup. volnosti $n-k$ ($k$ je počet parametrů)

In [30]:
from scipy import stats as st
tval=abs(pars2)/errs2
print(tval)
print("pravdepodobnost t>tval [v %]")
st.t(8).sf(tval)*100

[2.14302811 3.15605253 2.01946026]
pravdepodobnost t>tval [v %]
Out[30]:
array([3.22358166, 0.6737459 , 3.90619814])
In [31]:
# mez (kvantil) pro 5% riziko chyby 1. druhu 
st.t(8).ppf(0.95)
Out[31]:
1.8595480375228424
In [32]:
pl.errorbar(m,lv,np.sqrt(1/w),fmt='s')
x=np.r_[11:20:1]
#a,b,c=pars2
pl.plot(x,np.polyval(pars2,x))
#pl.plot(x,a*x)
a,b=pars
pl.plot(x,a*x+b)
Out[32]:
[<matplotlib.lines.Line2D at 0x190560bdc88>]
In [118]:
res=lv-amat2@pars2
S2=res@wmat@res
S2/(len(lv)-4)
Out[118]:
0.5543906478246166

exponenciala

In [104]:
cof=np.r_[0.1:2:0.2]
chi2=[]
for c in cof:
    amat3=np.array([np.exp(-m*c),np.ones(len(m))]).T
    hess3=amat3.T@wmat@amat3
    cov3=np.linalg.inv(hess3)
    pars3=cov3@amat3.T@wmat@lv
    print(pars3)
    res=lv-amat3@pars3
    S3=res@wmat@res
    chi2.append(S3/(len(lv)-2))

[-8.6340463   5.50623196]
[-38.69802676   3.98340174]
[-294.6757235     3.68048197]
[-2575.45926447     3.554196  ]
[-2.39028696e+04  3.48741171e+00]
[-2.29039873e+05  3.44759150e+00]
[-2.23631111e+06  3.42200006e+00]
[-2.20934308e+07  3.40464577e+00]
[-2.19959117e+08  3.39237535e+00]
[-2.20140204e+09  3.38340075e+00]
In [106]:
np.array(chi2)*8
Out[106]:
array([ 4.11514555,  4.7742524 , 10.46227904, 17.28580147, 23.74982234,
       29.43287414, 34.30694734, 38.45851768, 41.99089845, 44.9966265 ])
In [97]:
res=lv-amat3@pars3
S3=res@wmat@res
S3/(len(lv)-2)
Out[97]:
3.3369231054617794
In [100]:
from scipy import stats as st
st.chi2(len(lv)-2).ppf(0.95)
Out[100]:
15.50731305586545
In [98]:
S3
Out[98]:
26.695384843694235