pymc3 multinomial mixture gets stuck

[antjoseme]

I am trying use PYMC3 to implement an example where the data comes from a mixture of multinomials. The goal is to infer the underlying state_prob vector (see below). The code runs, but the Metropolis sampler gets stuck at the initial state_prior vector. Also, for some reason I have not been able to get NUTS to work.

import numpy as np
from pymc3 import Model, Multinomial, Dirichlet, Uniform
import pymc3 as pm
import theano.tensor as tt

N = 20
Nstates = 3
Nlabel = Nstates
state_prob = np.array([0.3, 0.1, 0.6]) #need to infer this
state_data = np.random.choice(Nstates, N, p=state_prob)
state_label_tran = np.array([[0.3, 0.2, 0.5],
                               [0.1, 0.5, 0.4],
                               [0.0, 0.05, 0.95]])
label_prob_given_state = state_label_tran[state_data]
def rand_mult(row_p):
    return np.random.multinomial(1, row_p)
label_data = np.apply_along_axis(rand_mult, 1, label_prob_given_state)
#label_data = np.reshape(label_data, (-1,N))[0]
print 'Unobserved state data'
print state_data
print 'Observed label data'
print label_data

with Model() as simple_dag:
    state = Dirichlet('state', np.array([1.1, 1.1, 1.1]), shape=3)
    label_dist = [pm.Multinomial.dist(p=state_label_tran[i], n=1, shape=3) for i in range(3)]
    label = pm.Mixture('label', w=state, comp_dists=label_dist, observed=label_data)
    
with simple_dag:
    trace = pm.sample(1000, cores=2, chains=2, tune=1000, progressbar=True)

[fonnesbeck]

You are probably better off reporting the NITS issues, and getting that to work. Metropolis won’t work as well, and is liable to get stuck, as it appears to be doing.

[antjoseme]

@fonnesbeck
When I use NUTS (using the following commands)

with simple_dag:
    step = pymc3.NUTS(vars=[state])
    trace = pymc3.sample(10000, cores=2, chains=2, tune=500, step=step, progressbar=True)

I get a "Bad initial energy" error.

...
methods.hmc.nuts.NUTS object>, start={'state_stickbreaking__': array([0., 0.])}, trace=None, chain=0, tune=500, model=<pymc3.model.Model object>, random_seed=315794368)
    647         step.tune = bool(tune)
    648         for i in range(draws):
    649             if i == tune:
    650                 step = stop_tuning(step)
    651             if step.generates_stats:
--> 652                 point, states = step.step(point)
        point = {'state_stickbreaking__': array([0., 0.])}
        states = undefined
        step.step = <bound method NUTS.step of <pymc3.step_methods.hmc.nuts.NUTS object>>
    653                 if strace.supports_sampler_stats:
    654                     strace.record(point, states)
    655                 else:
    656                     strace.record(point)

...........................................................................
/Users/ajoseph/anaconda2/lib/python2.7/site-packages/pymc3/step_methods/arraystep.py in step(self=<pymc3.step_methods.hmc.nuts.NUTS object>, point={'state_stickbreaking__': array([0., 0.])})
    217     def step(self, point):
    218         self._logp_dlogp_func.set_extra_values(point)
    219         array = self._logp_dlogp_func.dict_to_array(point)
    220 
    221         if self.generates_stats:
--> 222             apoint, stats = self.astep(array)
        apoint = undefined
        stats = undefined
        self.astep = <bound method NUTS.astep of <pymc3.step_methods.hmc.nuts.NUTS object>>
        array = array([0., 0.])
    223             point = self._logp_dlogp_func.array_to_full_dict(apoint)
    224             return point, stats
    225         else:
    226             apoint = self.astep(array)

...........................................................................
/Users/ajoseph/anaconda2/lib/python2.7/site-packages/pymc3/step_methods/hmc/base_hmc.py in astep(self=<pymc3.step_methods.hmc.nuts.NUTS object>, q0=array([0., 0.]))
    112         start = self.integrator.compute_state(q0, p0)
    113 
    114         if not np.isfinite(start.energy):
    115             self.potential.raise_ok()
    116             raise ValueError('Bad initial energy: %s. The model '
--> 117                              'might be misspecified.' % start.energy)
        start.energy = nan
    118 
    119         adapt_step = self.tune and self.adapt_step_size
    120         step_size = self.step_adapt.current(adapt_step)
    121         self.step_size = step_size

ValueError: Bad initial energy: nan. The model might be misspecified.

I also printed the logp values. The emission_mix logp values is nan. (see below)

for RV in simple_dag.basic_RVs:
    print RV.name, ':', RV.logp(simple_dag.test_point)

gives,

state_stickbreaking__ : -2.488704650930508
label : nan

[fonnesbeck]

Just let it do the default initialization. So, don’t manually assign a step method. Also, you don’t need 10k iterations. Try just:

sample(1000, tune=2000)

[antjoseme]

@fonnesbeck It still doesn't work. I have pushed a notebook
https://github.com/antjoseme/explore-DAG/blob/master/DAG-simple-multinomial-mixture.ipynb
that has the whole example .

However, since we last chatted I was able to make it work using categorical distribution instead of a multinomial distribution. Note, I was using multinomial with n=1, so I could substitute it with a categorical distribution.
https://github.com/antjoseme/explore-DAG/blob/master/DAG-simple-categorical-mixture.ipynb

Any idea why one works and the other does not?

[gh-owestesson]

@antjoseme did you ever resolve this? I am working on something similar and any insight would be greatly valued. The switch from Multinomial to Categorical is interesting, but the notebooks linked are no longer there. Could you repost them?

[antjoseme]

Hi @breadbaron the notebooks should be available now

[gh-owestesson]

@antjoseme thanks!

[ricardoV94]

This was accidentally fixed by #3059

It had to do with with the last row of the state_label_tran which includes a zero-probability for the first class.

import numpy as np
import pymc3 as pm

state_label_tran = np.array([
    [0.3, 0.2, 0.5],
    [0.1, 0.5, 0.4],
    [0.0, 0.05, 0.95],  # Would work fine without this row
])
shape = len(state_label_tran)

with pm.Model() as m:
    label_dist = [
        pm.Multinomial.dist(p=state_label_tran[i], n=1, shape=shape)
        for i in range(len(state_label_tran))
    ]
    label = pm.Mixture('label', w=np.ones(shape)/shape, comp_dists=label_dist, observed=[0, 1, 0])

print(label.logp())  # nan before PR #3059, and -1.386... afterwards

In fact it had nothing to do with the mixture, but with the Multinomial:

with pm.Model() as m:
    label = pm.Multinomial('m', n=1, p=[0, 0.05, 0.95], observed=[0, 1, 0])
print(label.logp())  # nan before, and -2.9957 afterwards

None of our current Multinomial tests check for a probability vector that includes zeros (except for test_batch_dirichlet_multinomial which sometimes does so randomly). Maybe it would be good to add one?

[brandonwillard]

None of our current Multinomial tests check for a probability vector that includes zeros (except for test_batch_dirichlet_multinomial which sometimes does so randomly). Maybe it would be good to add one?

Definitely

[ricardoV94]

Closing in favor of #4487