由于积分表示意味着\Pr(\mathbf{x}\mid\boldsymbol{\alpha})是\Pr( \mathbf{x},\mathbf{p}\mid\boldsymbol{\alpha})。因此，生成\mathbf{x}的结果是生成\mathbf p，然后以\mathbf p为条件生成\mathbf{x}。这意味着生成机制是

$$\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}$$ $\Pr(\mathbf{x}\mid\boldsymbol{\alpha})$$\Pr(\mathbf{x},\mathbf{p}\mid\boldsymbol{\alpha})$$\mathbf{x}$$\mathbf p$$\mathbf{x}$$\mathbf p$

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

希望这两个步骤都很简单。

例如，这里是对应于上述的 R 命令：

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))

或者，当使用 rmultinom 时：

x=rmultinom(1,n,p)

举个例子

> 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