狄利克雷过程混合模型（Dirichlet Process Mixture Model, DPMM）

狄利克雷过程混合模型（Dirichlet Process Mixture Model, DPMM）是一种非参数贝叶斯模型，它可以理解为一种聚类方法，但是不需要指定类别数量，它可以从数据中推断簇的数量。这篇博客将描述该模型及其求解过程。

一、狄利克雷过程混合模型模型介绍

我们有一组数据$$y_1,y_2,...,y_n$$，它们是相互独立的、来自于一些未知的分布中。$$y_i$$可以是多元的，它的值可以是实数的，也可以是类别型的。假设这些观测值来自于一个混合分布，其形式是$$F(\theta)$$，其参数$$\theta$$来自于狄利克雷过程（Dirichlet Process, DP）$$G$$，即该分布的参数的先验是一个狄利克雷过程，器集中参数为$$\alpha$$，基分布是$$G_0$$。关于狄利克雷过程的内容可以参见Dirichlet Distribution（狄利克雷分布）与Dirichlet Process（狄利克雷过程）。那么，狄利克雷过程混合模型（Dirichlet Process Mixture Model, DPMM）有如下：

y_i | \theta_i \sim F(\theta_i)

\theta_i | G \sim G

G \sim DP(G_0,\alpha)

这里简单解释一下为什么$$\theta$$的先验来自于狄利克雷过程。我们的数据是来自于一个混合分布，即$$n$$个数据是来自于$$k$$个分布，其中$$k \leq n$$。因此，对于来自同一个分布的观测值，其分布的参数$$\theta_i$$应该是一样的。但，如果我们的先验$$G$$是连续分布，那么，抽样结果不可能出现相同的值，而由于狄利克雷过程的特性，我们可以抽取相同的参数$$\theta$$且使得其与基分布相似，具体可参照我们另外一篇博客。同时，由于该模型数据增长速度远大于混合模型个数的增长，因此，对于一定的样本数据，我们可以得到远小于样本的类别个数，因此该模型可以作为非参数聚类方法。 对于$$\theta$$的抽取，我们可以看作是一个从多项式分布中抽取的过程，即

\theta_i | \theta_1,...,\theta_{i-1} \sim \frac{1}{i-1+\alpha} \sum_{j=1}^{i-1}\delta(\theta_j)+\frac{1}{i-1+\alpha}G_0

其中$$\delta(\theta_j)$$与$$G_0$$可以理解为$$i$$个概率密度，各个概率密度前面的系数加和为1。因此，$$\theta_i$$可理解为来自于$$i-1$$个已存在的类别与一个新的的类别$$G_0$$。

因此，前面的模型形式也可以写成如下形式（第二种模型）：

y_i | z_i,\phi \sim F(\phi_{z_i})

z_i | \bold{p} \sim Multi(p_1,p_2,...,p_K)

\phi_z \sim G_0

\bold{p} \sim Dirichlet(\alpha/K,...,\alpha/K)

这里的$$z_i$$是因变量，是观测值$$y_i$$所属的类别。对于每个类，$$z$$都有一个对应的参数$$\phi_z$$，决定了这个类的分布情况。各个类别的混合比例是$$\bold{p}=(p_1,...,p_k)$$，来自于参数是$$\alpha/K$$的狄利克雷先验分布。这里的类别数量$$K$$可以趋向于无穷大。

我们看一下类别的概率。假设已有的观测值分别属于$$i-1$$个类别，那么接下来改观测值属于某个类别的概率为：

\begin{aligned}
P(z_i=z|z_1,...,z_{i-1}) &= P(z_1,...,z_i) / P(z_1,...,z_{i-1}) \\
&= \frac{ \int \prod_{k=1}^{K-1} p_{z_i} \cdot  p_c \cdot \bold{B}({\alpha / K }) \cdot \prod_{k=1}^{K}p_{k}^{\alpha/K-1} dp}{\int \prod_{k=1}^{K-1} p_{z_i} \cdot \bold{B}({\alpha / K }) \cdot \prod_{k=1}^{K}p_{k}^{\alpha/K-1} dp} \\

& = \frac{n_{i,z}+\alpha/K}{i-1+\alpha}
\end{aligned}

其中，$$\bold{B}({\alpha / K }) = \frac{\Gamma(\alpha)}{\Gamma(\alpha/K)^K}$$，这里的$$n_{i,z}$$是类别为$$z$$的观测值的数量，$$i$$为所有观测值数量。当我们的$$K \to \infty$$时，有：

P(z_i = z | z_1,...,z_{i-1}) \to \frac{n_{i,z}}{i-1+\alpha}

P(z_i \neq z_j | z_1,...,z_{i-1}) \to \frac{\alpha}{i-1+\alpha}

因此，上述概率对于定义狄利克雷过程混合模型来说是必要的。因为只有这样的概率定义才能使得该模型为非参数的贝叶斯模型。

二、使用吉布斯抽样求解狄利克雷过程混合模型

对狄利克雷过程混合模型的后验期望直接进行计算非常困难，这里我们可以使用吉布斯抽样的方法（参见机器学习中MCMC方法介绍）。假设我们有来自后验分布$$\theta=(\theta_1^{(t)},...,\theta_{n}^{(t)})$$的样本T。最直接的方法是我们重复的在给定数据和其他$$\theta_j$$（这里$$j \neq i$$，记作$$\theta_{-i}$$）的情况下，抽取$$\theta_i$$。这个条件概率可以通过似然函数$$F(y_i,\theta_i)$$和$$\theta_{-i}$$先验的乘积获得，$$\theta_i$$的条件先验如下：

\theta_{i} | \theta_{-i} \sim \frac{1}{n-1+\alpha} \sum_{j \neq i} \delta(\theta_j) + \frac{\alpha}{n-1+\alpha}G_{0}

那么，给定数据和其他参数的情况下的概率为：

\begin{aligned}
\theta_i | \theta_{-i},y_i &= \frac{ (\theta_i | \theta_{-i}) \cdot (\theta_{-i}) \cdot F(y_i,\theta_i) }{(\theta_{-i},y_i)} \\
&\sim C \cdot (\theta_i|\theta_{-i}) \cdot F(y_i,\theta_i) \\
&=C_1 \cdot \sum_{j \neq i} \delta(\theta_j) \cdot F(y_i,\theta_i) + C_2 \cdot G_0 \cdot F(y_i,\theta_i)
\end{aligned}

因此，我们有：

\theta_i | \theta_{-i} \sim \sum_{j \neq i} q_{i,j} \delta(\theta_j) + r_i H_i

该公式称为原始类别确定公式。其中，$$H_i$$是先验为$$G_0$$，单观测值为$$y_i$$的$$\theta$$的后验，其似然函数为$$F(y_i,\theta)$$。$$q_{i,j}$$与$$r_i$$的值分别如下：

q_{i,j} = b \cdot F(y_i,\theta_j)

r_i =b \cdot \alpha \int F(y_i,\theta_i) dG_0(\theta)

其中，$$b$$是常数，使得结果满足$$\sum_{j\neq i} q_{i,j} + r_i =1$$。如果要使得吉布斯抽样算法可以工作，必须使得$$r_i$$的积分可以求解，且从$$H_i$$中抽样是可行的。这里我们得到了求解算法1。

求解算法1

对于$$i=1,...,n$$：从上述原始类别确定公式$$\theta_i|\theta_{-i},y_i$$抽取一个新值。

但是该过程收敛速度太慢，使得该算法效率很低。这是因为相同$$\theta$$情况下，各个观测值是相互联系的。但是，该算法无法同时改变超过一个观测值的$$\theta$$。 该问题可以通过更改抽样对象来解决，我们可以从上述第二种模型来进行吉布斯抽样。我们可以考虑将多项式概率$$\bold{p}$$积分掉，然后进行抽样。当$$K \to \infty$$时，每次吉布斯抽样可以变成在给定$$y_i$$、$$\phi_z$$和$$z_j$$时，为$$z_i$$选择一个新的值。然后，当$$z_i=z$$，在给定$$y_i$$的情况下，从条件分布中为每个$$\phi_z$$选择一个新的值。这里$$z_i$$的条件概率很容易计算：

P(z_i=z|z_{-i},y_i,\phi) = b \cdot F(y_i,\phi_z) \cdot \frac{ n_{-i,z} + \alpha/K }{ n-1+\alpha }

这里的$$b$$是正则化常数。 当$$K \to \infty$$时，我们当然无法写出所有的$$\phi_z$$。这里我们可以只对和观测值有关的类别坐吉布斯抽样。因此，当$$z$$是来自于现有的类别时，即$$z=z_j$$时：

P(z_i=z|z_{-i},y_i,\phi) = b \cdot \frac{n_{-i,z}}{n-1+\alpha} \cdot F(y_i,\phi_z)

当$$z$$来自一个新的类别时，即$$z_i \neq z_j$$：

P(z_i=z|z_{-i},y_i,\phi) = b \cdot \frac{\alpha}{n-1+\alpha} \cdot \int F(y_i,\phi_z) dG_0(\phi)

这两个式子我们这里称为类别确定公式。该过程可描述如下求解算法2。

求解算法2

假设马尔可夫链是$$\bold{z}=(z_i,...,z_{n})$$，$$\phi=(\phi_z : z \in \{z_1,...,z_n\})$$。重复抽样过程如下： 1、对于$$i=1,...,n$$：如果当前类别不包含任何观测值，那么去掉$$\phi_z$$。从上述类别确定公式$$z_i|z_{-i},y_i,\phi$$中抽取一个类别，如果新类别与已有类别都不一样，则从$$H_i$$中抽取一个$$\phi_{z_{i}}$$。 2、对于$$z\in \{z_1,...,z_n\}$$：对于$$z_{i}=z$$，从$$\phi_{z} | \text{all} \space y_i$$中抽取一个新的值。

通常情况下，我们都可以积分掉$$\phi_{z}$$，把它们从算法中消除，最后的马尔科夫链就变成了只与$$z_{i}$$有关，也就是我们使用如下的条件概率来进行吉布斯抽样： 当$$z$$是来自于现有的类别时，即$$z=z_j$$时：

P(z_i=z|z_{-i},y_i) = b \cdot \frac{n_{-i,z}}{n-1+\alpha} \cdot \int F(y_i,\phi_z) dH_{-i,c}(\phi)

当$$z$$来自一个新的类别时，即$$z_i \neq z_j$$：

P(z_i=z|z_{-i},y_i) = b \cdot \frac{\alpha}{n-1+\alpha} \cdot \int F(y_i,\phi_z) dG_0(\phi)

这里我们称为新类别确定公式，因此，我们得到了算法3。

求解算法3

马尔可夫链为$$\bold{z}=(z_i,...,z_{n})$$，重复如下过程的抽样： 对于$$i=1,...,n$$：从新类别确定公式$$z_i|z_{-i}$$中抽取一个新值。