Hierarchical Gaussian Mixture Regression Model

A couple little things:

  1. For the standard deviations, use pm.HalfNormal in place of pm.Normal

  2. You shouldn’t need an extra level for error; the variances are already part of the GMM

  3. Right now, the slopes and intercepts are being treated as a mixture; instead of the y-values

  4. 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))