How to infer components of mixture?

Questions
1

Hello,
I’m new in pymc3 and I try to infer gaussian mixture components (weight and modes) with ranking data. I have used the following ressource: Order statistics in PyMC3 to do that. The inference is well, the bernouilli is updated but my mixture component doesn’t update…
Here is my code:

with pm.Model() as mixture_model:
        theta = pm.Bernoulli('theta', p = 0.5, shape=K)
        theta_stacked = tt.stack([1-theta, theta], axis=1)
        component0  = pm.Normal('component0', 30, 6, shape=K).distribution
        component1 =  pm.Normal('component1', 10, 4, shape=K).distribution
        mu_hat = pm.Mixture('mu_hat', w=theta_stacked , comp_dists=[component0,component1] , shape=K)
        mu = tt.concatenate([mu_hat])
        latent = pm.Normal('latent',
                           mu=mu[y_argsort],
                           sd=2, 
                           transform=Ordered2D(), 
                           shape=(K,J),
                           testval=np.repeat(np.arange(K)[:,None], J, axis=1))

Anyone to help me please ?

2

Hello, I Have changed my code because it wasn’t really what I wanted to infer…
Here is the new model, let’s imagine that I’m trying to compare mixture rank. To do that I’m using the following ressource Order statistics in PyMC3 . I got two components for each mixture in my case, and here is the new code that I’m using, K is the number of Mixture that I want to compare between each other:

y_argsort = np.argsort(y, axis=0)
component_1_dic = {}
component_2_dic = {}
mu_dic = {}
weight_dic = {}
with pm.Model() as mixture_model:
  for k in range(K):
    component_1_dic['component_1' + str(k)]  = pm.Normal('component_1' + str(k), mu = 30, sigma = 8,  shape=1)
    component_2_dic['component_2' + str(k)]  = pm.Normal('component_2' + str(k), mu = 15, sigma = 6,  shape=1)
    weight_dic['weight_' + str(k)] = pm.Bernoulli('weight_'+ str(k), p = 0.5, shape=1)
    theta_stacked = tt.stack([1-weight_dic['weight_' + str(k)], weight_dic['weight_' + str(k)]], axis=1)
    component_stacked = [component_1_dic['component_1' + str(k)].distribution, component_2_dic['component_2' + str(k)].distribution]
    mu_dic['mixture_' + str(k)] = pm.NormalMixture('mixture_' + str(k), w = theta_stacked , mu  = tt.concatenate([component_1_dic['component_1' + str(k)], component_2_dic['component_2' + str(k)]]), sigma = pm.floatX([3, 4]), shape = 1)
  mu = tt.concatenate([mu_dic['mixture_' + str(0)]])
  for k in range(1, K):
    mu = tt.concatenate([mu, mu_dic['mixture_' + str(k)]])
  latent = pm.Normal('difference',
                     mu=mu[y_argsort],
                     sd=1, 
                     transform=Ordered2D(), 
                     shape=(K,J),
                     testval=np.repeat(np.arange(K)[:,None], J, axis=1))
  trace = pm.sample(100000, tune=10000)

My point is already the same, the components of my mixture doesn’t change a lot and python return each time : The number of effective samples is smaller than 10% for some parameters. I have tryied to rise the number of sample to the million (I have acces to a very good computer), and is already the same… Anyone could help me please ?

3

I am still not very clear about your use case - are you assuming that the latent rank follows a mixture?

4

Thanks for your answers. So let’s imagine the following model:

B_i \sim \mathcal{B}(p_i) \\ S_i \sim p_i\ \mathcal{N}(\mu_{i1}, \sigma_{i1} + c_1) + (1-p_i)\mathcal{N}(\mu_{i2}, \sigma_{i2} + c_2) \\ P_i \sim \mathcal{N}(S_i, 1) \\

And my observation is P_1 < P_2 <..<P_K (a rank). I want to infer p_i, \mu_{i1}, \ \mu_{i2}, \sigma_{i1}, \sigma_{i2} for i = 1,..K. c_1 and c_2 are know, in my code I set c_1 = 3, c_2 = 4.
I hope that the code reflect the model. I Have tried a lot of differents number of samples, but I always got the same message:
The number of effective samples is smaller than 10% for some parameters.

5

In fact, because I’m only interested by some parameters, a solution could be to specifity in differents steps the variables on which inference should be focus in first with something like that:

  list_interest_continous = []
  list_interest_discrete = []
  for k in range(K):
    list_interest_continous.append(component_1_dic['component_1' + str(k)])
    list_interest_continous.append(component_2_dic['component_2' + str(k)])
    list_interest_discrete.append(weight_dic['weight_' + str(k)] )

  step1 = pm.NUTS( vars = list_interest_continous)
  step2 = pm.BinaryGibbsMetropolis(vars=list_interest_discrete)
  trace = pm.sample(100000,[step1, step2])
6

Where is B_i came in in your model? are they observed?
I still not sure if I understand your intent clearly, as from my experience the rank model is already some kind of latent mixture, and mixture of mixture is pretty unidentifiable - maybe it is easier to demonstrate if you could write down your data simulation process?

7

The B_i are not observed. I observe only rank, and I assume that rank are generated indirectly by mixture. For moment my data (as in Order statistics in PyMC3) is :

K = 4 # number of items being ranked
J = 1 # number of raters
yreal = np.argsort(np.random.randn(K, 1), axis=0)
print(yreal)
y = np.argsort(yreal + 2*np.random.randn(K, J), axis=0)
print(y)

But in fact I want to be more fleaxible with the data simulation process, that’s why I’ve choosen latent Mixture instead of latent Normal.