Looking for C/C++ implementations of sampling from multinomial and Dirichlet distributions

Question

I'm looking for C/C++ implementations of functions that return random variates multinomial and Dirichlet distributions. This is in the context of a calculation for posterior predictive p-values, part of which includes a MCMC step. I've been using Python, but am rewriting parts of my code in C++ for speed reasons. I've been using the implementations from numpy.random thus far.

My current options include

1. Pull the implementations from numpy.random, which are most likely written in C.

2. Use the R versions of these functions, which are conveniently available for Debian in the r-mathlib library. This library has been my default choice for such things for a long time, because the R people know their stuff when it comes to probabilistic computing.

However, I'd be open to other versions of these functions. Suggestions?

EDIT: r-mathlib doesn't appear to have a function for sampling from the Dirichet distribution, though according to Wikipedia's entry on the Dirichlet distribution I can sum Gamma variates. I wonder if that is a good way of doing this.

EDIT2: If possible, please comment on why you think the implementation you suggest is a good choice.

Answer 1

I would recommend the GNU Scientific Library (GSL).

GSL has both gsl_ran_dirichlet and gsl_ran_multinomial for the Dirichlet and multinomial distributions respectively. See the manual for the Dirichlet distribution is here and the manual page for the multinomial distribution is here. Full documentation can be found here or online here.

Answer 2

The reason to choose the following is because it does not depend on additional libraries. So, this makes it easier to support multiple targets such as arm or blackfin processors and reduces dependencies and ABI trouble.

Multinomial

In <random> there is a class:

template <class IntType = int> class discrete_distribution;

The term discrete distribution is actually incorrect, but we cannot expect STL programmers to be experts in statistics. :-)

So, this is a multinomial distribution without the preceding $1/B(x)$ factor.

$$f(x) = \prod_i p_i^{x_i}$$

The class accepts values $w_i$ which will be normalized such they form an actual probability $\sum_i p_i = 1$. Quite convenient in many cases.

Dirichlet

I haven't found a header file for the Dirichlet distribution. However, the Dirichlet is of the form:

$$f(x) = \frac{1}{B(\alpha)} \prod_i x_i^{\alpha_i-1}$$

So, you might want to check if you also can use the above class for this. As indicated in one of the comments. If the Dirichlet is used as a prior for a multinomial distribution, there are closed forms available for the posterior.