Try pm.sample(init='adapt_diag'), and also providing a test_value to mu_b: mu_b = Normal('mu_b', mu=0., sd=0.1, shape=k, testval=np.asarray([-.1,.1]))
Those two notebooks are mixture regressions, sorry there is not many comments but you can see different ways of modelling it.