I am implementing a linear regression model in pymc3 where the unknown vector of weights is constrained to be a probability mass function, hence modelled as a Dirichlet distribution, as in the following code:

with pm.Model() as model:
    #prior on precision of normal likelihood
    tau = pm.Gamma('tau', alpha=1, beta=1)

    phi = np.empty(ncountries, dtype=object)
    y = np.empty((nyears-1, ncountries), dtype=object)
    for icountry, country in enumerate(countries):
        #prior Dirichlet allocation for each country
        phi[icountry] = pm.Dirichlet('mix_{c}'.format(c=country),
                                     np.roll(mix, icountry),

        for iyear, year in enumerate(years[1:]):
            suffix = '_{y}-{c}'.format(y=year, c=country)
            previous_pop = Xs[iyear, :]
            y[iyear, icountry] = pm.Normal('obs' + suffix,
                            'mu' + suffix,
                            dot(phi[icountry], previous_pop)),
                        observed=Ys[iyear, icountry])

After sampling the posterior by running:

    start = pm.find_MAP()
    step = pm.Metropolis()
    nsteps = 1000
    trace = pm.sample(nsteps, step, start=start)

I analysed the trace of the Dirichlet variables, and found that their values do not add to one (below is an example):

array([[ 0.01029745,  0.00627394,  0.00996922, ...,  1.83955829,
     0.00962185,  0.01020659],
   [ 0.01029745,  0.00627394,  0.00996922, ...,  1.83955829,
     0.00962185,  0.01020659],
   [ 0.01029745,  0.00627394,  0.00996922, ...,  1.83955829,
     0.00962185,  0.01020659],
   [ 0.02050308,  0.01685555,  0.01976797, ...,  1.92278065,
     0.03956622,  0.00473735],
   [ 0.01993214,  0.01632033,  0.01994876, ...,  1.92487858,
     0.04078728,  0.00481424],
   [ 0.01900882,  0.01528191,  0.02100671, ...,  1.92485693,
     0.0395159 ,  0.00524575]])

I am not familiar with theano variables, and found it difficult to explore how a Dirichlet RV is expressed in pymc3... Am I doing anything wrong, or should I just normalise the values returned in the trace so that they sum to one?

Quick update It looks like the function pm.find_MAP() employes a sort of gradient descent optimisation. This does not take into account the constraint resulting from the fact that a vector representing a draw from a Dirichlet distribution is a probability mass function (its values should be positive and their sum should be one). This constraint is apparently not enforced at the sampling stage of the algorithm either, and causes convergence problems as the precision of the likelihood distribution drifts towards zero.


We currently do not automatically transform the Dirichlet to be on the simplex. This is done explicitly using the simplex transform (from https://github.com/pymc-devs/pymc/blob/3eb2237a8005286fee32776c304409ed9943cfb3/pymc/examples/dirichlet.py#L10):

 p, p_m1 = pm.model.TransformedVar(
     'p', pm.Dirichlet.dist(a, shape=k),

Note though, that we have discussed to make this automatic: https://github.com/pymc-devs/pymc/issues/315


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.