Non-centred parameterisation for Dirichlet distribution


I am trying to fit a model of the form:

for(i in 1:N)
   X[i] ~ multinomial(theta);
theta ~ dirichlet(exp(alpha));
alpha ~ normal(0,1);

where alpha is a univariate parameter. Unsurprisingly, this model – being non-centred – suffers from divergent iterations. These do not go away if I increase adapt_delta nor if I decrease the step size. I think I need a non-centred parameterisation here.

I’ve had a look at can’t seem to find one that people have used in the past for Dirichlet distributions. Does anyone have an idea how to ‘non-centre’ the above?



[edit: escaped code]


I think you want to integrate out theta in this situation, i.e. use a multinomial-Dirichlet likelihood for the data.


Thanks Ben, I hadn’t thought of that. It looks like Bob’s done the hard work for me of coding up the log probability:!topic/stan-users/uhrIEddZ2I8


In case it’s helpful, I’ve just had to implement Bob’s solution, so this is what your functions block would probably look like if you implemented a probability mass function for the Dirichlet-Multinomial:

functions {
  real dirichlet_multinomial_lpmf(int[] y, vector alpha) {
    real alpha_plus = sum(alpha);
    return lgamma(alpha_plus) + sum(lgamma(alpha + to_vector(y))) 
                - lgamma(alpha_plus+sum(y)) - sum(lgamma(alpha));


Thanks for posting. This is about as efficiently as it can be implemented without custom derivatives in C++.

You don’t need parens around return values. And it’s conventional to start new lines with operators rather than end them that way (easier to scan) in languages that allow it (R doesn’t as it isn’t whitespace neutral like Stan).


Thanks Bob, I’ve implemented your style suggestions.