I’m coming in late to the thread, and have only really read the most recent post. I agree with you, Kyle, that power analysis is perfectly well-motivated in a Bayesian setting. One could re-cast it as a "how many samples do I need for the posterior variance of the estimate of the difference in rates/means to be smaller than \theta" – and then note that for a fixed difference \Delta_\mu this is equivalent to a statement about p-values.
As far as speeding up the model – even though the models are independent, using NUTS to sample from a (1000, 1000) space is (I suspect) slower than sampling from a (100, 100) space 100 times.
In this specific case where you have used a conjugate distribution (beta/binomial); the posterior has a closed-form solution. I’d recommend using pm.Deterministic to compute the parameters and pass them into a posterior = pm.Beta(...); that way you can just use pm.sample_prior_predictive and forgo the use of NUTS entirely.