Gaussian Mixture of regression

The correct way to do in this case is likelihood = pm.Mixture('likelihood', π, pm.Normal.dist(mu, sd), observed=Y)