In [7]:
from pymc3 import Categorical, Model, DensityDist
from pymc3.distributions.dist_math import bound
from theano.tensor import switch, log, le, constant
from theano import shared, gradient
import numpy as np

with Model() as simple_model:
    
    Age = Categorical('Age', p=[0.3, 0.5], testval=1)
    
    Sex = Categorical('Sex', p=[0.6, 0.4], testval=0)
    
    def edu_logp(value):
        p = constant(np.array([[0.25, 0.28, 0.12], [0.36, 0.30, 0.10]]))
        return bound(switch(value, log(p[Sex, Age]), log(1 - p[Sex, Age])),
            value >= 0, value <= 1)

    Education = DensityDist('Education', edu_logp, dtype='int64', testval=0)
    
    def res_logp(value):
        p = constant(np.array([0.75, 0.8]))
        return bound(switch(value, log(p[Education]), log(1 - p[Education])),
            value >= 0, value <= 1)

    Residence = DensityDist('Residence', res_logp, dtype='int64', testval=0)
    
    
    def occ_logp(value):
        p = constant(np.array([0.04, 0.08]))
        return bound(switch(value, log(p[Education]), log(1 - p[Education])),
            value >= 0, value <= 1)

    Occupation = DensityDist('Occupation', occ_logp, dtype='int64', testval=0)
    
    
    def travel_logp(value):
        p =  constant(np.array([[[0.48, 0.42, 0.10],  # Small and employer
                              [0.56, 0.36, 0.08]], # Small and self-employed
                             [[0.58, 0.24, 0.18],  # Big and employer
                              [0.70, 0.21, 0.09]]])) # Big and self-employed

        return bound(log(p[Residence, Occupation, value]),
            value >= 0,
            value <= 2)

    Travel = DensityDist('Travel', travel_logp, dtype='int64', testval=0)
In [8]:
from pymc3 import find_MAP, sample

with simple_model:
    tr = sample(1000)
Assigned BinaryMetropolis to Age
Assigned BinaryMetropolis to Sex
Assigned Metropolis to Education
Assigned Metropolis to Residence
Assigned Metropolis to Occupation
Assigned Metropolis to Travel
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
/Users/fonnescj/anaconda3/lib/python3.5/site-packages/theano/compile/function_module.py in __call__(self, *args, **kwargs)
    858         try:
--> 859             outputs = self.fn()
    860         except Exception:

IndexError: index out of bounds

During handling of the above exception, another exception occurred:

IndexError                                Traceback (most recent call last)
<ipython-input-8-c293546f01d7> in <module>()
      2 
      3 with simple_model:
----> 4     tr = sample(1000)

/Users/fonnescj/Github/pymc3/pymc3/sampling.py in sample(draws, step, start, trace, chain, njobs, tune, progressbar, model, random_seed)
    153         sample_args = [draws, step, start, trace, chain,
    154                        tune, progressbar, model, random_seed]
--> 155     return sample_func(*sample_args)
    156 
    157 

/Users/fonnescj/Github/pymc3/pymc3/sampling.py in _sample(draws, step, start, trace, chain, tune, progressbar, model, random_seed)
    162     progress = progress_bar(draws)
    163     try:
--> 164         for i, strace in enumerate(sampling):
    165             if progressbar:
    166                 progress.update(i)

/Users/fonnescj/Github/pymc3/pymc3/sampling.py in _iter_sample(draws, step, start, trace, chain, tune, model, random_seed)
    244         if i == tune:
    245             step = stop_tuning(step)
--> 246         point = step.step(point)
    247         strace.record(point)
    248         yield strace

/Users/fonnescj/Github/pymc3/pymc3/step_methods/compound.py in step(self, point)
     12     def step(self, point):
     13         for method in self.methods:
---> 14             point = method.step(point)
     15         return point

/Users/fonnescj/Github/pymc3/pymc3/step_methods/compound.py in step(self, point)
     12     def step(self, point):
     13         for method in self.methods:
---> 14             point = method.step(point)
     15         return point

/Users/fonnescj/Github/pymc3/pymc3/step_methods/arraystep.py in step(self, point)
    116         bij = DictToArrayBijection(self.ordering, point)
    117 
--> 118         apoint = self.astep(bij.map(point))
    119         return bij.rmap(apoint)
    120 

/Users/fonnescj/Github/pymc3/pymc3/step_methods/metropolis.py in astep(self, q0)
    123 
    124 
--> 125         q_new = metrop_select(self.delta_logp(q,q0), q, q0)
    126 
    127         if q_new is q:

/Users/fonnescj/anaconda3/lib/python3.5/site-packages/theano/compile/function_module.py in __call__(self, *args, **kwargs)
    869                     node=self.fn.nodes[self.fn.position_of_error],
    870                     thunk=thunk,
--> 871                     storage_map=getattr(self.fn, 'storage_map', None))
    872             else:
    873                 # old-style linkers raise their own exceptions

/Users/fonnescj/anaconda3/lib/python3.5/site-packages/theano/gof/link.py in raise_with_op(node, thunk, exc_info, storage_map)
    312         # extra long error message in that case.
    313         pass
--> 314     reraise(exc_type, exc_value, exc_trace)
    315 
    316 

/Users/fonnescj/anaconda3/lib/python3.5/site-packages/six.py in reraise(tp, value, tb)
    683             value = tp()
    684         if value.__traceback__ is not tb:
--> 685             raise value.with_traceback(tb)
    686         raise value
    687 

/Users/fonnescj/anaconda3/lib/python3.5/site-packages/theano/compile/function_module.py in __call__(self, *args, **kwargs)
    857         t0_fn = time.time()
    858         try:
--> 859             outputs = self.fn()
    860         except Exception:
    861             if hasattr(self.fn, 'position_of_error'):

IndexError: index out of bounds
Apply node that caused the error: Subtensor{int64, int64, int64}(TensorConstant{[[[ 0.48  ..1  0.09]]]}, ScalarFromTensor.0, ScalarFromTensor.0, ScalarFromTensor.0)
Toposort index: 9
Inputs types: [TensorType(float64, 3D), Scalar(int64), Scalar(int64), Scalar(int64)]
Inputs shapes: [(2, 2, 3), (), (), ()]
Inputs strides: [(48, 24, 8), (), (), ()]
Inputs values: ['not shown', 2, 0, 1]
Outputs clients: [[Elemwise{Composite{((Switch((GE(i0, i1) * LE(i0, i2)), Switch(i0, log(i3), Composite{log((i0 - i1))}(i4, i3)), i5) + Switch(Composite{(GE(i0, i1) * LE(i0, i2))}(i6, i1, i7), log(i8), i5)) - (Switch((GE(i9, i1) * LE(i9, i2)), Switch(i9, log(i3), Composite{log((i0 - i1))}(i4, i3)), i5) + Switch(Composite{(GE(i0, i1) * LE(i0, i2))}(i6, i1, i7), log(i10), i5)))}}(Subtensor{int64}.0, TensorConstant{0}, TensorConstant{1}, Subtensor{int64}.0, TensorConstant{1.0}, TensorConstant{-inf}, Travel_shared, TensorConstant{2}, Subtensor{int64, int64, int64}.0, Subtensor{int64}.0, Subtensor{int64, int64, int64}.0)]]

Backtrace when the node is created:
  File "<ipython-input-7-21690581ab52>", line 42, in travel_logp
    return bound(log(p[Residence, Occupation, value]),

HINT: Use the Theano flag 'exception_verbosity=high' for a debugprint and storage map footprint of this apply node.