In [2]:
%pylab inline
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
In [4]:
from scipy import stats
d=r_[1:6]
print "odds:",zip(d,1/stats.chi(1).sf(d).astype("single"))
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 [7]:
#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[7]:
In [12]:
#riziko !nahodneho" vyskytu alespon l binu pres 2 sigma
["%.3g"%g for g in (1-cumsum(prob20))]
Out[12]: