Index error in hierarchical model

Hi,

I’m trying to use Pymc3 to fit the following hierarchical model below:

with pm.Model() as model:
    
    # Hyperpriors
    mu_a = pm.Normal('mu_alpha', mu=17, sd=10)
    sigma_a = pm.HalfCauchy('sigma_alpha', beta=1)
    mu_b = pm.Normal('mu_beta', mu=0, sd=5)
    sigma_b = pm.HalfCauchy('sigma_beta', beta=1)
    
    # Intercept for each category, distributed around group mean mu_a
    a = pm.Normal('alpha', mu=mu_a, sd=sigma_a, shape=2)
    # Slope for each category, distributed around group mean mu_a
    b = pm.Normal('beta', mu=mu_b, sd=sigma_b, shape=2)
    
    # Model error
    eps = pm.HalfCauchy('eps', beta=1)
    
    # Expected value
    NbrCar_est = a[b_liens] + b[b_liens] * b_x
    
    # Data likelihood
    y_like = pm.Normal('y_like', mu=NbrCar_est, sd=eps, observed=b_y) 

Since my data has only two categories I specified shape=2 but I get the following index error:

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-57-aa462638de42> in <module>()
     32 
     33     # Expected value
---> 34     NbrCar_est = a[b_liens] + b[b_liens] * b_x
     35 
     36     # Data likelihood

C:\ProgramData\Anaconda3\lib\site-packages\theano\tensor\var.py in __getitem__(self, args)
    570                 return self.take(args[axis], axis)
    571             else:
--> 572                 return theano.tensor.subtensor.advanced_subtensor(self, *args)
    573         else:
    574             if np.newaxis in args:

C:\ProgramData\Anaconda3\lib\site-packages\theano\gof\op.py in __call__(self, *inputs, **kwargs)
    672                 thunk.outputs = [storage_map[v] for v in node.outputs]
    673 
--> 674                 required = thunk()
    675                 assert not required  # We provided all inputs
    676 

C:\ProgramData\Anaconda3\lib\site-packages\theano\gof\op.py in rval(p, i, o, n)
    890             # default arguments are stored in the closure of `rval`
    891             def rval(p=p, i=node_input_storage, o=node_output_storage, n=node):
--> 892                 r = p(n, [x[0] for x in i], o)
    893                 for o in node.outputs:
    894                     compute_map[o][0] = True

C:\ProgramData\Anaconda3\lib\site-packages\theano\tensor\subtensor.py in perform(self, node, inputs, out_)
   2188         out, = out_
   2189         check_advanced_indexing_dimensions(inputs[0], inputs[1:])
-> 2190         rval = inputs[0].__getitem__(inputs[1:])
   2191         # When there are no arrays, we are not actually doing advanced
   2192         # indexing, so __getitem__ will not return a copy.

IndexError: index 3 is out of bounds for axis 1 with size 2

I can’t figure out where my mistake is.

Colin

Are you sure b_liens only contains 0 and 1?