This might have to do with the mu/sd parametrization.
Your mu is on the correct interval, but if I see it correctly, sigma_a might occasionally get out of the interval 0 < sd < sqrt(mu * (1 - mu))
You can try two things:
- use a bounded sigma
- check if the model runs with alpha/beta parametrization (optionally use
pm.Deterministic()to calculate \mu on the way)
The second option should be more robust.