Testování hypotéz

Z předchozího měření (nebo teorie) činíme hypotézu o hodnotě neznámého parametru $\theta_0$ (odhad $\hat\theta_0 \pm \delta$ s pravděp. obsahem p) , ev. o rozdělení měřených hodnot. Následující měření umožní test hypotézy H0:

měření dává statistiku $t(y_1,y_2,..,y_N)$, úkolem je stanovit, jaká je pravděpodobnost pozorování $t$ za předpokladu platnosti/neplatnosti H0

určujeme tedy, s jakym rizikem nastane jedna ze 2 možných chyb

  • chyba 1. druhu = hypotéza platí, ale H0 zamítneme: $P(t \in K\ |\ H_0)=\alpha$
  • chyba 2. druhu = hypotéza neplatí (platí alternativa H1, pokud je jediná), ale H0 přijmeme: $P(t \notin K\ |\ \mathrm{not}\ H_0)=\beta$

kde K je kritická oblast "nepřijatelných" hodnot $t$ (pokud $t$ zde leží, můžeme H0 zamítnout na hladině významnosti $\alpha$ - jinak je naše statistika nedostatečná a vyžaduje další měření). Pravděpodobnost $1-\beta$ je pak síla (mohutnost) testu (závisí na $\alpha$)...

In [4]:
%matplotlib inline
M,N=12,10 #velikosti vzorků
mu=2.
sig=1
from matplotlib import pyplot as pl
from scipy import stats
import numpy as np
x=np.r_[-3:8:0.01]
hyp0=stats.t(M+N-2,loc=0)
hyp1=stats.t(M+N-2,loc=mu)
y1=hyp0.pdf(x)
y2=hyp1.pdf(x)
pl.plot(x,y1)
pl.plot(x,y2)
prob=0.05
pos=hyp0.ppf(1-prob) # hodnoty mensi nez tato 
pl.axvline(pos,color='r')
pl.fill([pos]+list(x[x>pos]),[0]+list(y1[x>pos]),'b',alpha=0.4,hold=1)
pl.fill(list(x[x<pos])+[pos],list(y2[x<pos])+[0],'g',alpha=0.4,hold=1)
pl.legend(["H0 - nul.","H1","prob %.2f"%prob,"alpha","beta"])
Out[4]:
<matplotlib.legend.Legend at 0x7f41fe31c750>

Neyman-Pearsonův test

Volba kritické oblasti podle poměru $f_N(X|\theta_1)/f_N(X|\theta_0) \gt c_\alpha$, hodnota $c_\alpha$ zvolena podle dané hladiny významnosti; v principu jde o poměr věrohodností.

Test dobré shody

určení stupně shody získaných dat s hypotézou o jejich rozdělení

Pearsonův test

z naměřených hodnot sestavíme histogram o K buňkách, do každé z nich padne $n_i$ z měřených hodnot ($\sum n_i = N$). Hypotéza H0 stanoví distrib. funkci rozdělení F, tedy pro každou buňku pravděpodobnost $p_i=P[x \in \lt l_i,m_i)]=F(m_i)-F(l_i)$, kde $l_i, m_i$ jsou hranice i-té buňky (obvykle $l_i=m_{i-1}$). Náhodné proměnné $n_i$ mají binomické rozdělení a očekávanou hodnotu $E(n_i)=N p_i$.

Pearsonova statistika

$$T=\sum_i^k \frac{(n_i-N p_i)^2}{N p_i} = \frac{1}{N} \sum_i^k \left( \frac{n_i^2}{p_i} - 2 n_iN + N^2 p_i \right) = \frac{1}{N} \sum_i^k \frac{n_i^2}{p_i} - N$$

v limitě $N \to \infty$ (kdy rozdělení $n_i$ se blíží normálnímu) má $T$ rozdělení $\chi_{K-1}^2$ (Pearsonova věta).

Podmínkou je malé množství buňek, kde není dostatečné množství měřených hodnot (lze řešit úpravou mezí buňek / slučováním tak, aby $p_i$ byly vyrovnané).

Rozdělení může být závislé na parametrech $\theta_1, ..., \theta_r$ - pak hledáme hodnotu $\mathbf{\theta}$ dávající nejlepší shodu, testovaná statistika má rozdělení s méně stupni volnosti, někde mezi $\chi_{K-1}^2$ a $\chi_{K-r-1}^2$ (rozdíl je malý, pokud $K \gg r$).

In [55]:
N=1000
dx=0.5
ndf=5
dset=stats.chi2(ndf).rvs(N)
x=np.r_[0:15:dx]
ok=pl.hist(dset,x)
pred=stats.chi2(ndf).pdf(x+0.25)
pl.plot(x+0.25,pred*N*dx,'r',lw=2)
pl.title("Chi^2 n.d.f.=5 [1000 samp.]")
Out[55]:
<matplotlib.text.Text at 0x7f41fde5f210>
In [28]:
resid=ok[0]-pred[:-1]*N*dx
pl.plot(x[:-1]+dx/2.,resid,'d')
pl.axhline(0)
pl.grid()
pl.title("rezidua")
print("red. suma rezidui:",sum(resid**2./(pred[:-1]*N*dx)))
print("kvantily Chi2 CDF:",stats.chi2(len(ok[0])-1).isf([0.1,0.05,0.01]))

('red. suma rezidui:', 16.881612184700749)
('kvantily Chi2 CDF:', array([ 37.91592254,  41.33713815,  48.27823577]))
In [20]:
def testN(ndf,dat,dx=0.5):
    pred=stats.chi2(ndf).pdf(x[:len(dat)]+dx/2.)
    return sum((dat-pred*sum(dat)*dx)**2./(pred*sum(dat)*dx))
testN(4,ok[0]),testN(6,ok[0])
Out[20]:
(155.49945566272714, 149.73151462190449)

nulovou hypotézu (teor. křivka a histogram se neliší) nemůžeme zamítnout ani na hladině spolehlivosti 90%.

Kolgomorův test

sestavení distribuční funkce z N naměřených dat $S_N(X)$ = (počet měření<x) / N.

statistika $D_N=\max \|S_N(x)-F(x)\|$ pro velká N (>80) má rozdělení

$$F_{Klg}(\sqrt{N}D_N)=1-2 \sum_j^\infty (-1)^{j-1} \exp(-2j^2 N D_N^2)$$

H0 zamítáme s rizikem $\alpha$, pokud vyšlo $\sqrt{N}D_N > z_\alpha$

$\alpha$ $z_\alpha$
(Kolgomorov-
Smirnovo rozdělení)
0.01 1.63
0.05 1.36
0.10 1.22
In [66]:
edf=np.cumsum(ok[0])
edf/=float(N)
pred2=stats.chi2(ndf).cdf(x[1:])

dsetx=np.sort(dset)
pl.plot(dsetx,np.r_[:1:1j*N],'g.',label="EDF")
pl.bar(x[1:]-dx/2.,edf,width=dx,label="EDF z histogramu",alpha=0.5)
#pl.plot(x[1:]+0.25,edf,label="EDF")
pl.plot(x[1:],pred2,'r',label="Chi2 cumul.DF")
pl.xlim(0,15)
pl.grid()
Dn=max(abs(edf-pred2))
pl.legend(loc=4)
print("Komogorov D, Z alpha")
Dn,Dn*np.sqrt(N)

Komogorov D, Z alpha
Out[66]:
(0.033970166397133017, 1.074230983098563)

Opakujeme K-S test 1000krát - vyneseme rozdělení parametru $D \sqrt{N}$

In [23]:
N2=1000
multest=np.array([max(abs(np.cumsum(np.histogram(stats.chi2(ndf).rvs(N),x)[0])/float(N)-pred2)) for i in range(N2)])
multest*=np.sqrt(N)
pl.hist(multest,30)
pl.axvline(1.22,color='r')
pl.axvline(1.36,color='r')
print("mean %.4f st.d. %.4f "%(multest.mean(),multest.std()))

mean 0.7243 st.d. 0.2455
In [27]:
print('podil mereni nad kvantily 1.22, 1.36')
sum(multest>1.22)/float(N2),sum(multest>1.36)/float(N2)

podil mereni nad kvantily 1.22, 1.36
Out[27]:
(0.042000000000000003, 0.014)

Test homogenity

rovnost rozdělení 2 vzorků dat

Smirnovův test

srovnáváme empirické dist. funkce $S_n$, $S_m$: statistika $D_{nm}=\max \|S_n(x)-S_m(x)\|$ má asymptotické rozdělení jako v předchozím případě se $\sqrt{nm/(n+m)}D_{nm} < z_\alpha$

nebo jednostranné testy $D_{nm}^{\pm}=\max \{\pm (S_n(x)-S_m(x))\}$ má jednodušší rozdělení $K^\pm(z)=\lim P(\sqrt{nm/(n+m)}D_{nm}^{\pm} \lt z_\alpha)=1-\exp(-2 z^2)$ pro $n\rightarrow\infty, m\rightarrow\infty$

úpravou vzniká

Wayneův test

$$W_{nm}=\sqrt{nm/(n+m)}\ \max \{ \|S_n(x)-S_m(x)\|/S_m(x) \}$$

pro $x$ splňující $S(x)\geq a$ (kde $a>0$ je konstanta) má rozdělení

$$\lim P(W_{nm}^{\pm} \lt z_\alpha)=\sqrt{\frac{2}{\pi}} \int_0^{z\sqrt{a/(1-a)}} \exp\left(\frac{t^2}{2}\right) dt$$

pro $z \gt 0$.

Wilcoxonův (-Mann-Whitney) test

2 vzorky spojíme, (n+m) čísel uspořádáme podle velikosti a určíme pořadí (pozice) čísel $x_1, .., x_n$ z první sady: jejich součet je $T_n$ a výraz $$U_0=\frac{\frac{nm + n(n+1)}{2}-T_n}{\sqrt{\frac{nm}{12}(n+m+1)}}$$ má rozdělení N(0,1).

rovnost středních hodnot

předpokládáme normální rozdělení N($\mu_x,\sigma_x^2$), N($\mu_y,\sigma_y^2$)

hypotéza H0: $\mu_x=\mu_y$

Studentův test

zde neznáme disperzi, předpokládáme ale, že jsou $\sigma_x$ a $\sigma_y$ identické

testovací statistika $$t=\frac{\bar{x}-\bar{y}}{\sqrt{\frac{\left((n-1)\widehat{\sigma_x^2} + (m-1)\widehat{\sigma_y^2}\right)(n+m)}{nm(n+m-2)}}}$$

kde $\bar{x}=\sum_i{x_i}/n$, $\widehat{\sigma_x^2}=\sum_i{(x_i-\bar{x})^2}/(n-1)$ (resp. pro y a m) má Studentovo rozdělení s n+m-2 stupni volnosti.

Pozn.: Disperze rozdílu $\Delta = \bar{x} -\bar{y}$ je odhadnuta jako $$\widehat{\sigma_{\Delta}^2} = \widehat{\sigma_{\bar{x}}^2} + \widehat{\sigma_{\bar{y}}^2} = \frac{n+m}{nm} \widehat{\sigma^2},$$ kde $$\widehat{\sigma^2} = \frac{(n-1)\widehat{\sigma_x^2} + (m-1)\widehat{\sigma_y^2}}{n+m-2}$$ je společný odhad kombinující odhady disperze vzorků $\widehat{\sigma_x^2}$ a $\widehat{\sigma_y^2}$.

rovnost disperzí

Fischerův test

hypotéza H0: $\sigma_x=\sigma_y$

statistika $F=\widehat{\sigma_x^2}/\widehat{\sigma_y^2}$ (definice viz výše) má Fischer-Snedecorovo rozdělení $F_{n-1,m-1}$

In [42]:
from numpy import random
N,M=50,80
samp1=np.random.normal(3,1.,size=N)
samp2=np.random.normal(3,1.4,size=M)
sig1=sum((samp1-samp1.mean())**2)/(N-1)
sig2=sum((samp2-samp2.mean())**2)/(M-1)
sig1,sig2
Out[42]:
(0.83454589573730831, 2.0767641595179112)
In [53]:
print('"experiment": %.5f resp. %.5f'%(sig2/sig1,sig1/sig2))
print("kvantily F resp. 1/F (role vzorku prohozeny):")
x=np.r_[0.01:2.5:0.05]
pl.plot(x,stats.f(M,N).pdf(x),label="[%i,%i]"%(M,N))
pl.plot(x,stats.f(M,N).pdf(1/x),label="[%i,%u] 1/x"%(M,N))
pl.grid()
pl.title("Fisher. rozdel.")
pl.legend()
alphas=[0.1,0.05,0.025,0.01]
kvans=stats.f(N-1,M-1).isf(alphas)
for i in range(4):
    pl.axvline(kvans[i],color="b",ls=':')
for i in range(4):
    pl.axvline(1/kvans[i],color="g",ls=':')
stats.f(N-1,M-1).isf(alphas),1/stats.f(M-1,N-1).ppf(alphas)

"experiment": 2.48850 resp. 0.40185
kvantily F resp. 1/F (role vzorku prohozeny):
Out[53]:
(array([ 1.38026223,  1.51310274,  1.63831575,  1.79674087]),
 array([ 1.38026223,  1.51310274,  1.63831575,  1.79674087]))