simulujeme N=2000 opakování normálně rozdělených dat, vykreslíme ve scatter plotu i pro kernel density estimate (KDE)
In [1]:
%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,'.')
Out[1]:
In [2]:
#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])
Out[2]:
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
In [14]:
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
In [4]:
kernel.covariance
Out[4]:
In [5]:
# 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())
Out[5]:
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) $$
In [6]:
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()
Out[6]:
In [7]:
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
In [8]:
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")
Out[8]:
In [9]:
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")
Out[9]:
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})$$
In [26]:
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()
Out[26]:
In [25]:
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())
Out[25]:
směrodatná odchylka z odhadovaného $r$ ze vzorku
$$\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
In [45]:
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)$)
In [41]:
from scipy import stats
tlim=stats.t(n1-2).isf(0.05)
tlim,sum(t>tlim)/200.
Out[41]:
In [22]:
tlim2=stats.t(n2-2).isf(0.05)
tlim2,sum(t2>tlim2)/200.
Out[22]:
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.
In [42]:
ms=measure(20)
plot(ms[0],ms[1],'o')
Out[42]:
In [43]:
array(ms)
Out[43]:
In [44]:
rho=corrcoef(array(ms))[1,0]
print rho,rho*sqrt(18/(1-rho**2))
array(ms)[:,randint(0,20,20)]
Out[44]:
In [39]:
bset=[corrcoef(array(ms)[:,randint(0,20,20)])[1,0] for i in range(10000)]
In [43]:
hist(bset,20)
mean(bset),sum(r_[bset]<=0)
Out[43]:
In [48]:
tset=r_[bset]*sqrt(18/(1-r_[bset]**2))
hist(tset,20)
tset.mean(),tset.std()
Out[48]:
In [52]:
# srovnavaci test
ms2=measure(20,0.)
bset2=[corrcoef(array(ms2)[:,randint(0,20,20)])[1,0] for i in range(10000)]
In [55]:
hist(bset2,20)
rho2=corrcoef(array(ms2))[1,0]
print rho2,rho2*sqrt(18/(1-rho2**2))
In [56]:
tset2=r_[bset2]*sqrt(18/(1-r_[bset2]**2))
hist(tset2,20)
tset2.mean(),tset2.std()
xlabel("t-hodnota")
Out[56]:
předpokládáme stejnou nejistotu $\sigma$ v obou osách (jinak lze jednu osu přeškálovat - viz níže)
reziduum pro lineární fit je $$h_i=\frac{y-m x_i - c}{\sqrt{1+m^2}}$$
derivací výrazu $\sum{h_i^2}$ dostáváme pro m řešení
$$\hat m=A\pm\sqrt{A^2+1}$$kde $$A=\frac{V(y)-V(x)}{2 V(x,y)}$$
In [23]:
rho,s1,s2=(r1*r2).mean(),(r1**2).mean(),(r2**2).mean()
plot(n1,n2,'.m')
plot(r_[-2:2:.1],rho/s1*r_[-2:2:.1])
plot(r_[-1:1:.1]*rho/s2,r_[-1:1:.1])
A=(s2-s1)/2./rho
if rho>0: m=A+sqrt(A**2+1)
else: m=A-sqrt(A**2+1)
#plot(r_[-3:3:.1]/m,r_[-3:3:.1],'k')
plot(r_[-3:3:.1],m*r_[-3:3:.1],'k')
print(A,m,rho/s2,rho/s1)
legend(["regr.Y=f(X)","regr.X=f(Y)","orthogon."])
#xlim(-6,6)
grid()
osa mezi regresními přímkami $\tan (\phi_1+\phi_2)/2=f(t_1,t_2)$, kde $t_1=\tan \phi_1=\rho/s_1$, $t_2=\tan \phi_2=s_2/\rho$ $$t=\tan(\phi_1+\phi_2)=\frac{t_1+t_2}{1-t_1\ t_2}$$ $$f(t_1,t_2)=\frac{\sqrt{t^2+1}-1}{t}=\frac{\sqrt{(t_1+t_2)^2+(1-t_1\ t_2)^2}-1+t_1\ t_2}{t_1+t_2}=\frac{\sqrt{(1+t_1^2)(1+t_2^2)}-1+t_1\ t_2}{t_1+t_2}$$
funguje jen pro lineární fit - pro zakřivené modely je ortogonální vzdálenost podstatně komplikovanější!