A couple little things:
-
For the standard deviations, use
pm.HalfNormalin place ofpm.Normal -
You shouldn’t need an extra level for error; the variances are already part of the GMM
-
Right now, the slopes and intercepts are being treated as a mixture; instead of the y-values
-
The cluster weights are not being treated hierarchically
For sharing weights I might do something like:
def logp_alpha(a, l=0.5):
return tt.log(tt.exp(-l * a))
with pm.Model() as linear_gmm:
# prior on mixtures: try to be uninformative, but penalize very large values of ||a||
alpha_prior = pm.Flat('alpha_prior', shape=n_clusters)
pm.Potential('p(alpha)', logp_alpha(alpha))