Základní otázky otázka potřebnosti dalšího parametru - jaký nejmenší počet proměnných vysvětluje dostatečně data?
Matice $A$ popisuje transformaci (rotaci/inverzi) měřených veličin $$Y=A X$$
Proměnné $X$ (jednotlivé vektory) mají kovarianční matici $\Sigma$
Praktické řešení: najdeme vlastní čísla $\lambda_i$ a vlastní vektory $\pi_i$ kovar. matice, předpokládáme, že budou ortogonální (autom. splněno, pokud jsou vlastní čísla různá).
Matice $W$ vlastních vektorů matice $X^T X$ a sdružená matice $V$ vlastních vektorů matice $X X^T$ (ident. v případě čtvercové matice X) jsou transformačními maticemi singulární dekompozice matice X ve tvaru $X=W L V$, kde L je matice pouze s diagonálními nezápornými elementy.
Stopa kovar. matice je při transformaci zachována - součet variancí je součtem vlastních čísel. Vlastní vektory obvykle uspořádáme podle velikosti vlast. čísel.
[Francis] Paul J. Francis, Beverley J. Wills http://arxiv.org/abs/astro-ph/9905079 + code ref.
kombinujeme 5 funkcí: konstanta, lineární, sinusoida, gaussovka
%matplotlib inline
from matplotlib import pyplot as pl
import numpy as np
x=np.r_[0:8:100j]
fun=[np.ones_like(x),0.5*x]#0.05*x**2]
fun.append(0.3*np.sin(x*2.))
fun.append(np.exp(-(x-5)**2/2.))
fun.append(np.exp(-(x-2)**2))
nmeas=20
[pl.plot(x,f) for f in fun]
allpars=np.array([np.random.normal(a,0.4*abs(a),size=nmeas) for a in [3,-1,1.5,2,3.1]]) #sada ruznych koefientu jednotlivych komponent
testfunc=allpars.T.dot(fun)
for i in range(len(testfunc)):
testfunc[i]+=np.random.normal(0,0.3,size=len(x))
pl.plot(x,testfunc[3])
pl.plot(x,testfunc[1])
ncomp=4
renfunc4=testfunc[:ncomp]-testfunc[:ncomp].mean(1)[:,np.newaxis]
renfunc4/=np.sqrt(((renfunc4**2).sum(1)/(len(x)-1))[:,np.newaxis])
cormat4=renfunc4.dot(renfunc4.T)
eig4,vecs4=np.linalg.eig(cormat4)
eig4
Transformujeme standardizované proměnné $y=V z$
Nové proměnné jsou ortogonální, korelační matice $D=V^T\ R\ V$, kde $R=z^T\ z$ je kovarianční (korelační) matice standardizovaných prom.
Chceme maximalizovat diagonální elementy, za normalizační podmínky $v^T v=1$, pomocí Lagrangeova multiplikátoru řešíme
$$\partial (V^T\ R\ V -\lambda(v^T v-1) /\partial v=0$$ odkud vyjde rovnice pro vlast. čísla
$(R - \lambda\ \mathbb{I}) v = 0$
renfunc=testfunc-testfunc.mean(1)[:,np.newaxis]
renfunc/=np.sqrt((renfunc**2).sum(1))[:,np.newaxis]
#renfunc*=(len(x)-1)
cormat=renfunc.dot(renfunc.T)
eig,vecs=np.linalg.eig(cormat)
eig
a=renfunc[0]
(a**2).sum(),a.sum()
varlist=[v.dot(renfunc).std() for v in vecs.T]
cumlist=np.cumsum(varlist)
np.where(1-cumlist/cumlist[-1]<.10)[0][0] #kolik potrebujeme funkci k vysvetleni 90% variability
# ortonormalni vlastni vektory
np.all(np.isclose(vecs.T@vecs,np.eye(len(eig))))
qlim=9
recon=vecs.T[:,:qlim].dot(compon[:qlim])
#comnorm
compon=vecs.dot(renfunc)
comnorm=np.sqrt((compon**2).sum(1))
#compon/=comnorm[:,np.newaxis]
compon[0].dot(renfunc[0])
pl.plot(renfunc[0],renfunc[1],'.')
[pl.plot(x,v.dot(renfunc)) for v in vecs[:5]]
eig2,vecs2=np.linalg.eig(renfunc.T.dot(renfunc))
vecs2[:5].dot(renfunc[3])
# funkce z knihovny pro Machine learning
from sklearn.decomposition import PCA
zca=PCA(4)
op1=zca.fit(renfunc[:10])
zca.explained_variance_ratio_
[zca.components_.dot(r) for r in renfunc[:3]]
cormat10=renfunc[:10].dot(renfunc[:10].T)
eig10,vecs10=np.linalg.eig(cormat10)
eig10
prinfunc=vecs10.T.dot(renfunc[:10])
[pl.plot(p) for p in prinfunc[:4]]
#vazene funkce (hlavni faktory)
#prinfuncz=vecs10.T.dot(np.eye(10)/np.sqrt(eig10)).dot(renfunc[:10])
prinfuncw=((np.eye(10)/np.sqrt(eig10)).dot(prinfunc))
ix=np.r_[:100]
[pl.plot(ix,prinfuncw[:4][i]) for i in range(4)]
(prinfuncw**2).sum(1) #normalizovany rozptyl
selsrt=np.argsort(eig10)[::-1]
selsrt
[pl.plot(z+i*0.5) for i,z in enumerate(zca.components_)]
np.isclose(np.corrcoef(zca.components_),np.eye(4))
pok=zca.transform(renfunc4)
jina sada (vetsi rozptyl)
allparsb=np.array([np.random.normal(a,0.8*abs(a),size=nmeas) for a in [3,-1,1.5,2,3.1]])
testfuncb=allparsb.T.dot(fun)
for i in range(len(testfuncb)):
testfuncb[i]+=np.random.normal(0,0.3,size=len(x))
renfuncb=testfuncb-testfuncb.mean(1)[:,np.newaxis]
renfuncb/=np.sqrt((renfuncb**2).sum(1))[:,np.newaxis]
#renfunc*=(len(x)-1)
cormat10b=renfuncb[:10].dot(renfuncb[:10].T)
eig10b,vecs10b=np.linalg.eig(cormat10b)
eig10b
eig10
prinfuncb=vecs10b.T.dot(renfuncb[:10])
[pl.plot(p) for p in prinfuncb[:4]]
(eig10/eig10.sum())/(eig[:10]/eig.sum())
def regen(shift=0,fluct=0.4,noise=0.3,mabs=False,nmeas=10):
allparsb=np.array([np.random.normal(a+shift,fluct*abs(a),size=nmeas) for a in [3,-1,1.5,2,3.1]])
if mabs: allparsb=abs(allparsb)
testfuncb=allparsb.T.dot(fun)
for i in range(len(testfuncb)):
testfuncb[i]+=np.random.normal(0,noise,size=len(x))
renfuncb=testfuncb-testfuncb.mean(1)[:,np.newaxis]
renfuncb/=np.sqrt((renfuncb**2).sum(1))[:,np.newaxis]
#renfunc*=(len(x)-1)
cormat10b=renfuncb.dot(renfuncb.T)
eig10b,vecs10b=np.linalg.eig(cormat10b)
return eig10b,vecs10b.T.dot(renfuncb)
eigc,princ=regen(-2,0.8,mabs=True)