Dirichlet process mixtures bidimensional data

Questions
1

Hi folks, I’m studding Dirichlet mixtures process to evaluate the number of cluster of a bidimensional dataset.
In the tutorial example you asses the the number of clusters in a one dimensional data, whats wrong with this use to bidimensional data?

with pm.Model() as model:
    alpha = pm.Gamma('alpha', 1., 1.)
    beta = pm.Beta('beta', 1., alpha, shape=(K))
    w = pm.Deterministic('w', stick_breaking(beta))

    tau = pm.Gamma('tau', 1., 1., shape=(K,2))
    lambda_ = pm.Uniform('lambda', 0, 5, shape=(K,2))
    mu = pm.Normal('mu', 0, tau=lambda_ * tau, shape=(K,2))
    obs = pm.NormalMixture('obs', w, mu, tau=lambda_ * tau, comp_shape=(K,2),
                           observed=DATA[['x', 'y']].values)

Thank you all

2

what kind of error are you seeing?

3

dimension mismatch. Already tried (K,2) without success…

4

I see, I think you need to supply the comp_shape to NormalMixture.

5

It was supposed to be like the edit above?

6

Yes, and I take it doesnt work? what is the error you are seeing?

7

For k=35:

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

8

Workaround:

def stick_breaking(v):
    return v * tt.concatenate([tt.ones_like(v[..., :1]),
                               tt.extra_ops.cumprod(1 - v, axis=1)[..., :-1]],
                              axis=1)

K = 10
with pm.Model() as model:
    alpha = pm.Gamma('alpha', 1., 1.)
    beta = pm.Beta('beta', 1., alpha, shape=(1, K))
    w = pm.Deterministic('w', stick_breaking(beta))

    tau = pm.Gamma('tau', 1., 1., shape=(2, K))
    lambda_ = pm.Uniform('lambda', 0, 5, shape=(2, K))
    mu = pm.Normal('mu', 0, tau=lambda_ * tau, shape=(2, K))
    obs = pm.NormalMixture('obs', w, mu, tau=lambda_ * tau, comp_shape=(2, K),
                           observed=np.random.randn(100, 2, 1))

Explanation:
1, component of mixture should be the last dimension, which gets collapse.
2, the observation need to be evaluated on each component
–> solution: a place holder at the last dimension for observed, thus allow it to be evaluated on each component.

9

Got it! Thank you very much! You never disappoint!

10

I’m seeing this in the sampler run:

ValueError: No free random variables to sample.

or

ValueError: The model does not contain any free variables.

if I let pymc3 choose the sampler.

11

As shown in the error, you didnt defined any random variable. Usually that happen when you fixed all the random variables to a constant.

12

But as you can see, there are a bunch of hyperpriors (alpha, beta, tau, lambda_, mu). I’m running exactly as you wrote, because I’m trying to understand how it works before any further setting.

13 
def stick_breaking(v):
    return v * tt.concatenate([tt.ones_like(v[..., :1]),
                               tt.extra_ops.cumprod(1 - v, axis=1)[..., :-1]],
                              axis=1)

K = 10
with pm.Model() as model:
    alpha = pm.Gamma('alpha', 1., 1.)
    beta = pm.Beta('beta', 1., alpha, shape=(1, K))
    w = pm.Deterministic('w', stick_breaking(beta))
    w = w/w.sum()

    tau = pm.Gamma('tau', 1., 1., shape=(2, K))
    lambda_ = pm.Uniform('lambda', .0001, 5, shape=(2, K))
    mu = pm.Normal('mu', 0, tau=lambda_ * tau, shape=(2, K))
    obs = pm.NormalMixture('obs', w, mu, tau=lambda_ * tau, comp_shape=(2, K),
                           observed=np.random.randn(100, 2, 1))
    trace = pm.sample()

Could you try the above code?

14

Worked! Weren’t the weights already summing 1 in the stick_breaking function?

15

Nope, but when you have a large K usually that’s not a problem.

16

Thank you once more my friend!

17

Just a note to say, my problem wasn’t multidimensional observations, only that the observations occurred at a higher level so I couldn’t modify the observed data dimensions. But, adding a shape argument with an empty extra dimension seemed to do the trick.

In this example, it would have been:

obs = pm.NormalMixture('obs', w, mu, tau=lambda_ * tau, 
    comp_shape=(2, K), shape=(100, 2, ))