I think I understand what is happening, and it isn’t a GPU/CPU thing. Bambi is clever enough to use pm.invlogit() only for the unique values of beta_0 + mu_0[x], while my code ignored that there are duplicates in x and thus in beta_0 + mu_0[x].
In real case scenarios this is not a problem, because I can aggregate identical rows and use the binomial distribution instead of bernoulli to achieve the same optimisation. This was only a “pedagogical example”, and indeed I learnt something from it 