Interestingly, rerunning the last cell in the gist with sample(random_seed=1), it appears that sometimes the two chains disagree on which Gaussian is which:
>>> values = trace.get_values('mu', combine=False)
>>> values[0][-1], values[1][-1]
(array([-1.05009326, 2.59436901]), array([ 2.59135272, -1.01800517]))
This results in r_hat values on each of the n_shape=2 parameters (mu, w, sd) of 1.83 and warnings. Is there anything that can be done to avoid that?