Hello,

I am trying to fit a Hierarchical Softmax Regression Model, but I am struggling with specifying the beta coefficients.

If I create a list of betas for each participant (option #3) I can get the model to run, but if instead I try to specify the ‘equivalent’ shape inside a pm.Normal (option #1), or to stack a bunch of pm.Normal in a pm.Deterministic (option #2) I cannot get the model to run.

I have the feeling that either specification would correspond to the same model, but I must be wrong. The reason I would prefer options #1 and #2 is that it makes it easier to plot the posterior distributions as well as making the graphviz model tidier (for option #1).

```
with pm.Model() as m_pp: # partial pooling
beta = pm.Normal('beta', 0, 10, shape=(len(weights)))
spread = pm.Exponential('spread', .1, shape=len(weights))
# These do not work
# betas = pm.Normal('betas', beta, spread, shape=(participants, len(weights)))
# betas = pm.Deterministic('betas', tt.stack([pm.Normal(f'beta[{i}]', beta, spread, shape=(len(weights))) for i in range(participants)]))
# This does
betas = [pm.Normal(f'beta[{i}]', beta, spread, shape=(len(weights))) for i in range(participants)]
mu1 = pm.Deterministic('mu1', tt.concatenate([pm.math.dot(x1[i], betas[i]) for i in range(participants)]))
mu2 = pm.Deterministic('mu2', tt.concatenate([pm.math.dot(x2[i], betas[i]) for i in range(participants)]))
mu3 = pm.Deterministic('mu3', tt.concatenate([pm.math.dot(x3[i], betas[i]) for i in range(participants)]))
mus = pm.Deterministic('mus', tt.transpose(tt.stack((mu1, mu2, mu3))))
theta = pm.Deterministic('theta', tt.nnet.softmax(mus))
yl = pm.Multinomial('y', n=1, p=theta, observed=y)
trace_pp = pm.sample(500, target_accept=.85)
pm.model_to_graphviz(m_pp, )
```