Mixing/sampling coefficients from multiple linear regressions

Questions v5
1

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(
        "obs",
        mu=func,
        sigma=err,
        observed=data,
    )

    # 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
    pm.sample()

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