Welcome, @DaveD ! Good idea to follow examples: they are useful, but as you encountered, of limited scope.
The “shape” of variables is a typical complication in PyMC, but once you get a hang of it, there are what gives the creative options for model creation.
I tried to recover my own examples for it, and one is here (but note that this was an older version of PyMC).
Some pointers:
I tend to use the “dimensions”/“coordinates” (what you named “hd_vars”?) into the second dimension. For example, try
α = pm.Normal(‘α’, mu = 0.0, sigma = 10.0, shape=(1, len(hd_vars)))
Then, I tend to avoid the * operator, since it is so overloaded: try
aesara.tensor.dot(data[idx], slope)
so the outcome for each summand of your estimator will be shape: (data.shape[0], len(hd_vars)).
Thirdly, on some model components (intercepts, residual) you will have to use a column vector of ones to make them match that shape.
fid = pymc.Data('fake_intercept_data', numpy.ones((data.shape[0],1)), mutable = True)
intercept = aesara.tensor.dot(fid, pymc.Normal('intercept', ..., shape = (1, len(vars))))
Hope this helps!
Cheers