How can we build a mixture of mixtures?

Questions
21

Thanks a lot !! The problem is solved on the master…

But now (with the master) a problem occurs on the mixture of mixture…

import sys
import pymc3 as pm, theano.tensor as tt
import numpy as np

import matplotlib.pyplot as plt

# simulated data 
data = np.concatenate( (
    np.random.normal(loc=1,scale=1,size=1000),
    np.random.lognormal(mean=1,sigma=1,size=1000)) )


with pm.Model() as model:
    nbr = 5
    # mixtures components
    g_components = pm.Normal.dist(mu=pm.Exponential('mu_g', lam=1.0,shape=nbr)
    , sd=1,shape=nbr) #pm.HalfNormal('sd_g', sd=1.0,shape=nbr)
    l_components = pm.Lognormal.dist(mu=pm.Exponential('mu_l', lam=1.0,shape=nbr)
    , sd=1,shape=nbr) #pm.HalfNormal('sd_l', sd=1.0,shape=nbr)
    # weight vector for the mixtures
    g_w = pm.Dirichlet('g_w',a=np.array([0.0000001]*nbr))
    l_w = pm.Dirichlet('l_w',a=np.array([0.0000001]*nbr))
    mix_w = pm.Dirichlet('mix_w',a=np.array([1]*2))
    # mixtures
    g_mix = pm.Mixture.dist(w=g_w,comp_dists=g_components)
    l_mix = pm.Mixture.dist(w=l_w,comp_dists=l_components)
    # mixture of mixture
    mix = pm.Mixture('mix',w=mix_w,comp_dists=[g_mix,l_mix], observed=data)
    # MCMC
    trace = pm.sample(1000, tune=1000, live_plot=True)

    # graphes
    pm.traceplot(trace)
    plt.show()

It says (without any link to my code.) :

AttributeError: 'list' object has no attribute 'logp'

And another error :

ValueError: Input dimension mis-match. (input[0].shape[1] = 2, input[1].shape[1] = 5)

This error occurs on the mixture of mixture :

mix = pm.Mixture.dist('mix',w=mix_w,comp_dists=[g_mix,l_mix], observed=data)

None of these two errors occurs when using the pymc3.3 version.

How can I have a version that is both able to handle mixture of mixtures, and handle correctly when given an array to the logp function (The preivous question you solved)?

22

Mixture of Mixture is really stretching the intended use case of Mixture class. Try cherry picking the mixture logp function and write it explicitly. For example, see cell 2 in http://docs.pymc.io/notebooks/gaussian-mixture-model-advi.html

23

You said : “Mixture of Mixture is really stretching the intended use case of Mixture class”.

But it worked on version 3.3 ! (not the logp function of a mixture, but the rest was ok!)
So the master code looks like a regression on that point.

After all, a mixture is also a distribution, so there is no theoretical reason not to be able to build mixture of mixtures ?
It would be very sad to abandon this marvelous feature… It helps to build cleaner & more structured models.

Writing explicitely the logp function doesn’t seem like a good solution to me, because this can be considered like “duplicate code” : I’ll have to duplicate logp fomulas from the distribution components… (errors will be very hard to track)
(Thanks to the link to the code example you provided, I found it very interesting, but I think I will use another solution.)

In my case I’ll end up putting all components in one unique mixture, and
duplicate the model part where I need to access the logp function. It will be only model description duplicate, not inner code with math formulas. So I expect it will be less risky and easier to debug, but I admit that this is not totally satisfying…

So It would be reallly cool if the mixture of mixture feature work on the next pymc3 release.

Thanks for all your answers and suggestions.

24

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.

25

OK, I submitted a PR to make working with mixture of mixtures easier:

Using your mixture of mixtures example as a test case :wink:
What is your Github account btw? I would like to credit you as well in the PR.

26

Great !! Thanks you very much for you help and implication.

How could I have access to your modified version ?

(my Github account is : hwassner)

27

OK I just push the PR - if you update pymc3 to master your code should run with no problem now :wink:

1 Like
28

2 posts were split to a new topic: Mixture logp related

29

@junpenglao What’s the recommended approach for a mixture of multivariate Gaussians in 2021?