Your model works for me (though I didn’t test it by any means) if you call pm.sample(init='adapt_diag') instead of pm.sample(). The default is pm.sample(init='jitter+adapt_diag'), but the jittering can sometimes cause problems during initialization (but not sampling proper), e.g., when the scale of parameters is very small.