Managing indices/data shape in a multi-level linear regression for incremental model updating

Questions
1

I was hoping to do something similar to this https://github.com/pymc-devs/pymc3/blob/master/docs/source/notebooks/updating_priors.ipynb but for the model defined below (toy data is provided and variables have been semi-anonymized). I’m realizing that indexing etc is going to be kind of a headache though. Any suggestions?

import numpy as np
import pymc3 as pm

c_count = 2
coh_count = 58
c_c_indices = np.asarray([0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1,
       0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1,
       0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1])
coh_index = sorted(np.tile(np.arange(0,58), [52,]))
exposure = np.random.randint(10, 200, len(coh_index))
observed = np.random.randint(30, 400, len(coh_index))
x = np.tile(np.log(np.arange(1,53)), [58,])

with pm.Model() as model:
    sigma_1 = pm.Exponential('sigma_1', .1)
    sigma_2 = pm.Exponential('sigma_2', .1)
    sigma_3 = pm.Exponential('sigma_3', .1)
    sigma_4 = pm.Exponential('sigma_4', .1)

    global_m = pm.Normal('global_m',
                          mu=-1,
                          sigma=1)
    global_b = pm.Normal('global_b',
                          mu=-1,
                          sigma=1)
    c_m = pm.Normal('c_m',
                     mu=global_m,
                     sigma=sigma_1,
                     shape=c_count)
    c_b = pm.Normal('c_b',
                     mu=global_b,
                     sigma=sigma_2,
                     shape=c_count)
    c_c_m = pm.Normal('c_c_m',
                       mu=c_m[c_c_indices],
                       sigma=sigma_3,
                       shape=coh_count)
    c_c_b = pm.Normal('c_c_b',
                       mu=c_b[c_c_indices],
                       sigma=sigma_4,
                       shape=coh_count)

    theta = (c_c_b[coh_index] + x * c_c_m[coh_index] + pm.math.log(exposure))
    txn_like = pm.Poisson('txn_like', mu=np.exp(theta), observed=observed)

with model:
    trace = pm.sample(2000)