In [1]:
%pylab inline
Populating the interactive namespace from numpy and matplotlib

Identifikace píku

F. James, kap. 11.5

chceme rozhodnout, zda na pozadí $b(x,\theta)$ existuje nějaká "jemná struktura" (spektrální pík, změna četnosti v čase...) v intervalu AB

Pozadí je fitováno nějakou funkcí (s vyloučením int. AB) s odhadovanými parametry $\hat\theta$ (a kovar. maticí $V$), potom očekávaná hodnota pro nulovou hypotézu je $$\hat{b}_{AB} = \int_A^B b(x,\hat\theta) dx$$ s rozptylem $\sigma^2 _{AB} = D^T V D$, kde $$D_i=\frac{\partial \hat{b}_{AB}}{\partial \theta_i}$$

Testovací statistika

$$T=\frac{(n_{AB}-\hat{b}_{AB})^2}{V(n_{AB}-\hat{b}_{AB})}$$

kde varianci můžeme uvažovat podle nulové hypotézy a Poissonovy statistiky rovnu

$$V(n_{AB})=E(n_{AB})=\hat{b}_{AB}$$

Pokud $n_{AB}$ a $\hat{b}_{AB}$ jsou nekorelovány, platí $V(n_{AB}-\hat{b}_{AB})=\hat{b}_{AB} + \sigma^2 _{AB}$.

Pro velká $n_{AB}$ má $T$ rozdělení $\chi^2(1)$ s následujícími odds (převrácená hodnota pravděpodobnosti, že struktura dosáhne významnosti alespoň daného počtu sigma)

In [2]:
from scipy import stats
d=r_[1:6]
print("odds:",zip(d,1/stats.chi(1).sf(d).astype("single")))
odds: [(1, 3.1514871), (2, 21.977894), (3, 370.39835), (4, 15787.192), (5, 1744277.9)]

Pokud nevíme přesně, kde pík (událost) hledat, je pravděpodobnost náhodného nalezení signálu v daném intervalu (obsahujícím k možných intervalů) větší než výše uvedená hodnota $p$

$$q=1-(1-p)^k \approx kp $$

Pravděpodobnost 3-sigma události, pokud máme na výběr z 50 (časových, spektrálních atp.) binů, je najednou 13%. Pokud ovšem takových událostí zaregistrujeme více, můžeme se ptát, jaká je pravděpodobnost, že (stále při platnosti nulové hypotézy) překročí např. $l$ binů z $k$ významnost s rizikem $p$:

$${k \choose l} p^l (1-p)^{k-l}$$

a pravděpod., že jich bude alespoň $l$

$$\sum_{i=l}^{k} {k \choose l} p^l (1-p)^{k-l} = 1-\sum_{i=0}^{l-1} {k \choose l} p^l (1-p)^{k-l}$$

In [5]:
#prvnich 10 clenu rady
from scipy import misc
k=20
l=r_[:10]
p=1/22.
prob20=misc.comb(k,l)*pow(p,l)*pow(1-p,k-l)
prob20
Out[5]:
array([  3.94395797e-01,   3.75615045e-01,   1.69921092e-01,
         4.85488834e-02,   9.82536925e-03,   1.49719912e-03,
         1.78237991e-04,   1.69750468e-05,   1.31354528e-06,
         8.33997006e-08])
In [6]:
#riziko !nahodneho" vyskytu alespon l binu pres 2 sigma
zip(list(l),["%.3g"%g for g in (1-cumsum(prob20))])
Out[6]:
[(0, '0.606'),
 (1, '0.23'),
 (2, '0.0601'),
 (3, '0.0115'),
 (4, '0.00169'),
 (5, '0.000197'),
 (6, '1.84e-05'),
 (7, '1.4e-06'),
 (8, '8.8e-08'),
 (9, '4.56e-09')]