Sampling from Dirichlet-Multinomial

Question

How do I sample from a Dirichlet-Multinomial distribution. One which has this pmf:

http://en.wikipedia.org/wiki/Dirichlet-multinomial_distribution#For_a_multinomial_distribution_over_category_counts

Answer 1

Since $$\Pr(\mathbf{x}\mid\boldsymbol{\alpha})=\int_{\mathbf{p}}\Pr(\mathbf{x}\mid \mathbf{p})\Pr(\mathbf{p}\mid\boldsymbol{\alpha})\textrm{d}\mathbf{p}$$the integral representation implies that $\Pr(\mathbf{x}\mid\boldsymbol{\alpha})$ is the marginal of $\Pr(\mathbf{x},\mathbf{p}\mid\boldsymbol{\alpha})$. Generating $\mathbf{x}$ thus results from generating $\mathbf p$ and then generating $\mathbf{x}$ conditional on $\mathbf p$. This means a generating mechanism is

1. Generate $\mathbf p\sim\dfrac{1}{\mathrm{B}(\alpha)} \prod_{i=1}^K x_i^{\alpha_i - 1}$, the $\text{Dir}(\alpha_1,\ldots,\alpha_K)$ distribution;
2. Generate $\mathbf{x}|{\mathbf p} \sim \dfrac{n!}{x_1!\cdots x_K!}\,p_1^{x_1}\cdots p_K^{x_K}$, the $\text{Mult}(p_1,\ldots,p_K)$ distribution;

where both steps are hopefully straightforward.

For instance, here are the R commands that correspond to the above:

k=5  # dimension of the problem
n=75 # sample size
alpha=runif(k) # value of alpha, here chosen at random
p=rgamma(k,alpha) # pre-simulation of the Dirichlet
y=sample(1:k,n,prob=p/sum(p),rep=TRUE) # Multinomial
x=sum(y==1)
for (i in 2:k) x=c(x,sum(y==i))

or, when using rmultinom:

x=rmultinom(1,n,p)

with an example

> alpha
[1] 0.2216704 0.6642411 0.2528082 0.6309828 0.6942128
> p/sum(p)
[1] 0.02193043 0.07335277 0.23146885 0.00250276 0.67074519
> x
[1]  1  6 15  0 53