So I am not an expert on inner workings of the sampler, but I believe it basically comes down to gradients. Old sampling algorithms evaluate the log probability of any given set of parameters given the data and model. Newer algorithms require not just an evaluation of the log probability, but also of the gradient. That gives the sampler information about which directions in (high dimensional) parameter space lead to what changes in the log probability.
The reason why we can’t just write arbitrarily complex python code is that while it would provide the logp, it wouldn’t provide the gradient information. This is why we shift over to using Aesara (now PyTensor) or JAX or similar libraries… they provide gradient information that the sampler needs.
So getting to your question… I guess the answer is that while the ordered transformation ‘works’, it is perhaps hard to correctly do all the autodifferentiation etc to get the gradients. But by taking a slightly loner route around (code wise) with the Dirichlet and the cumulative sum, I assume that these operations are much more ‘friendly’ in terms of calculating the gradient information.
But I’m tagging @lucianopaz and @ricardoV94 to see if I’ve made any major errors in my explanation.