error with discrete observations in pyhsmm

[yarden]

I'm trying to run pyhsmm on discrete (categorical) observations. Sampling runs, but using the sampled models gives errors. Code:

Nmax = 25
# observations
obs_hypparams = {"K": 2,
                 "alpha_0": 1.}
# gamma prior on duration 
dur_hypparams = {'alpha_0': 1.,
                 'beta_0': 1.}
obs_distns = [distributions.Categorical(**obs_hypparams) for state in range(Nmax)]
dur_distns = [distributions.PoissonDuration(**dur_hypparams) for state in range(Nmax)]

posteriormodel = pyhsmm.models.WeakLimitHDPHSMM(
        alpha=6.,gamma=6.,
        init_state_concentration=6.,
        obs_distns=obs_distns,
        dur_distns=dur_distns)

# data
data = np.array([0, 1] * 10 + [0, 0] * 10)
posteriormodel.add_data(data)

# inference
models = []
for idx in progprint_xrange(150):
    posteriormodel.resample_model()
    if (idx+1) % 10 == 0:
        models.append(copy.deepcopy(posteriormodel))

# try to plot model (this fails)
fig = plt.figure()
for idx, model in enumerate(models):
    plt.clf()
    model.plot()
    plt.gcf().suptitle('HDP-HSMM sampled after %d iterations' % (10*(idx+1)))

# try to predict (this fails)
m = models[0]
m.predict(np.array([]), 5)

Prediction (m.predict) fails with the error:

    111     def predict(self,seed_data,timesteps,**kwargs):
--> 112         full_data = np.vstack((seed_data,np.nan*np.ones((timesteps,seed_data.shape[1]))))
    113         self.add_data(full_data,**kwargs)
    114         s = self.states_list.pop()

IndexError: tuple index out of range

and plotting fails with:

models.pyc in _plot_2d_data_scatter(self, ax, state_colors, plot_slice, update)
    308                 s._data_scatter.set_color(colorseq)
    309             else:
--> 310                 s._data_scatter = ax.scatter(data[:,0],data[:,1],c=colorseq,s=5)
    311             artists.append(s._data_scatter)
    312 

IndexError: too many indices for array

I don't care about using plotting; just using it to test if model worked. I mainly want to call predict on the sampled models. Any thoughts on how to run pyhsmm with discrete observations? Thanks!

[mattjj]

I think it's just the plotting that fails. The plotting pretty much only works with Gaussian emissions, so that code should really be separated out.

The model instance has a few plotting methods, and while plotting observations won't work with discrete emissions, plotting states should work. Can you try calling model.plot_statseseq(0)? (If that doesn't work, I can fix up that method to work independently, but you might want to write your own plotting code in the meantime.)

[mattjj]

Sorry, I read your question while distracted, so I missed the most important part: you want model.predict to work!

It looks like you need to pass it a 2D array. Can you try m.predict(np.array([[]]), 5)?

EDIT: 2D data not required anymore, and that was the wrong shape anyway, see below

[yarden]

Thanks, I tried playing with the array dimension and it didn't seem to work. When I try m.predict(np.array([[]]), 5) I get:

hmm_states.pyc in aBl(self)
    100             aBl = self._aBl = np.empty((data.shape[0],self.num_states))
    101             for idx, obs_distn in enumerate(self.obs_distns):
--> 102                 aBl[:,idx] = obs_distn.log_likelihood(data).ravel()
    103             aBl[np.isnan(aBl).any(1)] = 0.
    104 

multinomial.pyc in log_likelihood(self, x)
     95     def log_likelihood(self,x):
     96         err = np.seterr(divide='ignore')
---> 97         out = np.log(self.weights)[x]  # log(0) can happen, no warning
     98         np.seterr(**err)
     99         return out

IndexError: arrays used as indices must be of integer (or boolean) type

[mattjj]

Sorry, I guessed the wrong numpy call to make an empty array of the right shape. I think we want np.array([]).reshape(0, 1). I'll try to run your example to debug it.

EDIT: 2D data not required anymore

[mattjj]

Okay, got your script running. The problem is that Multinomial doesn't know to return zero potentials when called on np.nan arguments. The fix is in pybasicbayes (at least, one fix is). One sec...

[mattjj]

Can you try with pybasicbayes a2a4499 and pyhsmm 9b8b3f3? I just pushed those. The first should make Multinomial handle missing observations (nan observations) and the second should make HMM.predict handle 1D data.

[mattjj]

In particular, I tested your script with this call to predict:

m = models[0]
print m.predict(np.array([]), 100)

[yarden]

that works, thanks! For m.predict(np.array([]), 10) I get:

(array([ 0.,  0.,  1.,  0.,  0.,  0.,  0.,  1.,  0.,  1.]), array([20, 20,  9, 17, 17, 17, 17,  9,  9,  9], dtype=int32))

(I assume the second array is the index of the hidden states?)

[mattjj]

Exactly!

The predict code was a bit clunky because it was hastily written for a particular use case. It can probably be improved further (the seed data should be optional, for one). Let me know if you find other things that need fixing.