výpočet těžiště

váhy $w_i$ odpovídají hodnotám binů histogramu (nad jistý práh) nebo logaritmu hodnot

$x_{mean} = \sum x_i w_i / \sum w_i$

nejistota

  • $w_i=a N_i$

    • $\delta w_i=a \sqrt N_i$
    • $\delta \sum x_i w_i = \sqrt{\sum x_i^2 (\delta w_i)^2} = a \sqrt{\sum x_i^2 N_i} = a \sqrt{\bar {x^2}}\ \sqrt{\sum N_i}$
    • $$ \frac{\partial x_{mean}}{\partial w_i}= \frac{x_i}{\sum_j w_j} - \frac{\sum_j w_j x_j}{(\sum_j w_j)^2} = \frac{x_i-x_{mean}}{a \sum_j N_j}$$ tedy (váhy jsou nekorelované) $$\delta x_{mean}^2 = \frac{1}{(\sum_j N_j)^2} \sum N_i (x_i-x_{mean})^2 $$

  • $w_i=\log(a N_i)$

    • $\delta w_i=\frac{a \sqrt N_i}{a N_i} = 1/\sqrt N_i$

  • $w_i=N_i-l$ ($l$ - volba hladiny)

v našem případě je šum konstatní (vnější zdroj)

In [96]:
%pylab inline

from scipy import stats

dset1=stats.poisson(30).rvs(1000)
stids='151629,151633,151637,151639,151643,151644,151654,151656,151659'.split(',')
stids+='151660,151662,151666,151667,151668,144824,151669,151675,151677,151678,151679,151680,151683,151688'.split(',')
stids=open("../Prakt/slist2016.txt").read().split(', ')
len(stids)

Populating the interactive namespace from numpy and matplotlib

C:\Program Files\Anaconda3\lib\site-packages\IPython\core\magics\pylab.py:161: UserWarning: pylab import has clobbered these variables: ['norm']
`%matplotlib` prevents importing * from pylab and numpy
  "\n`%matplotlib` prevents importing * from pylab and numpy"
Out[96]:
23

demonstrace

simulovan asymetricky pik

In [2]:
x=r_[-10.:10:100j]
cln=exp(-(x-5)**2/4-3e-2*(x-5)**3-3e-3*(x-5)**4)
plot(x,cln)
Out[2]:
[<matplotlib.lines.Line2D at 0x2c864f956d8>]
In [100]:
x=r_[0:10:100j]
#cln=exp(-(x-5)**2/4-.3e-1*(x-5)**3-.1e-4*(x-5)**4)
cln=exp(-(x-5)**2/4-2e-2*(x-5)**3-3e-3*(x-5)**4)
nois=0.1
rln=cln+stats.norm(0,nois).rvs(size=cln.shape)
plot(x,rln)
Out[100]:
[<matplotlib.lines.Line2D at 0x2c869a3f400>]
In [4]:
a,b=rln,x

levs=r_[.1:.9:15j]
pos=[((a[a>l]-l)*b[a>l]).sum()/(a[a>l]-l).sum() for l in levs] # spoctene polohy teziste
plot(levs,pos)
xlabel("hladina")
ylabel("poloha teziste")
Out[4]:
<matplotlib.text.Text at 0x2c8628549b0>
In [5]:
plot(b,a)
plot(pos,levs,'r')
axhline(0.4)
sel=a>0.4
fill_between(b[sel],0.4,a[sel],alpha=0.2)
xlabel('x')
ylabel('y')
text(5,0.2,'poloha teziste',color='r')
Out[5]:
<matplotlib.text.Text at 0x2c862770d68>

možno nafitovat polynomem (2. stupně) nebo exponenciální závislostí (nelineární regrese) - zde dostaneme jasnou asymptotickou hodnotu $x_{mean}$

In [59]:
errorbar(levs,pos,0.1*mb/norm**2)
p=pres[0]
plot(levs,p[0]+p[1]*exp(-levs/p[2]))
p=pres2[0]
plot(levs,p[0]+p[1]*exp(-levs/p[2]))
Out[59]:
[<matplotlib.lines.Line2D at 0xaf7c76ec>]
In [60]:
pres[0][0],pres2[0][0]
Out[60]:
(4.9214886056260507, 4.9575933950220721)