Hello folks,

I am currently working on a modification of the model from the “Gaussian Mixture Model” example notebook (link).

The original model looks like this.

```
with pm.Model(coords={"cluster": range(k)}) as model:
μ = pm.Normal(
"μ",
mu=0,
sigma=5,
transform=pm.distributions.transforms.ordered,
initval=[-4, 0, 4],
dims="cluster",
)
σ = pm.HalfNormal("σ", sigma=1, dims="cluster")
weights = pm.Dirichlet("w", np.ones(k), dims="cluster")
pm.NormalMixture("x", w=weights, mu=μ, sigma=σ, observed=x)
```

I am now trying to get another dimension into the model and sample multiple mixtures at once.

Like this:

```
with pm.Model(coords={"cluster": range(k),"set":range(2)}) as model:
μ = pm.Normal(
"μ",
mu=0,
sigma=5,
shape=(2,3),
transform=pm.distributions.transforms.ordered,
initval=np.array([[-4, 0, 4],[-4,0,4]]),
dims=("set","cluster"),
)
σ = pm.HalfNormal("σ", sigma=1, dims=("set","cluster"))
weights = pm.Dirichlet("w", np.ones(k), dims=("set","cluster"))
pm.NormalMixture("x", w=weights, mu=μ, sigma=σ, observed=np.stack([x,x],axis=-1))
```

However, with the ordered transformation I am repeatedly running into shape issues, not if I use no transformation. Is this a bug? Does anybody know a way how I can make it work?

Thank you very much in advance