Pymc3 does not reformulate the model you provide it with under the hood. For example, if you have a string of iid Bernoulli observations, and you’re only interested in find p, pymc3 does not reformulate the model into a binomial, which is more memory friendly and also faster in this case.
Looking at your initial model, I think that you should be able to supply the total counts if you increase each multinomial’s n.