This is a cryptic one, so here is an equally cryptic solution.
Specify the shape of rho:
rho = pm.Beta("rho", 5, 1, dims=("lags",), shape=(1,))
Or more robust:
rho = pm.Beta("rho", 5, 1, dims=("lags",), shape=(lags,))
Coordinates with length 1 can be tricky because PyTensor treats unknown dimensions and dimensions of length 1 very differently, and doesn’t like to interchange between them.
When PyMC compiles the logp+dlogp, it swaps from unknown dimensions to the static shapes (since dims are not allowed to change during sampling), and for very in-the-weeds reasons the scan used in the logp does not like this change. Where it used to see an unknown length variable, it now sees a length 1, boo-oh!
If you had started with a model of two lags, this shouldn’t have raised an error.