How can we build a mixture of mixtures?

The Mixture class is designed for creating a mixture of 1D distributions, more complex use cases like multi-dimensional mixture or mixture of mixtures is not tested and not properly supported. There is no unit test for these usages.

I dont think using a custom logp would create that much difficulty for debugging, for example, using a DensityDist:

def mix_mixlogp(w, comp_dists):
    def logp_(value):
        comp_logp = tt.squeeze(tt.stack([comp_dist.logp(value) 
                                         for comp_dist in comp_dists], axis=1))
        return pm.math.logsumexp(tt.log(w) + comp_logp, axis=-1)
    return logp_

with pm.Model() as model:
    ...
    mix = pm.DensityDist('mix', mix_mixlogp(mix_w, [g_mix, l_mix]), observed=data)

or using a Potential

def mix_mixlogp(w, comp_dists, value):
    comp_logp = tt.squeeze(tt.stack([comp_dist.logp(value) 
                                     for comp_dist in comp_dists], axis=1))
    return pm.math.logsumexp(tt.log(w) + comp_logp, axis=-1)

with pm.Model() as model:
    ...
    mix = pm.Potential('mix', mix_mixlogp(mix_w, [g_mix, l_mix], data))

Both case only the mixture of mixture is using a custom density.