876 lines
24 KiB
Python
876 lines
24 KiB
Python
"""Module contain algorithms for low-rank models.
|
|
|
|
There is almost no typechecking of any kind here, just focus on speed
|
|
"""
|
|
|
|
import math
|
|
import warnings
|
|
from scipy.linalg import svd,inv
|
|
from scipy import dot,empty,eye,newaxis,zeros,sqrt,diag,\
|
|
apply_along_axis,mean,ones,randn,empty_like,outer,r_,c_,\
|
|
rand,sum,cumsum,matrix, expand_dims,minimum,where,arange,inner,tile
|
|
has_sym = True
|
|
has_arpack = True
|
|
try:
|
|
from symeig import symeig
|
|
except:
|
|
has_sym = False
|
|
try:
|
|
from scipy.sandbox import arpack
|
|
except:
|
|
has_arpack = False
|
|
|
|
|
|
def pca(a, aopt,scale='scores',mode='normal',center_axis=0):
|
|
""" Principal Component Analysis.
|
|
|
|
Performs PCA on given matrix and returns results in a dictionary.
|
|
|
|
:Parameters:
|
|
a : array
|
|
Data measurement matrix, (samples x variables)
|
|
aopt : int
|
|
Number of components to use, aopt<=min(samples, variables)
|
|
|
|
:Returns:
|
|
results : dict
|
|
keys -- values, T -- scores, P -- loadings, E -- residuals,
|
|
lev --leverages, ssq -- sum of squares, expvar -- cumulative
|
|
explained variance, aopt -- number of components used
|
|
|
|
:OtherParam eters:
|
|
mode : str
|
|
Amount of info retained, ('fast', 'normal', 'detailed')
|
|
center_axis : int
|
|
Center along given axis. If neg.: no centering (-inf,..., matrix modes)
|
|
|
|
:SeeAlso:
|
|
- pcr : other blm
|
|
- pls : other blm
|
|
- lpls : other blm
|
|
|
|
Notes
|
|
-----
|
|
Uses kernel speed-up if m>>n or m<<n.
|
|
|
|
If residuals turn rank deficient, a lower number of component than given
|
|
in input will be used. The number of components used is given in
|
|
results-dict.
|
|
|
|
Examples
|
|
--------
|
|
|
|
>>> import scipy,engines
|
|
>>> a=scipy.asarray([[1,2,3],[2,4,5]])
|
|
>>> dat=engines.pca(a, 2)
|
|
>>> dat['expvarx']
|
|
array([0.,99.8561562, 100.])
|
|
|
|
"""
|
|
m, n = a.shape
|
|
assert(aopt<=min(m,n))
|
|
if center_axis>=0:
|
|
a = a - expand_dims(a.mean(center_axis), center_axis)
|
|
if m>(n+100) or n>(m+100):
|
|
u, s, v = esvd(a, amax=None) # fixme:amax option need to work with expl.var
|
|
else:
|
|
u, s, vt = svd(a, 0)
|
|
v = vt.T
|
|
e = s**2
|
|
tol = 1e-10
|
|
eff_rank = sum(s>s[0]*tol)
|
|
aopt = minimum(aopt, eff_rank)
|
|
T = u*s
|
|
s = s[:aopt]
|
|
T = T[:,:aopt]
|
|
P = v[:,:aopt]
|
|
|
|
if scale=='loads':
|
|
T = T/s
|
|
P = P*s
|
|
|
|
if mode == 'fast':
|
|
return {'T':T, 'P':P, 'aopt':aopt}
|
|
|
|
if mode=='detailed':
|
|
E = empty((aopt, m, n))
|
|
ssq = []
|
|
lev = []
|
|
for ai in range(aopt):
|
|
E[ai,:,:] = a - dot(T[:,:ai+1], P[:,:ai+1].T)
|
|
ssq.append([(E[ai,:,:]**2).mean(0), (E[ai,:,:]**2).mean(1)])
|
|
if scale=='loads':
|
|
lev.append([((s*T)**2).sum(1), (P**2).sum(1)])
|
|
else:
|
|
lev.append([(T**2).sum(1), ((s*P)**2).sum(1)])
|
|
else:
|
|
# residuals
|
|
E = a - dot(T, P.T)
|
|
#E = a
|
|
SEP = E**2
|
|
ssq = [SEP.sum(0), SEP.sum(1)]
|
|
# leverages
|
|
if scale=='loads':
|
|
lev = [(1./m)+(T**2).sum(1), (1./n)+((P/s)**2).sum(1)]
|
|
else:
|
|
lev = [(1./m)+((T/s)**2).sum(1), (1./n)+(P**2).sum(1)]
|
|
# variances
|
|
expvarx = r_[0, 100*e.cumsum()/e.sum()][:aopt+1]
|
|
|
|
return {'T':T, 'P':P, 'E':E, 'expvarx':expvarx, 'levx':lev, 'ssqx':ssq, 'aopt':aopt, 'eigvals': e[:aopt,newaxis]}
|
|
|
|
def pcr(a, b, aopt, scale='scores',mode='normal',center_axis=0):
|
|
""" Principal Component Regression.
|
|
|
|
Performs PCR on given matrix and returns results in a dictionary.
|
|
|
|
:Parameters:
|
|
a : array
|
|
Data measurement matrix, (samples x variables)
|
|
b : array
|
|
Data response matrix, (samples x responses)
|
|
aopt : int
|
|
Number of components to use, aopt<=min(samples, variables)
|
|
|
|
:Returns:
|
|
results : dict
|
|
keys -- values, T -- scores, P -- loadings, E -- residuals,
|
|
levx -- leverages, ssqx -- sum of squares, expvarx -- cumulative
|
|
explained variance, aopt -- number of components used
|
|
|
|
:OtherParameters:
|
|
mode : str
|
|
Amount of info retained, ('fast', 'normal', 'detailed')
|
|
center_axis : int
|
|
Center along given axis. If neg.: no centering (-inf,..., matrix modes)
|
|
|
|
:SeeAlso:
|
|
- pca : other blm
|
|
- pls : other blm
|
|
- lpls : other blm
|
|
|
|
Notes
|
|
-----
|
|
|
|
Uses kernel speed-up if m>>n or m<<n.
|
|
|
|
If residuals turn rank deficient, a lower number of component than given
|
|
in input will be used. The number of components used is given in results-dict.
|
|
|
|
|
|
Examples
|
|
--------
|
|
|
|
>>> import scipy,engines
|
|
>>> a=scipy.asarray([[1,2,3],[2,4,5]])
|
|
>>> b=scipy.asarray([[1,1],[2,3]])
|
|
>>> dat=engines.pcr(a, 2)
|
|
>>> dat['expvarx']
|
|
array([0.,99.8561562, 100.])
|
|
|
|
"""
|
|
k, l = m_shape(b)
|
|
if center_axis>=0:
|
|
b = b - expand_dims(b.mean(center_axis), center_axis)
|
|
dat = pca(a, aopt=aopt, scale=scale, mode=mode, center_axis=center_axis)
|
|
T = dat['T']
|
|
weights = apply_along_axis(vnorm, 0, T)**2
|
|
if scale=='loads':
|
|
Q = dot(b.T, T*weights)
|
|
else:
|
|
Q = dot(b.T, T/weights)
|
|
|
|
if mode=='fast':
|
|
dat.update({'Q':Q})
|
|
return dat
|
|
if mode=='detailed':
|
|
F = empty((aopt, k, l))
|
|
for i in range(aopt):
|
|
F[i,:,:] = b - dot(T[:,:i+1], Q[:,:i+1].T)
|
|
else:
|
|
F = b - dot(T, Q.T)
|
|
expvary = r_[0, 100*((T**2).sum(0)*(Q**2).sum(0)/(b**2).sum()).cumsum()[:aopt]]
|
|
#fixme: Y-var leverages
|
|
dat.update({'Q':Q, 'F':F, 'expvary':expvary})
|
|
return dat
|
|
|
|
def pls(a, b, aopt=2, scale='scores', mode='normal', center_axis=-1, ab=None):
|
|
"""Partial Least Squares Regression.
|
|
|
|
Performs PLS on given matrix and returns results in a dictionary.
|
|
|
|
:Parameters:
|
|
a : array
|
|
Data measurement matrix, (samples x variables)
|
|
b : array
|
|
Data response matrix, (samples x responses)
|
|
aopt : int
|
|
Number of components to use, aopt<=min(samples, variables)
|
|
|
|
:Returns:
|
|
results : dict
|
|
keys -- values, T -- scores, P -- loadings, E -- residuals,
|
|
levx -- leverages, ssqx -- sum of squares, expvarx -- cumulative
|
|
explained variance of descriptors, expvary -- cumulative explained
|
|
variance of responses, aopt -- number of components used
|
|
|
|
:OtherParameters:
|
|
mode : str
|
|
Amount of info retained, ('fast', 'normal', 'detailed')
|
|
center_axis : int
|
|
Center along given axis. If neg.: no centering (-inf,..., matrix modes)
|
|
|
|
:SeeAlso:
|
|
- pca : other blm
|
|
- pcr : other blm
|
|
- lpls : other blm
|
|
|
|
Notes
|
|
-----
|
|
|
|
Uses kernel speed-up if m>>n or m<<n.
|
|
|
|
If residuals turn rank deficient, a lower number of component than given
|
|
in input will be used. The number of components used is given in results-dict.
|
|
|
|
Examples
|
|
--------
|
|
|
|
>>> import scipy,engines
|
|
>>> a=scipy.asarray([[1,2,3],[2,4,5]])
|
|
>>> b=scipy.asarray([[1,1],[2,3]])
|
|
>>> dat=engines.pls(a, b, 2)
|
|
>>> dat['expvarx']
|
|
array([0.,99.8561562, 100.])
|
|
|
|
"""
|
|
|
|
m, n = m_shape(a)
|
|
if ab!=None:
|
|
mm, l = m_shape(ab)
|
|
assert(m==mm)
|
|
else:
|
|
k, l = m_shape(b)
|
|
|
|
if center_axis>=0:
|
|
a = a - expand_dims(a.mean(center_axis), center_axis)
|
|
b = b - expand_dims(b.mean(center_axis), center_axis)
|
|
|
|
W = empty((n, aopt))
|
|
P = empty((n, aopt))
|
|
R = empty((n, aopt))
|
|
Q = empty((l, aopt))
|
|
T = empty((m, aopt))
|
|
B = empty((aopt, n, l))
|
|
tt = empty((aopt,))
|
|
|
|
if ab==None:
|
|
ab = dot(a.T, b)
|
|
for i in range(aopt):
|
|
if ab.shape[1]==1: #pls 1
|
|
w = ab.reshape(n, l)
|
|
w = w/vnorm(w)
|
|
elif n<l: # more yvars than xvars
|
|
if has_sym:
|
|
s, w = symeig(dot(ab, ab.T),range=[n,n],overwrite=True)
|
|
else:
|
|
w, s, vh = svd(dot(ab, ab.T))
|
|
w = w[:,:1]
|
|
else: # standard wide xdata
|
|
if has_sym:
|
|
s, q = symeig(dot(ab.T, ab),range=[l,l],overwrite=True)
|
|
else:
|
|
q, s, vh = svd(dot(ab.T, ab))
|
|
q = q[:,:1]
|
|
w = dot(ab, q)
|
|
w = w/vnorm(w)
|
|
r = w.copy()
|
|
if i>0:
|
|
for j in range(0, i, 1):
|
|
r = r - dot(P[:,j].T, w)*R[:,j][:,newaxis]
|
|
|
|
t = dot(a, r)
|
|
tt[i] = tti = dot(t.T, t).ravel()
|
|
p = dot(a.T, t)/tti
|
|
q = dot(r.T, ab).T/tti
|
|
ab = ab - dot(p, q.T)*tti
|
|
T[:,i] = t.ravel()
|
|
W[:,i] = w.ravel()
|
|
|
|
if mode=='fast' and i==aopt-1:
|
|
if scale=='loads':
|
|
tnorm = sqrt(tt)
|
|
T = T/tnorm
|
|
W = W*tnorm
|
|
return {'T':T, 'W':W}
|
|
|
|
P[:,i] = p.ravel()
|
|
R[:,i] = r.ravel()
|
|
Q[:,i] = q.ravel()
|
|
#B[i] = dot(R[:,:i+1], Q[:,:i+1].T)
|
|
|
|
|
|
|
|
qnorm = apply_along_axis(vnorm, 0, Q)
|
|
tnorm = sqrt(tt)
|
|
pp = (P**2).sum(0)
|
|
if mode=='detailed':
|
|
E = empty((aopt, m, n))
|
|
F = empty((aopt, k, l))
|
|
ssqx, ssqy = [], []
|
|
leverage = empty((aopt, m))
|
|
h2x = [] #hotellings T^2
|
|
h2y = []
|
|
for ai in range(aopt):
|
|
E[ai,:,:] = a - dot(T[:,:ai+1], P[:,:ai+1].T)
|
|
F[i-1] = b - dot(T[:,:i], Q[:,:i].T)
|
|
ssqx.append([(E[ai,:,:]**2).mean(0), (E[ai,:,:]**2).mean(1)])
|
|
ssqy.append([(F[ai,:,:]**2).mean(0), (F[ai,:,:]**2).mean(1)])
|
|
leverage[ai,:] = 1./m + ((T[:,:ai+1]/tnorm[:ai+1])**2).sum(1)
|
|
h2y.append(1./k + ((Q[:,:ai+1]/qnorm[:ai+1])**2).sum(1))
|
|
else:
|
|
# residuals
|
|
E = a - dot(T, P.T)
|
|
F = b - dot(T, Q.T)
|
|
sepx = E**2
|
|
ssqx = [sepx.sum(0), sepx.sum(1)]
|
|
sepy = F**2
|
|
ssqy = [sepy.sum(0), sepy.sum(1)]
|
|
# leverage
|
|
leverage = 1./m + ((T/tnorm)**2).sum(1)
|
|
h2x = []
|
|
h2y = []
|
|
# variances
|
|
tp= tt*pp
|
|
tq = tt*qnorm*qnorm
|
|
expvarx = r_[0, 100*tp/(a*a).sum()]
|
|
expvary = r_[0, 100*tq/(b*b).sum()]
|
|
|
|
if scale=='loads':
|
|
T = T/tnorm
|
|
W = W*tnorm
|
|
Q = Q*tnorm
|
|
P = P*tnorm
|
|
|
|
return {'Q':Q, 'P':P, 'T':T, 'W':W, 'R':R, 'E':E, 'F':F,
|
|
'expvarx':expvarx, 'expvary':expvary, 'ssqx':ssqx, 'ssqy':ssqy,
|
|
'leverage':leverage, 'h2':h2x}
|
|
|
|
def w_simpls(aat, b, aopt):
|
|
""" Simpls for wide matrices.
|
|
Fast pls for crossval, used in calc rmsep for wide X
|
|
There is no P or W. T is normalised
|
|
"""
|
|
bb = b.copy()
|
|
m, m = aat.shape
|
|
U = empty((m, aopt)) # W
|
|
T = empty((m, aopt))
|
|
H = empty((m, aopt)) # R
|
|
PROJ = empty((m, aopt)) # P?
|
|
|
|
for i in range(aopt):
|
|
q, s, vh = svd(dot(dot(b.T, aat), b), full_matrices=0)
|
|
u = dot(b, q[:,:1]) #y-factor scores
|
|
U[:,i] = u.ravel()
|
|
t = dot(aat, u)
|
|
t = t/vnorm(t)
|
|
T[:,i] = t.ravel()
|
|
h = dot(aat, t) #score-weights
|
|
H[:,i] = h.ravel()
|
|
PROJ[:,:i+1] = dot(T[:,:i+1], inv(dot(T[:,:i+1].T, H[:,:i+1])) )
|
|
if i<aopt:
|
|
b = b - dot(PROJ[:,:i+1], dot(H[:,:i+1].T,b) )
|
|
C = dot(bb.T, T)
|
|
|
|
return {'T':T, 'U':U, 'Q':C, 'H':H}
|
|
|
|
def w_pls(aat, b, aopt):
|
|
""" Pls for wide matrices.
|
|
Fast pls for crossval, used in calc rmsep for wide X
|
|
There is no P or W. T is normalised
|
|
|
|
aat = centered kernel matrix
|
|
b = centered y
|
|
"""
|
|
bb = b.copy()
|
|
k, l = m_shape(b)
|
|
m, m = m_shape(aat)
|
|
U = empty((m, aopt)) # W
|
|
T = empty((m, aopt))
|
|
R = empty((m, aopt)) # R
|
|
PROJ = empty((m, aopt)) # P?
|
|
|
|
for i in range(aopt):
|
|
if has_sym:
|
|
s, q = symeig(dot(dot(b.T, aat), b), range=(l,l),overwrite=True)
|
|
else:
|
|
q, s, vh = svd(dot(dot(b.T, aat), b), full_matrices=0)
|
|
q = q[:,:1]
|
|
u = dot(b , q) #y-factor scores
|
|
U[:,i] = u.ravel()
|
|
t = dot(aat, u)
|
|
|
|
t = t/vnorm(t)
|
|
T[:,i] = t.ravel()
|
|
r = dot(aat, t)#score-weights
|
|
#r = r/vnorm(r)
|
|
R[:,i] = r.ravel()
|
|
PROJ[:,: i+1] = dot(T[:,:i+1], inv(dot(T[:,:i+1].T, R[:,:i+1])) )
|
|
if i<aopt:
|
|
b = b - dot(PROJ[:,:i+1], dot(R[:,:i+1].T, b) )
|
|
C = dot(bb.T, T)
|
|
|
|
return {'T':T, 'U':U, 'Q':C, 'R':R}
|
|
|
|
def bridge(a, b, aopt, scale='scores', mode='normal', r=0):
|
|
"""Undeflated Ridged svd(X'Y)
|
|
"""
|
|
m, n = m_shape(a)
|
|
k, l = m_shape(b)
|
|
u, s, vt = svd(b, full_matrices=0)
|
|
g0 = dot(u*s, u.T)
|
|
g = (1 - r)*g0 + r*eye(m)
|
|
ag = dot(a.T, g)
|
|
u, s, vt = svd(ag, full_matrices=0)
|
|
W = u[:,:aopt]
|
|
K = vt[:aopt,:].T
|
|
T = dot(a, W)
|
|
tnorm = apply_along_axis(vnorm, 0, T) # norm of T-columns
|
|
|
|
if mode == 'fast':
|
|
if scale=='loads':
|
|
T = T/tnorm
|
|
W = W*tnorm
|
|
return {'T':T, 'W':W}
|
|
|
|
U = dot(g0, K) #fixme check this
|
|
Q = dot(b.T, dot(T, inv(dot(T.T, T)) ))
|
|
B = zeros((aopt, n, l), dtype='f')
|
|
for i in range(aopt):
|
|
B[i] = dot(W[:,:i+1], Q[:,:i+1].T)
|
|
|
|
if mode == 'detailed':
|
|
E = empty((aopt, m, n))
|
|
F = empty((aopt, k, l))
|
|
for i in range(aopt):
|
|
E[i] = a - dot(T[:,:i+1], W[:,:i+1].T)
|
|
F[i] = b - dot(a, B[i])
|
|
else: #normal
|
|
F = b - dot(a, B[-1])
|
|
E = a - dot(T, W.T)
|
|
|
|
if scale=='loads':
|
|
T = T/tnorm
|
|
W = W*tnorm
|
|
Q = Q*tnorm
|
|
|
|
return {'B':B, 'W':W, 'T':T, 'Q':Q, 'E':E, 'F':F, 'U':U, 'P':W}
|
|
|
|
|
|
def nipals_lpls(X, Y, Z, a_max, alpha=.7, mean_ctr=[2, 0, 1], scale='scores', verbose=False):
|
|
""" L-shaped Partial Least Sqaures Regression by the nipals algorithm.
|
|
|
|
(X!Z)->Y
|
|
:input:
|
|
X : data matrix (m, n)
|
|
Y : data matrix (m, l)
|
|
Z : data matrix (n, o)
|
|
|
|
:output:
|
|
T : X-scores
|
|
W : X-weights/Z-weights
|
|
P : X-loadings
|
|
Q : Y-loadings
|
|
U : X-Y relation
|
|
L : Z-scores
|
|
K : Z-loads
|
|
B : Regression coefficients X->Y
|
|
b0: Regression coefficient intercept
|
|
evx : X-explained variance
|
|
evy : Y-explained variance
|
|
evz : Z-explained variance
|
|
mnx : X location
|
|
mny : Y location
|
|
mnz : Z location
|
|
|
|
:Notes:
|
|
|
|
"""
|
|
if mean_ctr!=None:
|
|
xctr, yctr, zctr = mean_ctr
|
|
X, mnX = center(X, xctr)
|
|
Y, mnY = center(Y, yctr)
|
|
Z, mnZ = center(Z, zctr)
|
|
|
|
varX = (X**2).sum()
|
|
varY = (Y**2).sum()
|
|
varZ = (Z**2).sum()
|
|
|
|
m, n = X.shape
|
|
k, l = Y.shape
|
|
u, o = Z.shape
|
|
|
|
# initialize
|
|
U = empty((k, a_max))
|
|
Q = empty((l, a_max))
|
|
T = empty((m, a_max))
|
|
W = empty((n, a_max))
|
|
P = empty((n, a_max))
|
|
K = empty((o, a_max))
|
|
L = empty((u, a_max))
|
|
B = empty((a_max, n, l))
|
|
#b0 = empty((a_max, 1, l))
|
|
var_x = empty((a_max,))
|
|
var_y = empty((a_max,))
|
|
var_z = empty((a_max,))
|
|
|
|
MAX_ITER = 250
|
|
LIM = 1e-1
|
|
for a in range(a_max):
|
|
if verbose:
|
|
print "\nWorking on comp. %s" %a
|
|
u = Y[:,:1]
|
|
diff = 1
|
|
niter = 0
|
|
while (diff>LIM and niter<MAX_ITER):
|
|
niter += 1
|
|
u1 = u.copy()
|
|
w = dot(X.T, u)
|
|
w = w/sqrt(dot(w.T, w))
|
|
#w = w/dot(w.T, w)
|
|
l = dot(Z, w)
|
|
k = dot(Z.T, l)
|
|
k = k/sqrt(dot(k.T, k))
|
|
#k = k/dot(k.T, k)
|
|
w = alpha*k + (1-alpha)*w
|
|
#print sqrt(dot(w.T, w))
|
|
w = w/sqrt(dot(w.T, w))
|
|
t = dot(X, w)
|
|
c = dot(Y.T, t)
|
|
c = c/sqrt(dot(c.T, c))
|
|
u = dot(Y, c)
|
|
diff = dot((u-u1).T, (u-u1))
|
|
if verbose:
|
|
print "Converged after %s iterations" %niter
|
|
print "Error: %.2E" %diff
|
|
tt = dot(t.T, t)
|
|
p = dot(X.T, t)/tt
|
|
q = dot(Y.T, t)/tt
|
|
l = dot(Z, w)
|
|
|
|
U[:,a] = u.ravel()
|
|
W[:,a] = w.ravel()
|
|
P[:,a] = p.ravel()
|
|
T[:,a] = t.ravel()
|
|
Q[:,a] = q.ravel()
|
|
L[:,a] = l.ravel()
|
|
K[:,a] = k.ravel()
|
|
|
|
X = X - dot(t, p.T)
|
|
Y = Y - dot(t, q.T)
|
|
Z = (Z.T - dot(w, l.T)).T
|
|
|
|
var_x[a] = pow(X, 2).sum()
|
|
var_y[a] = pow(Y, 2).sum()
|
|
var_z[a] = pow(Z, 2).sum()
|
|
|
|
B[a] = dot(dot(W[:,:a+1], inv(dot(P[:,:a+1].T, W[:,:a+1]))), Q[:,:a+1].T)
|
|
#b0[a] = mnY - dot(mnX, B[a])
|
|
|
|
|
|
# variance explained
|
|
evx = 100.0*(1 - var_x/varX)
|
|
evy = 100.0*(1 - var_y/varY)
|
|
evz = 100.0*(1 - var_z/varZ)
|
|
if scale=='loads':
|
|
tnorm = apply_along_axis(vnorm, 0, T)
|
|
T = T/tnorm
|
|
W = W*tnorm
|
|
Q = Q*tnorm
|
|
knorm = apply_along_axis(vnorm, 0, K)
|
|
L = L*knorm
|
|
K = K/knorm
|
|
|
|
return {'T':T, 'W':W, 'P':P, 'Q':Q, 'U':U, 'L':L, 'K':K, 'B':B, 'evx':evx, 'evy':evy, 'evz':evz,'mnx': mnX, 'mny': mnY, 'mnz': mnZ}
|
|
|
|
|
|
|
|
def nipals_pls(X, Y, a_max, alpha=.7, ax_center=0, mode='normal', scale='scores', verbose=False):
|
|
"""Partial Least Sqaures Regression by the nipals algorithm.
|
|
|
|
(X!Z)->Y
|
|
:input:
|
|
X : data matrix (m, n)
|
|
Y : data matrix (m, l)
|
|
|
|
:output:
|
|
T : X-scores
|
|
W : X-weights/Z-weights
|
|
P : X-loadings
|
|
Q : Y-loadings
|
|
U : X-Y relation
|
|
B : Regression coefficients X->Y
|
|
b0: Regression coefficient intercept
|
|
evx : X-explained variance
|
|
evy : Y-explained variance
|
|
evz : Z-explained variance
|
|
|
|
:Notes:
|
|
|
|
"""
|
|
if ax_center>=0:
|
|
mn_x = expand_dims(X.mean(ax_center), ax_center)
|
|
mn_y = expand_dims(Y.mean(ax_center), ax_center)
|
|
X = X - mn_x
|
|
Y = Y - mn_y
|
|
|
|
varX = pow(X, 2).sum()
|
|
varY = pow(Y, 2).sum()
|
|
|
|
m, n = X.shape
|
|
k, l = Y.shape
|
|
|
|
# initialize
|
|
U = empty((k, a_max))
|
|
Q = empty((l, a_max))
|
|
T = empty((m, a_max))
|
|
W = empty((n, a_max))
|
|
P = empty((n, a_max))
|
|
B = empty((a_max, n, l))
|
|
b0 = empty((a_max, m, l))
|
|
var_x = empty((a_max,))
|
|
var_y = empty((a_max,))
|
|
|
|
t1 = X[:,:1]
|
|
for a in range(a_max):
|
|
if verbose:
|
|
print "\n Working on comp. %s" %a
|
|
u = Y[:,:1]
|
|
diff = 1
|
|
MAX_ITER = 100
|
|
lim = 1e-16
|
|
niter = 0
|
|
while (diff>lim and niter<MAX_ITER):
|
|
niter += 1
|
|
#u1 = u.copy()
|
|
w = dot(X.T, u)
|
|
w = w/sqrt(dot(w.T, w))
|
|
#l = dot(Z, w)
|
|
#k = dot(Z.T, l)
|
|
#k = k/sqrt(dot(k.T, k))
|
|
#w = alpha*k + (1-alpha)*w
|
|
#w = w/sqrt(dot(w.T, w))
|
|
t = dot(X, w)
|
|
q = dot(Y.T, t)
|
|
q = q/sqrt(dot(q.T, q))
|
|
u = dot(Y, q)
|
|
diff = vnorm(t1 - t)
|
|
t1 = t.copy()
|
|
if verbose:
|
|
print "Converged after %s iterations" %niter
|
|
#tt = dot(t.T, t)
|
|
#p = dot(X.T, t)/tt
|
|
#q = dot(Y.T, t)/tt
|
|
#l = dot(Z, w)
|
|
p = dot(X.T, t)/dot(t.T, t)
|
|
p_norm = vnorm(p)
|
|
t = t*p_norm
|
|
w = w*p_norm
|
|
p = p/p_norm
|
|
|
|
U[:,a] = u.ravel()
|
|
W[:,a] = w.ravel()
|
|
P[:,a] = p.ravel()
|
|
T[:,a] = t.ravel()
|
|
Q[:,a] = q.ravel()
|
|
|
|
X = X - dot(t, p.T)
|
|
Y = Y - dot(t, q.T)
|
|
|
|
var_x[a] = pow(X, 2).sum()
|
|
var_y[a] = pow(Y, 2).sum()
|
|
|
|
B[a] = dot(dot(W[:,:a+1], inv(dot(P[:,:a+1].T, W[:,:a+1]))), Q[:,:a+1].T)
|
|
b0[a] = mn_y - dot(mn_x, B[a])
|
|
|
|
# variance explained
|
|
evx = 100.0*(1 - var_x/varX)
|
|
evy = 100.0*(1 - var_y/varY)
|
|
|
|
if scale=='loads':
|
|
tnorm = apply_along_axis(vnorm, 0, T)
|
|
T = T/tnorm
|
|
W = W*tnorm
|
|
Q = Q*tnorm
|
|
|
|
return {'T':T, 'W':W, 'P':P, 'Q':Q, 'U':U, 'B':B, 'b0':b0, 'evx':evx, 'evy':evy,
|
|
'mnx': mnX, 'mny': mnY, 'xc': X, 'yc': Y}
|
|
|
|
|
|
########### Helper routines #########
|
|
|
|
def m_shape(array):
|
|
return matrix(array).shape
|
|
|
|
def esvd(data, amax=None):
|
|
"""SVD with the option of economy sized calculation
|
|
Calculate subspaces of X'X or XX' depending on the shape
|
|
of the matrix.
|
|
|
|
Good for extreme fat or thin matrices
|
|
|
|
:notes:
|
|
Numpy supports this by setting full_matrices=0
|
|
"""
|
|
m, n = data.shape
|
|
if m>=n:
|
|
kernel = dot(data.T, data)
|
|
if has_arpack:
|
|
if amax==None:
|
|
amax = n
|
|
s, v = arpack.eigen_symmetric(kernel,k=amax, which='LM',
|
|
maxiter=200,tol=1e-5)
|
|
if has_sym:
|
|
if amax==None:
|
|
amax = n
|
|
pcrange = None
|
|
else:
|
|
pcrange = [n-amax, n]
|
|
s, v = symeig(kernel, range=pcrange, overwrite=True)
|
|
s = s[::-1].real
|
|
v = v[:,::-1].real
|
|
else:
|
|
u, s, vt = svd(kernel)
|
|
v = vt.T
|
|
s = sqrt(s)
|
|
u = dot(data, v)/s
|
|
else:
|
|
kernel = dot(data, data.T)
|
|
if has_sym:
|
|
if amax==None:
|
|
amax = m
|
|
pcrange = None
|
|
else:
|
|
pcrange = [m-amax, m]
|
|
s, u = symeig(kernel, range=pcrange, overwrite=True)
|
|
s = s[::-1]
|
|
u = u[:,::-1]
|
|
else:
|
|
u, s, vt = svd(kernel)
|
|
s = sqrt(s)
|
|
v = dot(data.T, u)/s
|
|
# some use of symeig returns the 0 imaginary part
|
|
return u.real, s.real, v.real
|
|
|
|
def vnorm(x):
|
|
# assume column arrays (or vectors)
|
|
return math.sqrt(dot(x.T, x))
|
|
|
|
def center(a, axis):
|
|
# 0 = col center, 1 = row center, 2 = double center
|
|
# -1 = nothing
|
|
|
|
# check if we have a vector
|
|
is_vec = len(a.shape)==1
|
|
if not is_vec:
|
|
is_vec = a.shape[0]==1 or a.shape[1]==1
|
|
if is_vec:
|
|
if axis==2:
|
|
warnings.warn("Double centering of vecor ignored, using ordinary centering")
|
|
if axis==-1:
|
|
mn = 0
|
|
else:
|
|
mn = a.mean()
|
|
return a - mn, mn
|
|
# !!!fixme: use broadcasting
|
|
if axis==-1:
|
|
mn = zeros((1,a.shape[1],))
|
|
#mn = tile(mn, (a.shape[0], 1))
|
|
elif axis==0:
|
|
mn = a.mean(0)[newaxis]
|
|
#mn = tile(mn, (a.shape[0], 1))
|
|
elif axis==1:
|
|
mn = a.mean(1)[:,newaxis]
|
|
#mn = tile(mn, (1, a.shape[1]))
|
|
elif axis==2:
|
|
mn = a.mean(0)[newaxis] + a.mean(1)[:,newaxis] - a.mean()
|
|
return a - mn , a.mean(0)[newaxis]
|
|
else:
|
|
raise IOError("input error: axis must be in [-1,0,1,2]")
|
|
|
|
return a - mn, mn
|
|
|
|
def scale(a, axis):
|
|
if axis==-1:
|
|
sc = zeros((a.shape[1],))
|
|
elif axis==0:
|
|
sc = a.std(0)
|
|
elif axis==1:
|
|
sc = a.std(1)[:,newaxis]
|
|
else:
|
|
raise IOError("input error: axis must be in [-1,0,1]")
|
|
|
|
return a - sc, sc
|
|
|
|
|
|
|
|
## #PCA CALCS
|
|
|
|
## % Calculate Q limit using unused eigenvalues
|
|
## temp = diag(s);
|
|
## if n < m
|
|
## emod = temp(lv+1:n,:);
|
|
## else
|
|
## emod = temp(lv+1:m,:);
|
|
## end
|
|
## th1 = sum(emod);
|
|
## th2 = sum(emod.^2);
|
|
## th3 = sum(emod.^3);
|
|
## h0 = 1 - ((2*th1*th3)/(3*th2^2));
|
|
## if h0 <= 0.0
|
|
## h0 = .0001;
|
|
## disp(' ')
|
|
## disp('Warning: Distribution of unused eigenvalues indicates that')
|
|
## disp(' you should probably retain more PCs in the model.')
|
|
## end
|
|
## q = th1*(((1.65*sqrt(2*th2*h0^2)/th1) + 1 + th2*h0*(h0-1)/th1^2)^(1/h0));
|
|
## disp(' ')
|
|
## disp('The 95% Q limit is')
|
|
## disp(q)
|
|
## if plots >= 1
|
|
## lim = [q q];
|
|
## plot(scl,res,scllim,lim,'--b')
|
|
## str = sprintf('Process Residual Q with 95 Percent Limit Based on %g PC Model',lv);
|
|
## title(str)
|
|
## xlabel('Sample Number')
|
|
## ylabel('Residual')
|
|
## pause
|
|
## end
|
|
## % Calculate T^2 limit using ftest routine
|
|
## if lv > 1
|
|
## if m > 300
|
|
## tsq = (lv*(m-1)/(m-lv))*ftest(.95,300,lv,2);
|
|
## else
|
|
## tsq = (lv*(m-1)/(m-lv))*ftest(.95,m-lv,lv,2);
|
|
## end
|
|
## disp(' ')
|
|
## disp('The 95% T^2 limit is')
|
|
## disp(tsq)
|
|
## % Calculate the value of T^2 by normalizing the scores to
|
|
## % unit variance and summing them up
|
|
## if plots >= 1.0
|
|
## temp2 = scores*inv(diag(ssq(1:lv,2).^.5));
|
|
## tsqvals = sum((temp2.^2)');
|
|
## tlim = [tsq tsq];
|
|
## plot(scl,tsqvals,scllim,tlim,'--b')
|
|
## str = sprintf('Value of T^2 with 95 Percent Limit Based on %g PC Model',lv);
|
|
## title(str)
|
|
## xlabel('Sample Number')
|
|
## ylabel('Value of T^2')
|
|
## end
|
|
## else
|
|
## disp('T^2 not calculated when number of latent variables = 1')
|
|
## tsq = 1.96^2;
|
|
## end
|
|
|