-
Notifications
You must be signed in to change notification settings - Fork 26
Description
Trying to get a poisson lds working with inputs has me encountering a problem within the newton step. I am running the following code:
npr.seed(0)
Parameters
D_obs = 10
D_latent = 3
D_input = 2
T = 5000
sigma_init = np.random.rand(D_latent, D_latent)/10
sigma_init = sigma_init @ sigma_init.T
sigma_states = np.random.rand(D_latent, D_latent)/10
sigma_states = sigma_states @ sigma_states.T
Simulate from a Poisson LDS
inputs = np.abs(np.random.randn(T, D_input)/10)
truemodel = DefaultPoissonLDS(D_obs, D_latent, D_input=D_input, sigma_init=sigma_init, sigma_states=sigma_states)
data, stateseq = truemodel.generate(T, inputs=inputs)
sigma_init = np.random.rand(D_latent, D_latent)/10
sigma_init = sigma_init @ sigma_init.T
sigma_states = np.random.rand(D_latent, D_latent)/10
sigma_states = sigma_states @ sigma_states.T
model = DefaultPoissonLDS(D_obs, D_latent, D_input=D_input, sigma_init=sigma_init, sigma_states=sigma_states)#, A=A, B=B, C=C)
model.add_data(data, inputs=inputs)
model.states_list[0].gaussian_states *= 0
N_iters = 20
def em_update(model):
model.EM_step()
ll = model.log_likelihood()
return ll
lls = [em_update(model) for _ in progprint_xrange(N_iters)]
Compute baseline likelihood under Poisson MLE model
rates = data.mean(0)
baseline = 0
for n in range(D_obs):
baseline += poisson(rates[n]).logpmf(data[:,n]).sum()
The code fails with the following error:
LinAlgError: 6470-th leading minor not positive definite
I believe the error has something to do with the cholesky decomposition or getting the inverse of the hessian in the newton step.