In [ ]:
from __future__ import division
import numpy as np
import numpy.random as npr
import matplotlib.pyplot as plt
from pybasicbayes.distributions import Regression
from pybasicbayes.util.text import progprint_xrange
from autoregressive.distributions import AutoRegression
import pyximport
pyximport.install()
import os
os.chdir('../pylds')
from models import LDS, DefaultLDS
from distributions import Regression_diag, AutoRegression_input
from obs_scheme import ObservationScheme
from user_util import gen_pars, rand_rotation_matrix, init_LDS_model, collect_LDS_stats
#npr.seed(0)
def update(model):
model.EM_step()
return model.log_likelihood()
#########################
# set some parameters #
#########################
n = 3
p = 10
T = 10000
pars_true, _ = gen_pars(n, p, u_dim=0,
pars_in=None,
obs_scheme=None,
gen_A='diagonal', lts=np.linspace(0.95, 0.98, n),
gen_B='random',
gen_Q='identity',
gen_mu0='random',
gen_V0='identity',
gen_C='random',
gen_d='scaled',
gen_R='fraction',
diag_R_flag=True,
x=None, y=None, u=None)
sub_pops = (np.arange(0, p//2 + 2), np.arange(p//2 - 2, p))
obs_pops = np.array((0,1))
obs_time = np.array((T//2,T))
obs_scheme = ObservationScheme(p, T, sub_pops, obs_pops, obs_time)
###################
# generate data #
###################
truemodel = LDS(
dynamics_distn=AutoRegression(A=pars_true['A'].copy(),sigma=pars_true['Q'].copy()),
emission_distn=Regression_diag(A=np.hstack((pars_true['C'].copy(), pars_true['d'].copy().reshape(p,1))),
sigma=pars_true['R'].copy(), affine=True),
)
truemodel.mu_init = pars_true['mu0'].copy()
truemodel.sigma_init = pars_true['V0'].copy()
data, stateseq = truemodel.generate(T)
###################
# EM cycles #
###################
pars_init, _ = gen_pars(n, p, u_dim=0,
pars_in=None,
obs_scheme=obs_scheme,
gen_A='diagonal', lts=0.99 * np.ones((n,)),
gen_B='random',
gen_Q='identity',
gen_mu0='random',
gen_V0='identity',
gen_C='PCA_subpop',
gen_d='mean',
gen_R='fractionObserved',
diag_R_flag=True,
x=stateseq.T, y=data.T, u=None)
# get E-step results for init pars
model = init_LDS_model(pars_init, data, obs_scheme) # set to initialisation
model.E_step()
stats_init,_ = collect_LDS_stats(model)
# get EM-step results after one iteration
model.M_step()
model.E_step()
stats_first,pars_first = collect_LDS_stats(model)
# get EM-step results after 50 iterations
model = init_LDS_model(pars_init, data, obs_scheme) # reset to initialisation
print 'fitting'
likes = [update(model) for _ in progprint_xrange(25)]
stats_hat,pars_hat = collect_LDS_stats(model)
# get EM-step results from true parameters
model = init_LDS_model(pars_true, data, obs_scheme) # reset to true pars
model.E_step()
stats_true,_ = collect_LDS_stats(model)
model.M_step()
###################
# store results #
###################
save_file = '../../../results/LDS_save'
from scipy.io import savemat # store results for comparison with Matlab code
from scipy.linalg import solve_discrete_lyapunov as dtlyap # solve discrete-time Lyapunov equation
save_file_m = {'x': model.states_list[0].stateseq,
'y': model.states_list[0].data,
'u' : [],
'll' : likes,
'T' : model.states_list[0].T,
'Trial': len(model.states_list),
'elapsedTime' : 0,
'ifUseB':False,
'ifUseA':True,
'epsilon':0,
'ifRDiagonal':False,
'covConvEps':0,
'truePars':pars_true,
'initPars':pars_init,
'firstPars':pars_first,
'estPars': pars_hat,
'stats_0': stats_init,
'stats_1': stats_first,
'stats_h': stats_hat,
'stats_true': stats_true,
'Pi':dtlyap(pars_true['A'], pars_true['Q']),
'Pi_h':dtlyap(pars_hat['A'], pars_hat['Q']),
'Pi_t':pars_true['A'].dot(dtlyap(pars_true['A'], pars_true['Q'])),
'Pi_t_h':pars_hat['A'].dot(dtlyap(pars_hat['A'], pars_hat['Q'])),
'obsScheme' : obs_scheme}
savemat(save_file,save_file_m) # does the actual saving
In [ ]:
from __future__ import division
import numpy as np
import numpy.random as npr
import matplotlib.pyplot as plt
from pybasicbayes.distributions import Regression
from pybasicbayes.util.text import progprint_xrange
from autoregressive.distributions import AutoRegression
import pyximport
pyximport.install()
import os
os.chdir('../pylds')
from models import LDS, DefaultLDS
from distributions import Regression_diag, AutoRegression_input
from obs_scheme import ObservationScheme
from user_util import gen_pars, rand_rotation_matrix, init_LDS_model, collect_LDS_stats
#npr.seed(0)
def update(model):
model.EM_step()
return model.log_likelihood()
#########################
# set some parameters #
#########################
n = 10
from scipy.io import loadmat
data =loadmat('/home/mackelab/Desktop/Projects/Stitching/results/cosyne_poster/experiment_1/test_data')
data = data['data'].astype(np.float)
T,p = data.shape
#tmp = np.random.choice(np.arange(p),size=p,replace=False)
#sub_pops = (np.sort(tmp[:p//2 + 2]), np.sort(tmp[p//2 - 2:]))
sub_pops = (np.arange(p), np.arange(p))
obs_pops = np.array((0,1))
obs_time = np.array((T//2,T))
obs_scheme = ObservationScheme(p, T, sub_pops, obs_pops, obs_time)
###################
# EM cycles #
###################
pars_init, _ = gen_pars(n, p, u_dim=0,
pars_in=None,
obs_scheme=obs_scheme,
gen_A='diagonal', lts=0.99 * np.ones((n,)),
gen_B='random',
gen_Q='identity',
gen_mu0='random',
gen_V0='identity',
gen_C='PCA',
gen_d='mean',
gen_R='fractionObserved',
diag_R_flag=True,
x=None, y=data.T, u=None)
# get E-step results for init pars
model = init_LDS_model(pars_init, data, obs_scheme) # set to initialisation
model.E_step()
stats_init,_ = collect_LDS_stats(model)
# get EM-step results after one iteration
model.M_step()
model.E_step()
stats_first,pars_first = collect_LDS_stats(model)
# get EM-step results after 50 iterations
model = init_LDS_model(pars_init, data, obs_scheme) # reset to initialisation
print 'fitting'
likes = [update(model) for _ in progprint_xrange(50)]
stats_hat,pars_hat = collect_LDS_stats(model)
In [ ]:
from scipy.io import loadmat # store results for comparison with Matlab code
from __future__ import division
import numpy as np
import numpy.random as npr
import matplotlib.pyplot as plt
from pybasicbayes.distributions import Regression
from pybasicbayes.util.text import progprint_xrange
from autoregressive.distributions import AutoRegression
import pyximport
pyximport.install()
import os
os.chdir('../pylds')
from models import LDS, DefaultLDS
from distributions import Regression_diag, AutoRegression_input
from obs_scheme import ObservationScheme
from user_util import gen_pars, rand_rotation_matrix, init_LDS_model, collect_LDS_stats
#npr.seed(0)
def update(model):
model.EM_step()
return model.log_likelihood()
spikes = loadmat('/home/mackelab/Desktop/Projects/Stitching/results/cosyne_poster/gb_net/spikes_10trials_10msBins')
spikes= spikes['spikes_out']
p, T = 1000, 600
#data = spikes[0][0].T
idx_n = np.sort(np.random.choice(1000, size=p, replace=False))
data = np.vstack([spikes[i][0].T[:,idx_n] for i in range(spikes.size)]).astype(np.float)
T *= spikes.size
n = 10
#tmp = np.random.choice(np.arange(p),size=p,replace=False)
#sub_pops = (np.sort(tmp[:p//2 + 2]), np.sort(tmp[p//2 - 2:]))
sub_pops = (np.arange(p//2+50), np.arange(p//2-50, p))
obs_pops = np.array((0,1))
obs_time = np.array((T//2,T))
obs_scheme = ObservationScheme(p, T, sub_pops, obs_pops, obs_time)
###################
# EM cycles #
###################
loadfile = np.load('/home/mackelab/Desktop/Projects/Stitching/results/cosyne_poster/experiment_1/init/init_experiment.npz')
initkey = 'params_naive'
pars_init = loadfile[initkey].reshape(1,)[0]
tmp = {}
for j in range(len(pars_init.keys())):
tmp[str(pars_init.keys()[j])] = pars_init.values()[j]
pars_init = tmp
pars_init['R'] = np.diag(pars_init['R'])
pars_init['V0'] = pars_init['Q']
#D, V = np.linalg.eig(pars_init['A'])
#if np.any(np.abs(D) > 1):
# print(np.abs(D))
# D /= np.maximum(np.abs(D), 10)
# print(np.abs(D))
# pars_init['A'] = np.real(V.dot(np.diag(D).dot(np.linalg.inv(V))))
model = init_LDS_model(pars_init, data, obs_scheme) # set to initialisation
#stats_init,_ = collect_LDS_stats(model)
#print 'fitting'
#likes = [update(model) for _ in progprint_xrange(100)]
#stats_hat,pars_hat = collect_LDS_stats(model)
In [ ]:
spikes = loadmat('/home/mackelab/Desktop/Projects/Stitching/results/cosyne_poster/gb_net/spikes_10trials_10msBins')
spikes= spikes['spikes_out']
p, T = 1000, 6000
#data = spikes[0][0].T
idx_n = np.sort(np.random.choice(1000, size=p, replace=False))
data = np.vstack([spikes[i][0].T[:,idx_n] for i in range(spikes.size)]).astype(np.float)
T *= spikes.size
n = 10
loadfile = np.load('/home/mackelab/Desktop/Projects/Stitching/results/cosyne_poster/experiment_1/init/init_experiment.npz')
initkey = 'params_naive'
pars_init = loadfile[initkey].reshape(1,)[0]
tmp = {}
for j in range(len(pars_init.keys())):
tmp[str(pars_init.keys()[j])] = pars_init.values()[j]
pars_init = tmp
pars_init['R'] = np.diag(pars_init['R'])
pars_init['V0'] = pars_init['Q']
pars_hat = pars_init.copy()
pars_hat['R'] = np.diag(pars_hat['R'])
import scipy as sp
import matplotlib.pyplot as plt
#Pi = sp.linalg.solve_discrete_lyapunov(pars_hat['A'], pars_hat['Q'])
Pi = pars_hat['Pi']
plt.figure(1,figsize=(15,22))
try:
Pi_true = sp.linalg.solve_discrete_lyapunov(pars_true['A'], pars_true['Q'])
tmp1 = pars_true['C'].dot(Pi_true.dot(pars_true['C'].transpose())) + np.diag(pars_true['R'])
except:
cov_all = np.cov(np.hstack([data[1:], data[:-1]]).T)
tmp1 = cov_all[np.ix_(np.arange(p), np.arange(p))]
tmp2 = pars_hat['C'].dot(Pi.dot(pars_hat['C'].transpose())) + pars_hat['R']
m = np.min((tmp1-np.diag(np.diag(tmp1))).min(),(tmp2-np.diag(np.diag(tmp2))).min())
M = np.max((tmp1-np.diag(np.diag(tmp1))).max(),(tmp2-np.diag(np.diag(tmp2))).max())
plt.subplot(2,3,1)
plt.imshow(tmp1-np.diag(np.diag(tmp1)),interpolation='none')
plt.clim(m,M)
plt.colorbar()
plt.title('true instantaneous covs')
plt.subplot(2,3,2)
plt.imshow(tmp2-np.diag(np.diag(tmp2)), interpolation='none')
plt.clim(m,M)
plt.colorbar()
plt.title('estimated instantaneous covs')
plt.subplot(2,3,3)
plt.plot(tmp1[:], tmp2[:], '.')
plt.xlabel('true')
plt.ylabel('est')
try:
tmp1 = pars_true['C'].dot(pars_true['A']).dot(Pi_true.dot(pars_true['C'].transpose()))
except:
tmp1 = cov_all[np.ix_(np.arange(0, p), np.arange(p+1 , 2*p))]
tmp2 = pars_hat['C'].dot(pars_hat['A']).dot(Pi.dot(pars_hat['C'].transpose()))
m = np.min(tmp1.min(),tmp2.min())
M = np.max(tmp1.max(),tmp2.max())
plt.subplot(2,3,4)
plt.imshow(tmp1,interpolation='none')
plt.clim(m,M)
plt.colorbar()
plt.title('true time_lagged covs')
plt.subplot(2,3,5)
plt.imshow(tmp2, interpolation='none')
plt.clim(m,M)
plt.colorbar()
plt.title('estimated time-lagged covs')
plt.subplot(2,3,6)
plt.plot(tmp1[:], tmp2[:], '.')
plt.xlabel('true')
plt.ylabel('est')
plt.figure(2,figsize=(15,15))
plt.subplot(2,2,1)
plt.plot(pars_hat['d'])
try:
plt.plot(pars_true['d'])
except:
pass
plt.legend(['true', 'est'])
plt.title('d')
plt.subplot(2,2,2)
plt.plot(pars_hat['R'])
try:
plt.plot(pars_true['R'])
except:
pass
plt.legend(['true', 'est'])
plt.title('R')
plt.subplot(2,2,3)
plt.plot(np.sort(np.linalg.eig(pars_hat['A'])[0]))
try:
plt.plot(np.sort(np.linalg.eig(pars_true['A'])[0]))
except:
pass
plt.legend(['true', 'est'])
plt.title('eig(A)')
plt.title('eigenvalues of A')
plt.show()
""" second batch of figures"""
covy_h= np.dot( np.dot(pars_hat['C'], Pi),pars_hat['C'].transpose()) + pars_hat['R']
try:
covy_t= np.dot(np.dot(pars_true['C'], Pi_true),pars_true['C'].transpose()) + np.diag(pars_true['R'])
covy_tl_t=(np.dot(np.dot(pars_true['C'],np.dot(pars_true['A'], Pi_true)),pars_true['C'].transpose()))
plot_truth = True
except:
plot_truth = False
y_tl = np.zeros([2*p,T-1])
y_tl[range(p),:] = data[range(0,T-1),:].T
y_tl[range(p,2*p),:] = data[range(1,T),:].T
covy = np.cov(y_tl)
covy_e= covy[np.ix_(range(p),range(p))]
covy_tl_e= covy[np.ix_(range(p,2*p),range(0,p))]
covy_tl_h= np.dot(np.dot(pars_hat['C'], np.dot(pars_hat['A'],Pi)), pars_hat['C'].transpose())
idx_stitched = np.ones([p,p],dtype = bool)
for i in range(len(obs_scheme.sub_pops)):
if len(obs_scheme.sub_pops[i])>0:
idx_stitched[np.ix_(obs_scheme.sub_pops[i],obs_scheme.sub_pops[i])] = False
plt.imshow(idx_stitched,interpolation='none')
plt.figure(3, figsize=(20,10))
plt.subplot(1,3,1)
if plot_truth:
plt.plot(covy_e[np.invert(idx_stitched)], covy_t[np.invert(idx_stitched)], '.')
plt.title('obs. emp vs. obs. true')
plt.ylabel('instantaneous')
plt.subplot(1,3,2)
plt.plot(covy_e[np.invert(idx_stitched)], covy_h[np.invert(idx_stitched)], '.')
plt.title('obs. emp vs. obs. stitched')
plt.subplot(1,3,3)
if plot_truth:
plt.plot(covy_t[np.invert(idx_stitched)], covy_h[np.invert(idx_stitched)], '.')
plt.title('obs. true vs. obs. stitched')
plt.figure(4, figsize=(20,10))
plt.subplot(1,3,1)
if plot_truth:
plt.plot(covy_tl_e[np.invert(idx_stitched)], covy_tl_t[np.invert(idx_stitched)], '.')
plt.ylabel('time-lagged')
plt.title('obs. emp vs. obs. true')
plt.subplot(1,3,2)
plt.plot(covy_tl_e[np.invert(idx_stitched)], covy_tl_h[np.invert(idx_stitched)], '.')
plt.title('obs. emp vs. obs. stitched')
plt.subplot(1,3,3)
if plot_truth:
plt.plot(covy_tl_t[np.invert(idx_stitched)], covy_tl_h[np.invert(idx_stitched)], '.')
plt.title('obs. non-observed true vs. obs. stitched')
plt.figure(5, figsize=(20,10))
plt.subplot(1,3,1)
if plot_truth:
plt.plot(covy_e[idx_stitched], covy_t[idx_stitched], '.')
plt.title('non-obs. emp vs. non-obs. true')
plt.subplot(1,3,2)
plt.plot(covy_e[idx_stitched], covy_h[idx_stitched], '.')
plt.title('non-obs. emp vs. non-obs. stitched')
plt.subplot(1,3,3)
if plot_truth:
plt.plot(covy_t[idx_stitched], covy_h[idx_stitched], '.')
plt.title('non-obs. true vs. non-obs. titched')
plt.figure(6, figsize=(20,10))
plt.subplot(1,3,1)
if plot_truth:
plt.plot(covy_tl_e[idx_stitched], covy_tl_t[idx_stitched], '.')
plt.title('non-obs. emp vs. non-obs. true')
plt.subplot(1,3,2)
plt.plot(covy_tl_e[idx_stitched], covy_tl_h[idx_stitched], '.')
plt.title('non-obs. emp vs. non-obs. stitched')
plt.subplot(1,3,3)
if plot_truth:
plt.plot(covy_tl_t[idx_stitched], covy_tl_h[idx_stitched], '.')
plt.title('non-obs. true vs. non-obs. stitched')
print('corr(y_t,y_t) stitchted: ', np.corrcoef(covy_e[idx_stitched], covy_h[idx_stitched]))
print('corr(y_t,y_t) observed: ', np.corrcoef(covy_e[np.invert(idx_stitched)], covy_h[np.invert(idx_stitched)]))
print('corr(y_t,y_t-1) stitchted: ', np.corrcoef(covy_tl_e[idx_stitched], covy_tl_h[idx_stitched]))
print('corr(y_t,y_t-1) observed: ', np.corrcoef(covy_tl_e[np.invert(idx_stitched)], covy_tl_h[np.invert(idx_stitched)]))
In [ ]:
loadfile = np.load('/home/mackelab/Desktop/Projects/Stitching/results/cosyne_poster/experiment_1/fits/params_naive_p1000_iter700.npz')
pars_hat = loadfile['estPars'].reshape(1,)[0]
likes = loadfile['ll']
import scipy as sp
import matplotlib.pyplot as plt
Pi = sp.linalg.solve_discrete_lyapunov(pars_hat['A'], pars_hat['Q'])
plt.figure(1,figsize=(15,22))
try:
Pi_true = sp.linalg.solve_discrete_lyapunov(pars_true['A'], pars_true['Q'])
tmp1 = pars_true['C'].dot(Pi_true.dot(pars_true['C'].transpose())) + np.diag(pars_true['R'])
except:
cov_all = np.cov(np.hstack([data[1:], data[:-1]]).T)
tmp1 = cov_all[np.ix_(np.arange(p), np.arange(p))]
tmp2 = pars_hat['C'].dot(Pi.dot(pars_hat['C'].transpose())) + pars_hat['R']
m = np.min((tmp1-np.diag(np.diag(tmp1))).min(),(tmp2-np.diag(np.diag(tmp2))).min())
M = np.max((tmp1-np.diag(np.diag(tmp1))).max(),(tmp2-np.diag(np.diag(tmp2))).max())
plt.subplot(2,3,1)
plt.imshow(tmp1-np.diag(np.diag(tmp1)),interpolation='none')
plt.clim(m,M)
plt.colorbar()
plt.title('true instantaneous covs')
plt.subplot(2,3,2)
plt.imshow(tmp2-np.diag(np.diag(tmp2)), interpolation='none')
plt.clim(m,M)
plt.colorbar()
plt.title('estimated instantaneous covs')
plt.subplot(2,3,3)
plt.plot(tmp1[:], tmp2[:], '.')
plt.xlabel('true')
plt.ylabel('est')
try:
tmp1 = pars_true['C'].dot(pars_true['A']).dot(Pi_true.dot(pars_true['C'].transpose()))
except:
tmp1 = cov_all[np.ix_(np.arange(0, p), np.arange(p+1 , 2*p))]
tmp2 = pars_hat['C'].dot(pars_hat['A']).dot(Pi.dot(pars_hat['C'].transpose()))
m = np.min(tmp1.min(),tmp2.min())
M = np.max(tmp1.max(),tmp2.max())
plt.subplot(2,3,4)
plt.imshow(tmp1,interpolation='none')
plt.clim(m,M)
plt.colorbar()
plt.title('true time_lagged covs')
plt.subplot(2,3,5)
plt.imshow(tmp2, interpolation='none')
plt.clim(m,M)
plt.colorbar()
plt.title('estimated time-lagged covs')
plt.subplot(2,3,6)
plt.plot(tmp1[:], tmp2[:], '.')
plt.xlabel('true')
plt.ylabel('est')
plt.figure(2,figsize=(15,15))
plt.subplot(2,2,1)
plt.plot(pars_hat['d'])
try:
plt.plot(pars_true['d'])
except:
pass
plt.legend(['true', 'est'])
plt.title('d')
plt.subplot(2,2,2)
plt.plot(pars_hat['R'])
try:
plt.plot(pars_true['R'])
except:
pass
plt.legend(['true', 'est'])
plt.title('R')
plt.subplot(2,2,3)
plt.plot(np.sort(np.linalg.eig(pars_hat['A'])[0]))
try:
plt.plot(np.sort(np.linalg.eig(pars_true['A'])[0]))
except:
pass
plt.legend(['true', 'est'])
plt.title('eig(A)')
plt.title('eigenvalues of A')
plt.show()
""" second batch of figures"""
covy_h= np.dot( np.dot(pars_hat['C'], Pi),pars_hat['C'].transpose()) + pars_hat['R']
try:
covy_t= np.dot(np.dot(pars_true['C'], Pi_true),pars_true['C'].transpose()) + np.diag(pars_true['R'])
covy_tl_t=(np.dot(np.dot(pars_true['C'],np.dot(pars_true['A'], Pi_true)),pars_true['C'].transpose()))
plot_truth = True
except:
plot_truth = False
y_tl = np.zeros([2*p,T-1])
y_tl[range(p),:] = data[range(0,T-1),:].T
y_tl[range(p,2*p),:] = data[range(1,T),:].T
covy = np.cov(y_tl)
covy_e= covy[np.ix_(range(p),range(p))]
covy_tl_e= covy[np.ix_(range(p,2*p),range(0,p))]
covy_tl_h= np.dot(np.dot(pars_hat['C'], np.dot(pars_hat['A'],Pi)), pars_hat['C'].transpose())
idx_stitched = np.ones([p,p],dtype = bool)
for i in range(len(obs_scheme.sub_pops)):
if len(obs_scheme.sub_pops[i])>0:
idx_stitched[np.ix_(obs_scheme.sub_pops[i],obs_scheme.sub_pops[i])] = False
plt.imshow(idx_stitched,interpolation='none')
plt.figure(3, figsize=(20,10))
plt.subplot(1,3,1)
if plot_truth:
plt.plot(covy_e[np.invert(idx_stitched)], covy_t[np.invert(idx_stitched)], '.')
plt.title('obs. emp vs. obs. true')
plt.ylabel('instantaneous')
plt.subplot(1,3,2)
plt.plot(covy_e[np.invert(idx_stitched)], covy_h[np.invert(idx_stitched)], '.')
plt.title('obs. emp vs. obs. stitched')
plt.subplot(1,3,3)
if plot_truth:
plt.plot(covy_t[np.invert(idx_stitched)], covy_h[np.invert(idx_stitched)], '.')
plt.title('obs. true vs. obs. stitched')
plt.figure(4, figsize=(20,10))
plt.subplot(1,3,1)
if plot_truth:
plt.plot(covy_tl_e[np.invert(idx_stitched)], covy_tl_t[np.invert(idx_stitched)], '.')
plt.ylabel('time-lagged')
plt.title('obs. emp vs. obs. true')
plt.subplot(1,3,2)
plt.plot(covy_tl_e[np.invert(idx_stitched)], covy_tl_h[np.invert(idx_stitched)], '.')
plt.title('obs. emp vs. obs. stitched')
plt.subplot(1,3,3)
if plot_truth:
plt.plot(covy_tl_t[np.invert(idx_stitched)], covy_tl_h[np.invert(idx_stitched)], '.')
plt.title('obs. non-observed true vs. obs. stitched')
plt.figure(5, figsize=(20,10))
plt.subplot(1,3,1)
if plot_truth:
plt.plot(covy_e[idx_stitched], covy_t[idx_stitched], '.')
plt.title('non-obs. emp vs. non-obs. true')
plt.subplot(1,3,2)
plt.plot(covy_e[idx_stitched], covy_h[idx_stitched], '.')
plt.title('non-obs. emp vs. non-obs. stitched')
plt.subplot(1,3,3)
if plot_truth:
plt.plot(covy_t[idx_stitched], covy_h[idx_stitched], '.')
plt.title('non-obs. true vs. non-obs. titched')
plt.figure(6, figsize=(20,10))
plt.subplot(1,3,1)
if plot_truth:
plt.plot(covy_tl_e[idx_stitched], covy_tl_t[idx_stitched], '.')
plt.title('non-obs. emp vs. non-obs. true')
plt.subplot(1,3,2)
plt.plot(covy_tl_e[idx_stitched], covy_tl_h[idx_stitched], '.')
plt.title('non-obs. emp vs. non-obs. stitched')
plt.subplot(1,3,3)
if plot_truth:
plt.plot(covy_tl_t[idx_stitched], covy_tl_h[idx_stitched], '.')
plt.title('non-obs. true vs. non-obs. stitched')
print('corr(y_t,y_t) stitchted: ', np.corrcoef(covy_e[idx_stitched], covy_h[idx_stitched]))
print('corr(y_t,y_t) observed: ', np.corrcoef(covy_e[np.invert(idx_stitched)], covy_h[np.invert(idx_stitched)]))
print('corr(y_t,y_t-1) stitchted: ', np.corrcoef(covy_tl_e[idx_stitched], covy_tl_h[idx_stitched]))
print('corr(y_t,y_t-1) observed: ', np.corrcoef(covy_tl_e[np.invert(idx_stitched)], covy_tl_h[np.invert(idx_stitched)]))
In [ ]:
%matplotlib inline
from scipy.io import loadmat # store results for comparison with Matlab code
from __future__ import division
import numpy as np
import numpy.random as npr
import matplotlib.pyplot as plt
from pybasicbayes.distributions import Regression
from pybasicbayes.util.text import progprint_xrange
from autoregressive.distributions import AutoRegression
import pyximport
pyximport.install()
import os
os.chdir('../pylds')
from models import LDS, DefaultLDS
from distributions import Regression_diag, AutoRegression_input
from obs_scheme import ObservationScheme
from user_util import gen_pars, rand_rotation_matrix, init_LDS_model, collect_LDS_stats
#npr.seed(0)
def update(model):
model.EM_step()
return model.log_likelihood()
print('loading data')
spikes = loadmat('/home/mackelab/Desktop/Projects/Stitching/results/cosyne_poster/gb_net/spikes_10trials_10msBins')
spikes= spikes['spikes_out']
print('concatenating data')
p, T = 1000, 600
#data = spikes[0][0].T
idx_n = np.sort(np.random.choice(1000, size=p, replace=False))
data = np.vstack([spikes[i][0].T[:,idx_n] for i in range(spikes.size)]).astype(np.float)
T *= spikes.size
n = 10
###################
# EM cycles #
###################
print('loading parameters')
loadfile = np.load('/home/mackelab/Desktop/Projects/Stitching/results/cosyne_poster/experiment_1/fits/params_naive_p1000_iter500.npz')
pars_init = loadfile['estPars'].reshape(1,)[0]
pars_init['R'] = np.diag(pars_init['R'])
sub_pops = tuple(loadfile['sub_pops'].tolist())
obs_pops =loadfile['obs_pops']
obs_time = loadfile['obs_time']
obs_scheme = ObservationScheme(p, T, sub_pops, obs_pops, obs_time)
model = init_LDS_model(pars_init, data, obs_scheme) # set to initialisation
#stats_init,_ = collect_LDS_stats(model)
print('fitting')
likes = [update(model) for _ in progprint_xrange(200)]
stats_hat,pars_hat = collect_LDS_stats(model)
In [ ]:
likes
In [ ]:
from scipy.io import savemat
broken = False
save_file = '/home/mackelab/Desktop/Projects/Stitching/results/cosyne_poster/experiment_1/fits/params_naive_p1000_iter700'
eps = - np.inf
save_file_m = {'ifBroken':broken,
'll' : likes,
'T' : model.states_list[0].T,
'Trial': len(model.states_list),
'epsilon':eps,
'initPars':pars_init,
'estPars': pars_hat,
'stats_h': stats_hat,
'sub_pops' : obs_scheme.sub_pops,
'obs_pops' : obs_scheme.obs_pops,
'obs_time' : obs_scheme.obs_time}
savemat(save_file,save_file_m) # does the actual saving
np.savez(save_file,
broken=broken,
ll=likes,
T=model.states_list[0].T,
Trial=len(model.states_list),
epsilon=eps,
initPars=pars_init,
estPars =pars_hat,
stats_h = stats_hat,
sub_pops=obs_scheme.sub_pops,
obs_time=obs_scheme.obs_time,
obs_pops=obs_scheme.obs_pops)
In [ ]:
broken=False
In [ ]:
from __future__ import division
import numpy as np
import numpy.random as npr
import matplotlib.pyplot as plt
from pybasicbayes.distributions import Regression
from pybasicbayes.util.text import progprint_xrange
from autoregressive.distributions import AutoRegression
import pyximport
pyximport.install()
import os
os.chdir('../pylds')
from models import LDS, DefaultLDS
from distributions import Regression_diag, AutoRegression_input
from obs_scheme import ObservationScheme
from user_util import gen_pars, rand_rotation_matrix, init_LDS_model, collect_LDS_stats
#npr.seed(0)
def update(model):
model.EM_step()
return model.log_likelihood()
#########################
# set some parameters #
#########################
n = 3
p = 20
T = 10000
num_sets = 100;
for idx_d in range(num_sets):
print('generating set #', idx_d)
pars_true, _ = gen_pars(n, p, u_dim=0,
pars_in=None,
obs_scheme=None,
gen_A='full', lts=np.linspace(0.95, 0.98, n),
gen_B='random',
gen_Q='identity',
gen_mu0='random',
gen_V0='identity',
gen_C='random',
gen_d='scaled',
gen_R='fraction',
diag_R_flag=True,
x=None, y=None, u=None)
sub_pops = (np.arange(0, p//2 + 2), np.arange(p//2 - 2, p))
obs_pops = np.array((0,1))
obs_time = np.array((T//2,T))
obs_scheme = ObservationScheme(p, T, sub_pops, obs_pops, obs_time)
###################
# generate data #
###################
truemodel = LDS(
dynamics_distn=AutoRegression(A=pars_true['A'].copy(),sigma=pars_true['Q'].copy()),
emission_distn=Regression_diag(A=np.hstack((pars_true['C'].copy(), pars_true['d'].copy().reshape(p,1))),
sigma=pars_true['R'].copy(), affine=True),
)
truemodel.mu_init = pars_true['mu0'].copy()
truemodel.sigma_init = pars_true['V0'].copy()
data, stateseq = truemodel.generate(T)
save_file = '../../../results/cosyne_poster/debug/LDS_save_n' + str(n) + '_p' + str(p) + '_T' +str(T) + '_idx'+ str(idx_d)
from scipy.io import savemat # store results for comparison with Matlab code
from scipy.linalg import solve_discrete_lyapunov as dtlyap # solve discrete-time Lyapunov equation
save_file_m = {'x': truemodel.states_list[0].stateseq,
'y': truemodel.states_list[0].data,
'u' : [],
'T' : truemodel.states_list[0].T,
'Trial': len(truemodel.states_list),
'truePars':pars_true,
'obsScheme' : obs_scheme}
savemat(save_file,save_file_m) # does the actual saving
In [ ]: