Calculating pdf of Dirichlet distribution in python

Question

I'd like to calculate the pdf for the Dirichlet distribution in python, but haven't been able to find code to do so in any kind of standard library. scipy.stats includes a long list of distributions but doesn't seem to include the Dirichlet, and numpy.random.mtrand lets one sample from it, but doesn't give the pdf.

Since the Dirichlet is pretty common, I'm wondering whether there's some other name I should be searching for it by in scipy.stats or similar, or whether I've just missed it somehow.

Answer 1

I couldn't find one in numpy, but it looked enough to implement. Here's an ugly little one-liner. (I followed the function given on Wikipedia, except you have to provide x = [x1, ..., xk] and alpha = [a1, ..., ak]).

import math
import operator

def dirichlet_pdf(x, alpha):
  return (math.gamma(sum(alpha)) / 
          reduce(operator.mul, [math.gamma(a) for a in alpha]) *
          reduce(operator.mul, [x[i]**(alpha[i]-1.0) for i in range(len(alpha))]))

Warning: I haven't tested this. Let me know if it works.

Answer 2

As of scipy version 0.15, you can use scipy.stats.dirichlet.pdf (see here)