po převedení na standardiz. proměnné $X_m=(X-E(X))/\sqrt{V(X)}$
dostáváme korelační koef. $r=X_m Y_m$
udává i sklon regresní přímky (bez standardizace je to poměr $cov(X,Y)/V(X)$). Nejistota $V(\hat r)=V(Y)/(N V(X))$
simulujeme N=2000 opakování normálně rozdělených dat, vykreslíme ve scatter plotu i pro kernel density estimate (KDE)
%pylab inline
import numpy as np
from scipy import stats
def measure(n,frac=0.5):
"Measurement model, return two coupled measurements."
l1 = np.random.normal(size=n)
l2 = np.random.normal(scale=0.5, size=n)
return l1+frac*l2, frac*l1-l2
m1, m2 = measure(2000)
xmin = m1.min()
xmax = m1.max()
ymin = m2.min()
ymax = m2.max()
plot(m1,m2,'.')
#Perform a kernel density estimate on the data:
X, Y = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]
positions = np.vstack([X.ravel(), Y.ravel()])
values = np.vstack([m1, m2])
kernel = stats.gaussian_kde(values)
kdeZ = np.reshape(kernel(positions).T, X.shape)
imshow(kdeZ,origin="lower",extent=[xmin, xmax, ymin, ymax])
stats.gaussian_kde?
porovnání korelační matice a odhadu z KDE
teoretický rozptyl
$V(m_1)=V(l_1)+f^2 V(l_2)$ = 1+0.25 * 0.25=1.063
$V(m_2)=f^2 V(l_1) + V(l_2)$ = 0.25 + 0.25=0.5
print corrcoef(m1,m2)
r1=m1-m1.mean()
r2=m2-m2.mean()
print (r1*r2).sum()/sqrt((r1**2).sum()*(r2**2).sum()),
print "zname stred:",(m1*m2).sum()/sqrt((m1**2).sum()*(m2**2).sum()),
print "zname sigma:",(m1*m2).sum()/sqrt(0.5*1.0625)/2000
kernel.covariance
# mensi stupen korelace
n1,n2=measure(2000,0.1)
r1=n1-n1.mean()
r2=n2-n2.mean()
(r1*r2).sum()/sqrt((r1**2).sum()*(r1**2).sum())
Jaký bude výsledek při opakování pokusu?
teoretické rozdělení (Fisher 1944) $$prob(r|\rho,H) =\frac{(1-\rho^2)^{(N-1)/2} (1-r^2)^{(N-4)/2} }{(1-\rho r)^{N-3/2}} \left( 1+ \frac{1}{N-1/2} \frac{1+\rho r}{8} + \dots \right) $$
poloha maxima z 1. derivace $$-(N/2-2)\frac{(1-\rho^2)^{(N-1)/2} 2 r (1-r^2)^{N/2-3} }{(1-\rho r)^{N-3/2}}+\rho(N-3/2)\frac{(1-\rho^2)^{(N-1)/2} (1-r^2)^{N/2-2} }{(1-\rho r)^{N-1/2}}=0$$
tefun=lambda r,rho,N:pow(1-rho**2,(N-1)/2.)*pow(1-r**2,(N-4)/2.)/pow(1-rho*r,N-1.5)*(1+(1+rho*r)/8./(N-0.5))
n1,n2=10,20
cores=r_[[corrcoef(*measure(n1,0.5))[1,0] for i in range(200)]]
h1=hist(cores,20,alpha=0.5)
cores2=r_[[corrcoef(*measure(n2,0.5))[1,0] for i in range(200)]]
h2=hist(cores2,h1[1],alpha=0.5)
rd=np.r_[-0.5:1:0.05]
v1=tefun(rd,0.5,n1)
plot(rd,v1/v1.sum()*200)
v2=tefun(rd,0.5,n2)
plot(rd,v2/v2.sum()*200)
v1.max(),v2.max()
from scipy.integrate import quad
qlist=[.25,.5,.75]
cnorm=[[quad(lambda x:tefun(x,q,i),-0.999,0.999)[0] for i in range(2,30)] for q in qlist]
cmids=[[quad(lambda x:tefun(x,q,i),0,0.999)[0] for i in range(2,30)] for q in qlist]
podle teoret. rozdělení určíme pravděpodobnost, že hodnota $r$ překročí jistou mez (0, 0.5) pro danou hodnotu $\rho$
cquant=array(cmids)/array(cnorm)
[plot(range(2,30),c) for c in cquant]
legend(qlist,loc=4)
xlabel("velikost vzorku")
ylabel("prst. nezapornosti r")
clev1=[[quad(lambda x:tefun(x,q,i),0.5,0.999)[0] for i in range(2,30)] for q in qlist]
[plot(range(2,30),c) for c in array(clev1)/array(cnorm)]
legend(qlist,loc=4)
xlabel("velikost vzorku")
ylabel("prst. r>0.5")
podle Wolfram Mathworld
$$E(r)=\rho - \frac{\rho(1-\rho^2)}{2n}$$$$V(r)=\frac{(1-\rho^2)^2}{n}(1+\frac{11 \rho^2}{2n})$$n1,n2=10,20
rm1,rm2=0.5-0.5*(1-0.25)/(2*n1),0.5-0.5*(1-0.25)/(2*n2)
print "střední hodnota - teorie: ",rm1,rm2
print "simulace:"
sum(rd*v1)/v1.sum(),sum(rd*v2)/v2.sum()
print "směrod. odchylka - teorie: ",sqrt((1-0.25)**2/n1*(1+11*0.25/2./n1)),sqrt((1-0.25)**2/n2*(1+11*0.25/2./n2))
print "simulace:"
sqrt(sum((rd-rm1)**2*v1)/v1.sum()),sqrt(sum((rd-rm2)**2*v2)/v2.sum())
Máme-li jen jedno měření, lze odhadnout směrodatná odchylka určovaného $r$ ze vzorku podle vztahu $$\sigma_r=\frac{1-r^2}{\sqrt{n-1}}$$
Při testování nulové hypotézy ($\rho=0$) lze použít transformaci na proměnnou $$t=r \frac{\sqrt{n-2}}{\sqrt{1-r^2}}$$ která má Studentovo rozdělení s $n-2$ stupni volnosti (původní proměnná $r$ má pro $\rho=0$ rozdělení $prob(r)= (1-r^2)^{N/2-2}$).
n1,n2=10,20
t=cores/sqrt(1-cores**2)*sqrt(n1-2)
t2=cores2/sqrt(1-cores2**2)*sqrt(n2-2)
qa,qb,cx=hist(t,20,alpha=0.5)
qa2,qb,cx2=hist(t2,qb,alpha=0.5)
print '10 samp.:',t.mean(),t.std()
print '20 samp.:',t2.mean(),t2.std()
Vynesené histogramy samozřejmě neodpovídají Studentovu rozdělení, které má střední hodnotu 0 (mohou být srovnána s necentrálním Studentovým rozdělením, které má složitější tvar).
provádíme Studentův test (isf = inverse survival function dává kvantil $P(\alpha)$)
from scipy import stats
tlim=stats.t(n1-2).isf(0.05)
tlim,sum(t>tlim)/200.
tlim2=stats.t(n2-2).isf(0.05)
tlim2,sum(t2>tlim2)/200.
V prvním případě (korelace určené z 10 bodů) můžeme 'nulovou hypotézu' vyloučit v polovině případů, při určování korelací z 20 bodů už v 88 procentech.
bootstrap = náhodné výběry s opakováním ze stále téhož vzorku (simulace opakování měření)
ms=measure(20)
plot(ms[0],ms[1],'o')
array(ms)
rho=corrcoef(array(ms))[1,0]
print rho,rho*sqrt(18/(1-rho**2))
array(ms)[:,randint(0,20,20)]
bset=[corrcoef(array(ms)[:,randint(0,20,20)])[1,0] for i in range(10000)]
hist(bset,20)
mean(bset),sum(r_[bset]<=0)
převádí blíže normálnímu rozdělení
tset=r_[bset]*sqrt(18/(1-r_[bset]**2))
hist(tset,20)
tset.mean(),tset.std()
# srovnavaci test
ms2=measure(20,0.)
bset2=[corrcoef(array(ms2)[:,randint(0,20,20)])[1,0] for i in range(10000)]
hist(bset2,20)
rho2=corrcoef(array(ms2))[1,0]
print rho2,rho2*sqrt(18/(1-rho2**2))
tset2=r_[bset2]*sqrt(18/(1-r_[bset2]**2))
hist(tset2,20)
tset2.mean(),tset2.std()
xlabel("t-hodnota")
alternativně Spearman rank test
$$ r_s= 1-\frac{6 \sum d_i^2}{n(n^2-1)}$$
bližší konceptu maximální korelace