téma

jak co nejpřesněji určit polohu maxima - a s jakou nejistotou

podle teorému implicitní funkce

lze nejistotu určení maxima $\mu$ funkce $a + bx + cx^2 = q$ stanovit z diferenciálu první derivace

$$b + 2cx = 0$$$$\mathrm{d} b+ 2x \mathrm{d} c + 2c \mathrm{d} x = 0$$$$\frac{\partial x}{\partial b}= -\frac{1}{2c}$$$$\frac{\partial x}{\partial c}= -\frac{x}c = \frac{b}{2c^2}$$

s použitím explicitního řešení $x=-b/2c$. Totéž dostaneme přímým derivováním explicitní funkce.

Pro vyšší řád $$\frac{\partial^2 x}{\partial c \partial b}= \frac{1}{2c^2}$$.

Tayloruv rozvoj do 2. řádu

$$x=x_0 + \frac{\partial x}{\partial a} \Delta a + \frac{\partial x}{\partial b} \Delta b + \frac{1}2 \left (\frac{\partial^2 x}{\partial a^2} \Delta a^2 + \frac{\partial^2 x}{\partial b^2} \Delta b^2 + 2\frac{\partial^2 x}{\partial a\partial b} \Delta a \Delta b \right) + \dots$$
$$E(x)=x_0 + \frac{1}2 \left (\frac{\partial^2 x}{\partial a^2} E(\Delta a^2) + \frac{\partial^2 x}{\partial b^2} E(\Delta b^2) + 2\frac{\partial^2 x}{\partial a\partial b} E(\Delta a \Delta b) \right) + \dots$$

pokud $E(\Delta a)=E(\Delta b)=0$, ozn. kovarianci $E(\Delta a \Delta b)=\sigma_{ab}$ a $E(\Delta a^2)=\sigma_a^2$, $E(\Delta b^2)=\sigma_b^2$

$$E(x^2)=x_0^2 + x_0 \left (\frac{\partial^2 x}{\partial a^2} \sigma^2_{a} + \frac{\partial^2 x}{\partial b^2} \sigma^2_{b} + 2\frac{\partial^2 x}{\partial a\partial b} \sigma_{ab} \right) + \left(\frac{\partial x}{\partial a}\right)^2 \sigma^2_{a} + \left(\frac{\partial x}{\partial b}\right)^2 \sigma^2_{b} + 2\frac{\partial x}{\partial a} \frac{\partial x}{\partial b} \sigma_{ab} + \dots$$

nakonec (opět do 2. řádu v $\Delta a, \Delta b$) se velké závorky odečtou

$$\sigma^2_x = E(x^2)-E(x)^2 = \left(\frac{\partial x}{\partial a}\right)^2 \sigma^2_{a} + \left(\frac{\partial x}{\partial b}\right)^2 \sigma^2_{b} + 2\frac{\partial x}{\partial a} \frac{\partial x}{\partial b} \sigma_{ab} + \dots$$

Konkrétně pro fit paraboly

$$ \sigma_x^2=\frac{1}{4c^2} \sigma_b^2 + \frac{b^2}{4c^4} \sigma_c^2 + \frac{b}{2c^3} \sigma_{bc}$$

vyšší polynom

Často pík nelze dobře proložit parabolou - vykazuje jistou míru asymetrie.

Uvažujeme-li fit 4. řádu s podmínkou pro maximum $$b + 2cx + 3 ex^2 + 4 fx^3 = 0$$ a diferenciálem $$\d b+ 2x \d c + 2c \d x + 3 x^2 \d e + 6 e x \d x + 4 x^3 \d f + 12 fx^2 \d x= \d b+ 2x \d c + 3 x^2 \d e + 4 x^3 \d f + (2 c+ 6 e x + 12 fx^2) \d x= 0$$

$$\frac{\partial x}{\partial b}= -\frac{1}{m}$$$$\frac{\partial x}{\partial c}= -\frac{2x}{m}$$$$\frac{\partial x}{\partial e}= -\frac{3x^2}{m}$$$$\frac{\partial x}{\partial f}= -\frac{4x^3}{m}$$

kde $m=2 c+ 6 e x + 12 fx^2$, $x$ je třeba vyřešit (např. numericky).

In [1]:
#modelova funkce
%matplotlib inline
from matplotlib import pyplot as pl
import numpy as np
x=np.r_[:3:0.1]
mat=np.array([x**i for i in range(5)]) #polynom 4. stupne
p0=np.r_[1,7.2,-4.,0.4,0.05]
y0=p0.dot(mat)
y=y0+np.random.normal(size=x.shape[0])*0.1
pl.plot(x,y)
Out[1]:
[<matplotlib.lines.Line2D at 0x7fc8b1b44940>]
In [2]:
#nevazena linearni regrese polynomem
cov=np.linalg.inv(mat.dot(mat.T))
sol=cov.dot(mat.dot(y))
print("chi2:",((sol.dot(mat)-y)**2).sum()*100/(len(x)-5))
sol

chi2: 1.10415500218
Out[2]:
array([ 1.09156457,  6.70893116, -3.34866818,  0.07813442,  0.10271255])
In [3]:
allsol=np.roots((sol*np.r_[:5])[::-1])
sel=(allsol>0.5)*(allsol<2)
xm=allsol[sel][0] #vybereme relevantni koren v intervalu (0.5,2)
xm
Out[3]:
1.137216337552478
In [4]:
print("maxim. hodnota",y.max())
.37/(np.sqrt(2)*0.2)

maxim. hodnota 4.84636210855
Out[4]:
1.3081475451951126
In [4]:
# spocteme parcialni derivace
m=np.polyval([12*sol[4],6*sol[3],2*sol[2]],xm)
ders=np.r_[0,-1/m,-2*xm/m,-3*xm**2/m,-4*xm**3/m]
ders
Out[4]:
array([ 0.        ,  0.21932186,  0.49963058,  0.85364515,  1.29644237])
In [5]:
np.sqrt(ders.dot(cov.dot(ders)))*0.1
Out[5]:
0.015280005546266531
In [6]:
errs=np.sqrt(cov.diagonal())*0.1
errs
Out[6]:
array([ 0.07571108,  0.37485728,  0.53889655,  0.28150165,  0.04813487])

numerický experiment - 1000 náhodných sad parametrů podle hodnot sol a příslušných nejistot

In [8]:
ids=np.array([sol[i]+errs[i]*np.random.normal(size=1000) for i in range(5)]).T
allmax=np.array([np.polyval(k[::-1],xm) for k in ids])
pl.hist(allmax,20)
allmax.std()
Out[8]:
0.91911573343335662
In [14]:
#poloha maxima urcena z derivace 
allpos=[np.roots((k*np.r_[:5])[::-1])[2] for k in ids]
In [15]:
np.array(allpos).real.mean()
Out[15]:
1.1901803626699805
In [16]:
np.array(allpos).real.std()
Out[16]:
0.31102091875588916
In [18]:
ook=pl.hist(np.array(allpos).real,20)
In [15]:
(cov/errs[:,np.newaxis]/errs[np.newaxis,:]+0.5).astype(int)
Out[15]:
array([[100, -81,  67, -57,  52],
       [-81, 100, -95,  91, -84],
       [ 67, -95, 100, -97,  96],
       [-57,  91, -97, 100, -98],
       [ 52, -84,  96, -98, 100]])
In [78]:
cval,cvec=np.linalg.eig(cov)
rvals=cvec.dot(np.sqrt(cval)[:,np.newaxis]*0.1*np.random.normal(size=(len(sol),3000)))
rvals.std(1)
Out[78]:
array([ 0.07618767,  0.37310925,  0.53460714,  0.2785802 ,  0.04754585])
In [79]:
dvals=rvals-rvals.mean(1).reshape(5,1)
dvals.dot(dvals.T)/30./cov
Out[79]:
array([[ 1.01262941,  1.00075793,  0.99690938,  0.99416659,  0.99149703],
       [ 1.00075793,  0.99069537,  0.98721037,  0.98426587,  0.98157077],
       [ 0.99690938,  0.98721037,  0.98414411,  0.98155859,  0.97921345],
       [ 0.99416659,  0.98426587,  0.98155859,  0.97935147,  0.97735764],
       [ 0.99149703,  0.98157077,  0.97921345,  0.97735764,  0.97567592]])
In [86]:
k=rvals[:,2]
np.roots(((k+sol)*np.r_[:5])[:0:-1]),k
Out[86]:
(array([-9.31374046,  3.1218292 ,  1.12744772]),
 array([-0.00314229,  0.20909618, -0.42477534,  0.25725955, -0.04767332]))
In [95]:
almax=np.array([np.roots(((sol+k)*np.r_[:5])[:0:-1])[2] for k in rvals.T[:2000]])
almax.mean(),almax.std()
Out[95]:
(1.135317759198921, 0.01510632918891267)
In [96]:
ook=pl.hist(almax,20)