basic_model = pm.Model()
with basic_model:

    # Priors for unknown model parameters
    a = pm.Normal('a', mu=10, sd=10, shape=4)
    b = pm.Normal('b', mu=10, sd=10, shape=4)
    c = pm.Normal('c', mu=0, sd=10)
    sigma = pm.HalfNormal('sigma', sd=1)

    alpha = np.tile(a, (4, 1)).T
    beta = np.tile(b, (4, 1))
    offset = np.tile(c, (4, 4))

    # Expected values of outcome
    mu = sigmoid(alpha, beta, offset)  # 4x4 matrix here

    # Likelihood (sampling distribution) of observations
    # (data also 4x4 matrix)
    Y_obs = pm.Normal('Y_obs', mu=mu, sd=sigma, observed=data)  # ERROR!

with basic_model:
    # draw 500 posterior samples
    trace = pm.sample(500)

alpha = pm.Normal('alpha', mu=0, sd=10)
beta = pm.Normal('beta', mu=0, sd=10, shape=2)
sigma = pm.HalfNormal('sigma', sd=1)
mu = alpha + beta[0] * X1 + beta[1] * X2  # mu = Elemwise{add,no_inplace}.0