Mixing/sampling coefficients from multiple linear regressions

Hi all,
I’m currently investigating how bayesian statistics and pymc can be used in an analysis of extreme weather events. In broad strokes the current analysis involves fitting a distribution to some pooled data (m x n) using a bootstrap and performing a linear regression between each member series in the data and a predictor (n), i.e. m linear regressions against the predictor.

From the examples in the documentation (Thanks for these!) I have figured out how to fit a distribution and perform a linear regression in pymc, so far so good. In the case of the linear regression this means I have priors for the slope and intercept with shape=m, so each regression have individual distributions for slope and intercept.

In a later stage of the analysis I want to sample all the slopes from the regressions. As I understand it a mixture with equal weights could be used for this, but I’m not sure if this is the right way to go? Testing this leads to a lot of divergences. My code for the linear regression looks like this

# data is an 50 x 100 array
# predictor is an 50 x 1 array broadcasted to 50 x 100
with pm.Model():
    # Coeffs
    a = pm.Normal("slopes", mu=0, sigma=10, shape=100)
    b = pm.Normal("intercepts", mu=0, sigma=10, shape=100)
    err = pm.HalfCauchy("sigma", beta=10, initval=1.0)

    # Predictor
    x_ = pm.ConstantData("predict", predictor)
    # The linear regression
    func = a * x_ + b
    # Regression
    obs = pm.Normal(

    # Running without this part works with no divergences
    w = pm.Dirichlet('w', a=np.ones(100))
    pooled_slopes = pm.NormalMixture("pooled slopes", w=w, mu=a, sigma=20)

    # Sample the model

Very new to the world of bayesian statistics so grateful for any input!