Fit many regressions simultaneously

Questions
3

Got it, thanks. I had naively assumed that doing this would cause the model to assume it’s a multivariate normal (which is a different class in pymc3…).

The following seems to work:

import numpy as np
import pymc3 as pm

J = 13
N = 500
mu_true = 10
beta_true = 1.5/N
sigma_true = 1.5

X = np.ones(shape=(N,J)) * np.arange(0, N)[:,np.newaxis]
y_obs = np.ones(shape=(N,J))
for j in np.arange(J):
    y_obs[:,j] = mu_true + beta_true * X[:,j] + np.random.normal(loc=0, scale=sigma_true, size=N)
y_obs

with pm.Model() as basic_model:
    # PRIORS
    mu = pm.Normal('mu', 10, 10, shape=J)
    beta = pm.Normal('beta', 0, 2, shape=J)
    sigma = pm.HalfNormal('sigma', 2, shape=J)
    # FOR CONCISION
    mu = mu + beta*X
    # DATA MODEL
    y = pm.Normal('y', mu=mu, sd=sigma, observed=y_obs)

with basic_model:
    trace = pm.sample(draws=1000)

Is there a way to speed up by changing the definition of mu above?

Thanks! (This isn’t Stack Overflow so hope it’s OK to say that…)