Debugging pymc probability calculations

Question

I've tried to model a mixture of exponentials by copying the mixture-of-Gaussians example given here. The code is below. I know there are some funky aspects to the inference here, but my question is more about how to debug the calculations in models like this.

The idea is that it's a mixture of three exponentials, with scale parameters taken from the Gamma assigned to scales. However, all observations get assigned to the zeroth exponential during the ElemwiseCategoricalStep. You can see that the assignments of the observations to the exponential components are initially diverse by looking at initial_assignments, and you can see that all observations are assigned to the zeroth component on all interations from the fact that set(tr['exp'].flatten()) contains only 0.

I assume this is because all of the values assigned to p in the expression array([logp(v * self.sh) for v in self.values]) in ElemwiseCategoricalStep.astep are minus infinity. I would like to know why that is and how to correct it, but even more, I would like to know what tools are available to debug this kind of thing. Is there any way for me to step through the calculation of logp(v * self.sh) to see how the result is determined? If I try to do it using pdb, I think I get stymied at outputs = self.fn() in theano.compile.function_module.Function.__call__, which I guess I can't step into because it's a native function.

Even knowing how to compute the pdf for a given set of model parameters would be a useful start.

import numpy as np
import pymc as pm
from pymc import Model, Gamma, Normal, Dirichlet, Exponential
from pymc import Categorical
from pymc import sample, Metropolis, ElemwiseCategoricalStep

durations = np.concatenate(
    [np.random.exponential(1/lam, 10)
     for lam in [1e-3,7e-5,2e-6]])

initial_assignments = np.random.randint(0, 3, len(durations))

print 'initial_assignments', initial_assignments

with Model() as model:
    scales = Gamma('hp', 1, 1, shape=3)
    props = Dirichlet('props', a=np.array([1., 1., 1.]), shape=3)
    category = Categorical('exp', p=props, shape=len(durations))
    points = Exponential('obs', lam=scales[category], observed=durations)
    step1 = pm.Metropolis(vars=[props,scales])
    step2 = ElemwiseCategoricalStep(var=category, values=[0,1,2])
    start = {'exp': initial_assignments,
             'hp':  np.ones(3),
             'props': np.ones(3),}
    tr = sample(3000, step=[step1, step2], start=start)

print set(tr['exp'].flatten())

Answer 1

Excellent question. One thing you can do is look at the pdf for each of the components.

The Model and each of the variables should have both a .logp and a .elemwise_logp property and them which returns a function that can take a point or parameter values.

Thus you can say something like print scales.logp(start) or print model.logp(start) or print scales.dlogp()(start).

For now, I think you unfortunately have to specify all the parameter values (even ones that don't affect the result for a particular variable).

Model, FreeRV and ObservedRV all inherit from Factor which provides this functionality and has a few other methods. You'll probably want the non fast versions since those are more forgiving in the kinds of arguments they accept.

Does that help? Please let me know if you have other ideas for things that might help you in debugging. This is one area where we know pymc3 and theano needs some work.