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.

  • $\begingroup$ Dirichlet and multinomial random generation? Smells like LDA... $\endgroup$
    – mbq
    Commented Dec 13, 2011 at 8:08
  • $\begingroup$ @mbq: LDA? What is that? $\endgroup$ Commented Dec 13, 2011 at 8:16
  • $\begingroup$ Latent Dirichlet allocation. $\endgroup$
    – mbq
    Commented Dec 13, 2011 at 10:10
  • $\begingroup$ @mbq: Yes, something like that. In what context have you seen/used LDA? $\endgroup$ Commented Dec 13, 2011 at 17:10
  • 2
    $\begingroup$ @mbq: Yes, that (summing gammas) is what I decided to do. Thanks. $\endgroup$ Commented Dec 14, 2011 at 23:16

2 Answers 2


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.


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.


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.


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.

  • 1
    $\begingroup$ The multinomial distribution is called discrete, since $x$ only takes discrete values. You do not answer how to sample from the distributions. $\endgroup$
    – Olivier
    Commented Aug 1, 2016 at 0:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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