Dirichletmultinomial mixture model doesn't mix

Questions v5
1

Hello.
I’m trying to model a system where the observable variable is the cumulative change of Nr variables multiplied by Nr factors that depend on time (Nf observations). Each of these variables have 3 possibilities for each timestep: decrease, increase or remain unchanged, each of these options have a prior proportion and a specific distribution on their magnitude. My current model looks like this:

    model = pm.Model()
    with model:
        def negWB(alpha, beta, size):
            return -1.0 * pm.Weibull.dist(alpha=alpha, beta=beta, size=size)

        w = pm.DirichletMultinomial('w', n=1, a=ws_, shape=(Nf, Nr-1, 3))
        #Mixture doesn't like elements on the boundaries
        eps = 1E-6
        ww = (w + eps)/(1.0 + 3.0*eps)
        change_c = pm.Weibull.dist(alpha=alphas_c_, beta=betas_c_, shape=(Nf, Nr-1))
        change_v = pm.CustomDist.dist(alphas_v_, betas_v_, class_name="negWB", dist=negWB)
        neutral = pm.Normal.dist(mu=np.zeros_like(alphas_c_), sigma=1E-6, shape=(Nf, Nr-1))

        components = [change_v, neutral, change_c]
        delta0 = pm.Mixture('delta0',
                           w=ww, 
                           comp_dists=components, 
                           shape=(Nf, Nr-1))

        delta = pm.Deterministic("delta",
                    pm.math.concatenate(
                        [delta0, -delta0.sum(axis=1, keepdims=True)], 
                        axis=1))
        delta_ac = pm.Deterministic("delta_ac",
                    pt.tensor.cumsum(delta, axis=0))
 
        r_ac = pm.Deterministic("R_ac",
                pt.tensor.mul(factors, delta_ac).sum(axis=1))
        rerror = pm.Normal("rerror", r_ac, 1E-6, observed=error)
        trace = pm.sample(cores=10)

Here ws_ is a matrix of (Nr-1, 3) with the proportions for each variable. alphas_x_ and betas_x_ are the corresponding parameters for the weibull distributions. Each has length (Nr-1).
I construct the last column as the sum of the previous (Nr-1) so the total change is 0.
factors is a (Nf, Nr) matrix with values that are given.
error is the vector of length Nf with the observed values.

The problem I’m facing is that

trace['posterior']['w'].mean(dim=['chain','draw','w_dim_0']))

gives a vector with 1 in the component corresponding to the higher value on ws_ and 0 otherwise, so the sample never touches other possibilities that are not the most likely one. Given the data and some other tests I’ve done, this should not be the case, there should be other possible combinations that explain the error term.

I’ll be very grateful on any suggestion on what I may be doing wrong.