"Bayesian Source Separation: Beyond PCA and ICA" (Mohammad-Djafari, 2006) http://djafari.free.fr/pdf/esann06.pdf
Bayesian nonlinear independent component analysis by multi-layer perceptrons (Lappalainen and Honkela) https://www.hiit.fi/u/ahonkela/papers/ch7.pdf
Variational methods for Bayesian Independent Component Analysis (Choudrey, 2002, Oxford PhD thesis) http://www.robots.ox.ac.uk/~parg/pubs/theses/RizwanChoudrey_thesis.pdf
In [1]:
import numpy as np
import numpy.random as npr
import pylab as pl
%matplotlib inline
n = 100000
t = np.arange(n)
x1 = np.sin(t*0.001)
x1_noise = npr.randn(n)*0.1
x2 = np.cos(t*0.0001)
x2_noise = npr.randn(n)*0.01
S = np.c_[x1+x1_noise,x2+x2_noise]
S/= S.std(axis=0)
mix = np.array([[0.2,0.8,0.1],[0.1,0.8,0.2]])
obs = np.dot(S,mix)
In [142]:
y = S[:,0]+S[:,1]
In [150]:
pl.plot(y)
pl.figure()
pl.plot(S[:,0])
pl.figure()
pl.plot(S[:,1])
Out[150]:
In [121]:
from sklearn.decomposition import FastICA
In [156]:
ica = FastICA(n_components=2)
In [157]:
S_ = ica.fit_transform(obs)
In [158]:
mix_ = ica.mixing_
In [159]:
mix_
Out[159]:
In [160]:
pl.plot(S_[:,0])
pl.figure()
pl.plot(S_[:,1])
pl.figure()
pl.plot(S_[:,2])
In [66]:
In [296]:
p = 10
q = 5
N = 1000
X = npr.randn(N,p)
mu = np.mean(X,axis=0)
V = np.linalg.svd(X)[2][:q].T
lam = np.ones(q)
In [297]:
def reconstruction_error(X,mu,V,lam):
return np.sum(X - mu - V.dot(lam))
In [298]:
reconstruction_error(X,mu,V,lam)
Out[298]:
We can optimize for $\mu$ and $\lambda_i$ to obtain $$ \begin{align} \hat{\mu} & = \bar{x}\\ \hat{\lambda}_i & = \mathbf{V}_q^T (x_i - \bar{x}) \end{align} $$
In [299]:
lam = np.sum(V.dot(x - np.mean(x)),axis=0)
In [300]:
lam.shape
Out[300]:
In [301]:
U,D,V =np.linalg.svd(X.T)
In [302]:
D.shape
Out[302]:
In [303]:
PCs = (U*D).T
SVs = D
In [307]:
[pl.plot(PC) for PC in PCs];
In [305]:
pl.plot(SVs)
Out[305]:
In [247]:
np.sum(U.T.dot(U))
Out[247]:
In [333]:
U,D,Vt = np.linalg.svd(X.T,full_matrices=0)
D = np.diag(D)
N = len(X)
S = np.sqrt(N)*U
At = D*Vt / np.sqrt(N)
X = S * At
In [334]:
print(U.shape)
print(D.shape)
print(Vt.shape)
In [332]:
X.shape
Out[332]:
In [340]:
D.dot(Vt)
Out[340]:
In [350]:
U.conj().T == np.linalg.inv(U)
Out[350]:
In [321]:
D
Out[321]:
ktICA notes
tICA notes:
Also equivalent to the generalized eigenvalue problem: $$ C^{(\tau)} \vec{v} = \lambda \Sigma \vec{v}$$
In [437]:
def time_lag_corr_cov(X,tau=1):
mu = (X[:-tau].mean(0) + X[tau:].mean(0)) / 2
X_ = X-mu
M = len(X) - tau
dim = len(X.T)
corr = np.zeros((dim,dim))
cov = np.zeros((dim,dim))
for i in range(M):
corr += np.outer(X_[i],X_[i+tau]) + np.outer(X_[i+tau],X_[i])
cov += np.outer(X_[i],X_[i]) + np.outer(X_[i+tau],X_[i+tau])
return corr / (2.0*M),cov / (2.0*M)
In [438]:
X = np.ones((10000,100))
In [439]:
for i in range(10,len(X)):
X[i] = X[i-10] + npr.randn(X.shape[1])
In [440]:
pl.imshow(np.outer(X[10],X[11]),interpolation='none')
pl.colorbar()
Out[440]:
In [441]:
corr,cov = time_lag_corr_cov(X,1)
pl.imshow(corr,interpolation='none')
pl.colorbar()
Out[441]:
In [442]:
corr,cov = time_lag_corr_cov(X,1)
pl.imshow(corr,interpolation='none')
pl.colorbar()
pl.figure()
pl.imshow(cov,interpolation='none')
pl.colorbar()
Out[442]:
In [443]:
np.sum(abs(corr-cov))
Out[443]:
In [446]:
diff = []
for i in range(1,100):
corr,cov = time_lag_corr_cov(X,i)
diff.append(np.sum(abs(corr-cov)))
if i % 10 == 0:
print(i)
pl.plot(diff)
Out[446]:
In [447]:
pl.plot(np.log(diff))
Out[447]:
In [449]:
a,b = time_lag_corr_cov(X,100)
In [450]:
from scipy import linalg
In [456]:
w,v = linalg.eigh(a,b)
w.shape,v.shape
w = w[::-1]
v = v[:,::-1]
In [457]:
pl.plot(w)
Out[457]:
In [463]:
a,b = time_lag_corr_cov(df_megatraj_full[0],100)
w,v = linalg.eigh(a,b)
print(w.shape,v.shape)
w = w[::-1]
v = v[:,::-1]
In [464]:
pl.plot(w)
Out[464]:
In [466]:
pl.scatter(np.ones(len(w)),w,linewidth=0,alpha=0.5)
Out[466]:
In [488]:
pca = PCA(2)
X_pca = pca.fit_transform(df_megatraj_full[0])
pl.scatter(X_pca[:,0],X_pca[:,1],c=range(len(X_pca)),linewidths=0)
Out[488]:
In [467]:
from sklearn.decomposition import FastICA
In [479]:
ica = FastICA(3,fun='exp')
X_ = ica.fit_transform(df_megatraj_full[0])
In [480]:
pl.scatter(X_[:,0],X_[:,1],c=range(len(X_)),linewidths=0)
Out[480]:
Duh-- idea: perform time-aware clustering. You can visually see distinct clusters when they're colored by time
In [484]:
X_.shape
Out[484]:
In [487]:
pl.plot(abs(np.sum(X_[1:]-X_[:-1],1)))
Out[487]:
In [160]:
def gaussian(x,y,bandwidth=1.0):
return np.exp(-1.0/(2*bandwidth**2)*np.sum((x-y)**2))
In [168]:
%timeit gaussian(X[0],X[1])
In [177]:
type(gaussian) == type(np.exp)
Out[177]:
In [3]:
def time_lagged_K(kernel_matrix,tau=1):
# to-do: accept either a kernel function or a precomputed kernel matrix
# for now, just accept a precomputed kernel matrix
M = len(kernel_matrix) - tau
K = np.zeros((2*M,2*M))
#for i in range(M):
#for j in range(M):
#K[i,j] = kernel(X[i],X[j])
#K[i,j+M] = kernel(X[i],X[j+tau])
#K[i+M,j] = kernel(X[i+tau],X[j])
#K[i+M,j+M] = kernel(X[i+tau],X[j+tau])
K[:M,:M] = kernel_matrix[:M,:M]
K[:M,M:] = kernel_matrix[:M,tau:]
K[M:,:M] = kernel_matrix[tau:,:M]
K[M:,M:] = kernel_matrix[tau:,tau:]
return K
In [382]:
from msmbuilder.example_datasets import AlanineDipeptide,FsPeptide
#from msmbuilder.featurizer import SuperposeFeaturizer
print(AlanineDipeptide.description())
dataset = AlanineDipeptide().get()
ala_trajectories = dataset.trajectories
t = ala_trajectories[0]
import mdtraj
n=5000
rmsd = np.zeros((n,n))
for i in range(n):
rmsd[i,i:] = mdtraj.rmsd(t[i:n],t[i])
#if i % 100 == 0:
# print(i)
rmsd_ala = rmsd + rmsd.T
dataset = FsPeptide().get()
fs_trajectories = dataset.trajectories
fs_t = fs_trajectories[0]
rmsd = np.zeros((n,n))
fs_t = fs_trajectories[0]
for i in range(n):
rmsd[i,i:] = mdtraj.rmsd(fs_t[i:n],fs_t[i])
rmsd_fs = rmsd + rmsd.T
In [227]:
from msmbuilder import featurizer
rpf = featurizer.RawPositionsFeaturizer()
rpft = rpf.fit_transform(t.center_coordinates())
In [230]:
len(rpft),rpft[0].shape
Out[230]:
In [237]:
print(rpft[0].shape)
rpft[0][:10]
Out[237]:
In [239]:
np.array(rpft).shape
Out[239]:
In [241]:
np.reshape(rpft,(9999,22,3)).shape
Out[241]:
In [228]:
raw_pos = np.reshape(rpft,(9999,22,3))
raw_pos.shape
Out[228]:
In [260]:
from numpy.linalg import det
BC = lambda x,y: det(x.T.dot(y))/ np.sqrt(det(x.T.dot(x)) * det(y.T.dot(y)))
In [261]:
BC(raw_pos[10],raw_pos[100])
Out[261]:
In [262]:
n = 5000
In [263]:
kernel_matrix = np.zeros((n,n))
for i in range(n):
for j in range(n):
kernel_matrix[i,j] = BC(raw_pos[i],raw_pos[j])
if i % 500 == 0:
print(i)
In [307]:
kernel_matrix[0,0]
Out[307]:
In [309]:
kernel_matrix.shape
Out[309]:
In [163]:
kernel_matrix = 1-rmsd_ala
In [164]:
tau=5
M = len(kernel_matrix)-tau
from time import time
t = time()
K = time_lagged_K(kernel_matrix,tau=tau)
print(time() - t)
K.shape
Out[164]:
In [200]:
In [169]:
# earlier way to do it, which scaled as O((2n)**2)
from time import time
t = time()
K = time_lagged_K(X[:1000])
print(time() - t)
K.shape
Out[169]:
In [165]:
def construct_R(M):
''' R is a 2x2 block matrix with MxM entries, [[0,I],[I,0]]'''
R = np.zeros((2*M,2*M))
R[M:,:M] = np.eye(M)
R[:M,M:] = np.eye(M)
return R
In [166]:
R = construct_R(M)
R.shape
Out[166]:
In [167]:
K.shape
Out[167]:
In [168]:
KRK = (K.dot(R)).dot(K)
In [169]:
KK = K.dot(K)
In [170]:
evals,evecs = linalg.eigh(KRK,KK+0.01*np.eye(2*M))
In [171]:
evals[0],evals[-1]
Out[171]:
In [172]:
pl.scatter(range(len(evals)),evals,linewidth=0,alpha=0.5)
Out[172]:
In [173]:
sum(evals>0.1)
Out[173]:
In [174]:
evecs[-1]
Out[174]:
In [175]:
projection = np.zeros((len(rmsd_ala),2))
comp1,comp2 = (1,2)
for j in range(len(rmsd_ala)):
phi_t1 = 0
phi_t2 = 0
for i in range(M):
phi_t1 += evecs[comp1][i]*kernel_matrix[j,i] + evecs[comp1][i+M]*kernel_matrix[j,i+tau]
phi_t2 += evecs[comp2][i]*kernel_matrix[j,i] + evecs[comp2][i+M]*kernel_matrix[j,i+tau]
projection[j] = (phi_t1,phi_t2)
In [176]:
pl.scatter(projection[:,0],projection[:,1],c=range(len(projection)),
alpha=0.5,linewidth=0)
pl.colorbar()
Out[176]:
In [ ]:
In [139]:
from sklearn.manifold import SpectralEmbedding
se = SpectralEmbedding()
se.
In [146]:
X = npr.randn(100,10)
In [158]:
from scipy.spatial import distance
k_ = distance.squareform(distance.pdist(X))
In [183]:
rmsd_ala.shape
Out[183]:
In [264]:
from sklearn.decomposition import KernelPCA
kpca = KernelPCA(3,kernel='precomputed')
In [265]:
kpca_X = kpca.fit_transform(1-rmsd_ala)
In [266]:
pl.scatter(kpca_X[:,0],kpca_X[:,1],c=range(len(kpca_X)),linewidths=0)
Out[266]:
In [267]:
pl.scatter(kpca_X[:,0],kpca_X[:,1],c=dha[:len(kpca_X),3],linewidths=0)
Out[267]:
In [268]:
from sklearn.decomposition import KernelPCA
kpca = KernelPCA(3,kernel='precomputed')
kpca_X = kpca.fit_transform(kernel_matrix)
In [269]:
pl.scatter(kpca_X[:,0],kpca_X[:,1],c=dha[:len(kpca_X),3],linewidths=0)
Out[269]:
In [272]:
kernel_matrix.shape
Out[272]:
In [373]:
def kpca(kernel_matrix,n_components=2,c=None):
M=len(kernel_matrix)-1
evals,evecs=linalg.eigh(kernel_matrix,eigvals=(M-n_components,M))
order = sorted(range(len(evals)),key=lambda i:-evals[i])
alphas = evecs[:,order]
lambdas = evals[order]
pl.plot(range(1,len(lambdas)+1),lambdas)
sqrt_lambdas = np.diag(np.sqrt(lambdas))
if c==None:
c=range(M+1)
projection = alphas.dot(sqrt_lambdas)
pl.figure()
pl.scatter(projection[:,1],projection[:,2],c=c[:M+1],linewidths=0,alpha=0.5)
pl.figure()
pl.scatter(projection[:,2],projection[:,3],c=c[:M+1],linewidths=0,alpha=0.5)
pl.figure()
pl.scatter(projection[:,1],projection[:,3],c=c[:M+1],linewidths=0,alpha=0.5)
return lambdas,alphas
In [374]:
M=len(kernel_matrix[:100,:100])-1
n_components=10
evals,evecs=linalg.eigh(kernel_matrix[:100,:100],eigvals=(M-n_components,M))
In [375]:
evecs.shape,evals.shape
Out[375]:
In [376]:
lambdas_bc,alphas_bc=kpca(kernel_matrix[:5000,:5000],n_components=n_components,c=dha[:,3])
In [377]:
lambdas_rmsd,alphas_rmsd=kpca((1-rmsd_ala)[:5000,:5000],n_components=n_components,c=dha[:,3])
In [383]:
lambdas_rmsd,alphas_rmsd=kpca((1-rmsd_fs)[:5000,:5000],n_components=n_components)
In [401]:
thinned = [fs_trajectories[i][::30] for i in range(len(fs_trajectories))]
megatraj = thinned[0]
for thin in thinned[1:]:
megatraj += thin
megatraj
Out[401]:
In [413]:
megatraj_full = fs_trajectories[0]
for full_traj in fs_trajectories[1:]:
megatraj_full += full_traj
megatraj_full
Out[413]:
In [404]:
n = len(megatraj)
rmsd = np.zeros((n,n))
for i in range(n):
rmsd[i,i:] = mdtraj.rmsd(megatraj[i:n],megatraj[i])
rmsd_fs_m = rmsd + rmsd.T
In [405]:
lambdas_rmsd,alphas_rmsd=kpca((1-rmsd_fs_m),n_components=n_components)
In [387]:
lambdas_rmsd,alphas_rmsd=kpca((1-rmsd_fs_m)[:5000,:5000],n_components=n_components)
In [415]:
df_megatraj = df.fit_transform([megatraj])
df_megatraj_full = df.fit_transform([megatraj_full])
In [492]:
from msmbuilder.decomposition import tICA
tica = tICA(2,500)
X_ = tica.fit_transform(df_megatraj)[0]
X_full = tica.fit_transform(df_megatraj_full)[0]
In [493]:
pl.scatter(X_[:,0],X_[:,1],c=range(len(X_)),linewidths=0,alpha=0.5)
Out[493]:
In [494]:
pl.scatter(X_full[:,0],X_full[:,1],c=range(len(X_full)),linewidths=0,alpha=0.5)
Out[494]:
In [ ]:
In [389]:
sqrt_lambdas = np.diag(np.sqrt(lambdas_rmsd))
projection = alphas_rmsd.dot(sqrt_lambdas)
In [390]:
projection.shape
Out[390]:
In [398]:
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(projection[:,1],projection[:,2],projection[:,3],
alpha=0.5,c=range(len(projection)),cmap='Blues',linewidths=0)
Out[398]:
In [378]:
np.cumsum(lambdas_bc)/np.sum(lambdas_bc)
Out[378]:
In [379]:
np.cumsum(lambdas_rmsd)/np.sum(lambdas_rmsd)
Out[379]:
In [ ]:
In [300]:
vec.shape,val.shape
Out[300]:
In [290]:
kernel_matrix.shape
Out[290]:
In [ ]:
In [ ]:
In [ ]:
In [230]:
from sklearn.decomposition import PCA
pca = PCA(2)
pca_X = pca.fit_transform(raw_pos.reshape(9999,66))
In [243]:
pl.scatter(pca_X[:,0],pca_X[:,1],c=dha[:len(pca_X),3],linewidths=0)
Out[243]:
In [244]:
kpca_rbf = KernelPCA(2,kernel='rbf')
rbf_X = kpca_rbf.fit_transform(raw_pos.reshape(9999,66))
In [245]:
pl.scatter(rbf_X[:,0],rbf_X[:,1],c=dha[:len(rbf_X),3],linewidths=0)
Out[245]:
In [197]:
df = featurizer.DihedralFeaturizer()
dha = df.transform(t)
In [202]:
dha = np.reshape(np.array(dha),(9999,4))
In [252]:
kpca_rbf = KernelPCA(2,kernel='rbf')
rbf_X = kpca_rbf.fit_transform(dha)
In [253]:
pl.scatter(rbf_X[:,0],rbf_X[:,1],c=dha[:len(rbf_X),3],linewidths=0)
Out[253]:
In [216]:
dha[:,3].shape
Out[216]:
In [217]:
kpca_X[:len(dha),0].shape
Out[217]:
In [241]:
pl.scatter(dha[:,3],dha[:,1],c=range(len(dha)),linewidths=0)
Out[241]:
In [ ]:
Ideas:
Kernel choice:
Scaling up:
Kernel learning:
In [526]:
from sklearn.cross_decomposition import CCA
In [589]:
cca = CCA(n_components=2)
In [590]:
df_megatraj[0].shape
Out[590]:
In [624]:
t = 1
length=1000
cca.fit(df_megatraj[0][:length-t],df_megatraj[0][t:length])
Out[624]:
In [625]:
cca.score(df_megatraj[0][:1000-t],df_megatraj[0][t:1000])
Out[625]:
In [626]:
transformed = cca.transform(df_megatraj[0][:1000-t])
In [627]:
pl.scatter(transformed[:,0],transformed[:,1],c=range(len(transformed)),linewidths=0)
Out[627]:
In [628]:
transformed = cca.transform(df_megatraj_full[0])
In [629]:
pl.scatter(transformed[:,0],transformed[:,1],c=range(len(transformed)),alpha=0.5,linewidths=0)
Out[629]:
In [632]:
tica = tICA(2,1)
tica.fit([df_megatraj[0][:2000]])
transformed = tica.transform(df_megatraj_full)[0]
In [633]:
pl.scatter(transformed[:,0],transformed[:,1],c=range(len(transformed)),alpha=0.5,linewidths=0)
Out[633]:
In [2]:
from msmbuilder.example_datasets import AlanineDipeptide
ala = AlanineDipeptide().get()
In [5]:
print(ala.DESCR)
In [8]:
ala.trajectories
Out[8]:
In [ ]: