如何理解狄利克雷过程（Dirichlet Process）

[TOC]

假设我们有$N$个数据$X=\\{x\_1, \cdots, x\_N\\}$，是来自于$K$个高斯分布的样本，那么所有的样本的联合概率分布是：

```math
\log P(X) = \sum_{i=1}^N \log \sum_{l=1}^K \alpha_l N(x_i|\mu_l,\Lambda_l)
```

在这里，模型的参数是

```math
\Theta= \{ \mu_1, \cdots, \mu_K, \Lambda_1,\cdots, \Lambda_K, \alpha_1, \cdots, \alpha_{K-1} \}
```

对于$K$的的值是需要我们自己手动去指定的。但是，这在实际情况中是很困难的，尤其是数据在高维的情况下。因此，我们希望让模型自动确定这个$K$的数量。这个时候，一种简单的方式是可以假设$K$未知，把$K$当成一个位置参数，使用类似极大似然估计的方法来求解。这也是不行的，因为通过计算我们可以发现，当$K=N$的时候，似然最大，显然这不是我们需要的结果。那么，该如何确定这个$K$呢？这时候，我们其实希望得到一个关于$K$的表达式，它是关于$N$的一个函数，即$K=f(N)$。这里的$f$的形式很多，而我们期望的是只要随着$N$的增加，$K$的增长速度很慢，这样就可以找出一个合适的K了。那么其实这个函数可以用log表示，也就是

```math
K \propto \log N
```

就行了（其实在Dirichlet Process中，$E(K)\propto \log N$）。这个就是我们期望的$K$的表达式。

下面我们从另一个角度来理解这个问题，假设我们有一组数据$X=\\{x\_1, \cdots, x\_N\\}$，每个数据都是来自于参数为$\theta$的一个分布，也就是说这些数据都对应着一组参数$\\{\theta\_1,\cdots, \theta\_N\\}$。那么，其实这些参数也是需要有一个分布的，比如说假设：

```math
\theta_i \sim H(\theta)
```

但是注意到，如果这个$H$是一个连续的分布，那么从$H$中抽取的样本是不可能有完全一样的，也就是说上面这$N$个$\theta$都不可能相等。这也不符合我们聚类的要求了。那这个时候怎么办呢？我们想能不能构造一个离散的分布$G$，我们的$\theta$是从$G$里面产生的，但是呢，我们希望这个$G$与$H$很相似，这时候可以近似得到上述结果。怎么构造这个离散的分布呢？这就是Dirichlet Process的作用了。我们可以采用如下的方法来构造G：

```math
G \sim DP(\alpha,H)
```

这里的$H$就是刚才那个$H$，称为基分布，另一个参数$\alpha$是一个正值，它的作用是决定$G$的离散程度，即$G$与$H$多么的接近。实际上，当$\alpha=0$的时候，$G$是最离散的情况，它就是一个值。当$\alpha$值很大的时候，它就不是很离散。如果$\alpha$是无穷大的话，那么$G=H$。这里有一个视频演示：

https://upload.wikimedia.org/wikipedia/commons/9/93/Chinese_Restaurant_Process_for_DP%280.5%2CH%29.webm#t=0,00:00:25.866

它展示的是当$\alpha=0.5$的时候，$DP(\alpha,H)$产生的样本的类别的个数，也就是G的离散情况了。

大家可以根据上述操作步骤分别尝试使用$\alpha=0$和$\alpha=100$的时候样本的离散情况，会发现第一种只有一个桌子一直存在，而第二种几乎每个样本都是一个新桌子。地址分别如下：

$\alpha=0$的情况： http://topicmodels.west.uni-koblenz.de/ckling/tmt/crp.html?parameters=0.8&dp=1
$\alpha=100$的情况： http://topicmodels.west.uni-koblenz.de/ckling/tmt/crp.html?parameters=100&dp=1

那么这里的$DP(\alpha,H)$就是我们常说的狄利克雷过程了。我们这里说明一下DP的样本是什么。一般情况下，我们从一个分布抽取一个样本的结果是一个值，例如，我们从高斯分布中抽取一个样本，那么它就是一个值（向量）。我们可以用如下的图表示样本的抽取结果（自己鼠标画的，略丑别介意）。水平线是样本的值，这里比如说是0、0.1、-1.0等。每个值上的竖线代表这个值的权重，从样本的抽取结果来看就是该样本出现的次数。

<center>
![](http://www.datalearner.com/resources/blog_images/2e8ec19c-a433-4930-9a87-2d36f35ee8c2.png)</center>

那么对与DP来说，其样本不是一个值了，它是一个分布。每个DP的样本都可以当做是包含了无限个值及其对应的权重。那么这无限个值及其权重其实就是可以理解成一个分布。因此，我们说DP的样本是一个分布，即狄利克雷过程是一个关于分布的分布。一般的，DP的样本称之为一个测度（measure），大家可以理解成一个划分。什么意思呢？就是说假设我们有一个矩形，DP的样本就是矩形的一种划分方法。显然，一个矩形的划分有无数多种情况。

以下图为例（来自英文版维基百科，所以说把维基百科封禁真是个傻逼决策）：

<center>![](https://upload.wikimedia.org/wikipedia/commons/thumb/d/d3/Dirichlet_process_draws.svg/800px-Dirichlet_process_draws.svg.png)</center>

上面是9个图来自$DP(\alpha,N(0,1))$的样本，总共有4行，每一行从上倒下对应的$\alpha$分别是1、10、100和1000的实验。每一行是三个样本（参数一样）。我们可以看到，对于$\alpha=1$来说，它的样本是比较离散的情况，因为就集中在那几个值上，而随着$\alpha$的增加，分布越来越不离散了。

接下来我们继续说说DP产生的$G$的性质。$G$是一个随机变量。如下图所示，假设红色的那条线是基分布$H$，我们有两个样本分别是$G\_1$和$G\_2$。

<center>![](http://www.datalearner.com/resources/blog_images/7ff26de9-bab7-4908-9af3-6bb4f07ced33.png)</center>

关于这个G的特性可以通过如下方式来看，首先我们用几条绿色的竖线将每个G都划分成不同的区域，可以划分成任意多个的区域，每个区域大小都随便。每个区域我们都可以起个名字，比如我们划分成$K$个区域，$a\_1,a\_2,\cdots,a\_i,\cdots,a\_K$（紫色字母）。由于G是一个随机测度，假设那么每一个区域的权重总和（也就是竖线的总长度）分别表示成$G(a\_1)$、$G(a\_2)$也是随机的测度，它们也是有一个概率的性质的。这个概率的性质就是他们在一起服从一个Dirichlet分布：

```math
(G(a_1),G(a_2),\cdots,G(a_d)) \sim DIR (\alpha H(a_1), \alpha H(a_1), \cdots, \alpha H(a_d) )
```

这就是Dirichlet Process的定义了。注意，这里的$G(a\_1)$表示$G$在$a\_1$区域的权重总和，$H(a\_1)$表示$H$在$a\_1$区域的权重总和。也就是说每个DP的样本下的某种划分都是服从一个Dirichlet分布的。

下面，我们回忆一下Dirichlet分布的一个性质，假设：
```math
P(X_1,\cdots,X_K) \sim DIR(\alpha_1,\cdots,\alpha_K)
```
那么，有：

```math
E[X_i] = \frac{\alpha_i}{\sum_{k=1}^K\alpha_k}
```

```math
Var[X_i] = \frac{\alpha_i (\sum_{k=1}^K \alpha_k - \alpha_i) }{ (\sum_{k=1}^K\alpha_k)^2 (\sum_{k=1}^K\alpha_k + 1)}
```

那么，上面的DP中，有：
```math
E[G(a_k)] = H(a_k)
```
推导很简单，下面求和之后等于$\alpha$，和上面的一起划掉之后就剩下了这个了。我们发现$\alpha$并没有在均值中出现，也就是说这个$\alpha$对均值没有任何影响。接下来我们看一下方差：

```math

Var[G(\alpha_k)] = \frac{ \alpha H(a_k) (\alpha - \alpha H(a_k)) } {\alpha^2(
\alpha+1)} = \frac{ H(\alpha_k)(1-H(\alpha_k)) }{\alpha+1}

```

接下来我们看两个极端情况：
1、当$\alpha=\infty$时候，方差为0。也就是说G在某个区域的测度就等于H在这个区域的测度。
2、当$\alpha=0$的时候，$Var[G(a_k)] = H(\alpha_k)(1-H(\alpha_k))$，这个时候G在某个区域的分布就是一个贝努利分布，也就是说要么是在这个区域内，要么不是在这个区域内，这时候G是H的最离散的情况。

#### Dirichlet Process的构造
首先解释一下什么叫构造（Construction）。对于统计中的分布，我们都有一个概率密度函数，我们从这个分布中采样都是一个点。但是在Dirichlet process中，每个样本都是一个随机测度，它包含无数个点和相应的权重。这种采样是非常困难的，而且从定义中，我们也无法获取采样的方法。因此，我们需要有一个构造方法可以生成这样的测度。
在DP中，有三种常见的构造方法。

##### Stick-Breaking构造方法
就是掰棍子模型。在前面我们说了，每一个DP的样本都是一个随机测度，它由无数个点及其相应的长度组成，在统计中，这个点称为原子（Atom）。那么，在DP的一次采样中，我们是通过抽取无数个Atom及其权重得到的。我们来一步一步看如何采样这个Atom及其权重。首先，如下图所示，我们有一个基分布是H，图形是红色的曲线。下面那条黑色的线就是这个曲线对应的原子的组成，首先我们从H中随机抽取一个值$\theta\_1$，落在图中那个点的位置，这个位置对应着的就是一个Atom，其权重就是位置上竖线的高度。那么，第一个样本的Atom我们知道了，就是$\theta\_1$：

```math
\theta_1 \sim H
```

接下来我们要确定那个竖线的高度，也就是这个Atom对应的权重。如何抽取这个权重呢，我们首先从一个Beta分布中抽取一个值，其参数是$(1,\alpha)$，注意这个参数就是这个，不能变：

```math

\beta_1 \sim Beta(1,\alpha)

```

如图中绿色横线表示的是一个长度为1的线段，左边起点是0右边终点是1。我们从Beta分布中抽取一个点就是在这个线段中落下到一个位置$\beta\_1$上。这个位置与0之间的长度对应着刚才那个竖线的长度，也就是这个Atom的权重了，即$\pi_1=\beta_1$，这里的$\pi\_1$就是第一个Atom的权重。那么，第一个Atom和它的权重就已经采完了。就是$(\theta\_1,\pi\_1)$了。

第二次采样的时候，我们也依然先随机抽取一个原子，是$\theta\_2$，它对应的权重抽取方法是，先从$Beta(1,\alpha)$中抽取一个结果，得到$\beta\_2$，然后它的权重是$\pi\_2 = (1-\pi\_1) \beta\_2 $。意思就是把刚才第一次抽取结束的$\beta\_1$的位置到1终点的位置当做新的线段，然后抽取一个新的位置，得到了新的权重。那么新的权重的位置一定是在$\beta\_1$到1之间。如图所示。也就是说第二个Atom抽的结果$(\theta_2,\pi_2)$如下：

```math

\theta_2 \sim H

```

```math
\pi_2 = (1-\pi_1)\beta_1
```

对于后面的样本（Atom）继续按照上面的方式抽取，最终得到的所有的Atom及其权重就是一个G的实现结果了（realization）。
我们继续看一下$\alpha$的影响：
1.当$\alpha=0$的时候，$E[\beta_1]=1$，也就是说所有的权重都在第一个Atom上了，其他的权重都是0了。那么，也就是最离散的情况，所有的结果都是在$\beta_1$上。
2.当$\alpha=\infty$的时候，$E[\beta_1]=0$，那么所有的Atom权重都接近于0，那么这就是说明每个原子权重都差不多且非常接近于0，也就是说这个$G=H$了。

因此，如果G是DP的一个样本，那么它就是由无数个原子及其权重组成，G可以写成$G=\sum\_{i=1}^\infty \pi\_i \delta (\theta\_i)$

#### 抽样

```math
\begin{aligned}

P(\theta_i|\theta_{-i}) &= \int _G p(\theta_i,G|\theta_{-i}) dG \\

&\\

& = \int_G p(\theta_i|G, \theta_{-i}) p(G|\theta_{-i}) dG \\

&\\

& = \int_G p(\theta_i|G) p(G|\theta_{-i}) dG

\end{aligned}
```

```math
\begin{aligned}

P(z_i=m|z_{-i}) &= \frac{p(z_i=m,z_{-i})}{p(z_{-i})} \\

&\\

& = \frac{ \int_{p_1,\cdots,p_K} p(z_i=m,z_{-i}| p_1,p_2,\cdots, p_K) p(p_1,\cdots,p_K) }{ \int_{p_1,\cdots,p_K} p(z_{-i}|p_1,\cdots,p_K)p(p_1,\cdots,p_K) } \\

&\\

& = \frac{ \int_{p_1,\cdots,p_K} p(z_i=m,z_{-i}| p_1,p_2,\cdots, p_K) DIR(\alpha/K,\cdots,\alpha/K) }{ \int_{p_1,\cdots,p_K} p(z_{-i}|p_1,\cdots,p_K)DIR(\alpha/K,\cdots,\alpha/K) } \\

&\\

& = \frac{\Gamma(n_{m,-i} + \alpha/K +1) \prod_{l=1,l\neq m}^K \Gamma(n_{l,-i})}{ \Gamma( \alpha + n) } \times \frac{ \Gamma(\alpha+n-1)}{ \prod_{l=1,}^K \Gamma(n_{l,-i}) } \\

&\\

&=\frac{n_{m,-i}+\alpha/K}{ n+ \alpha-1} \\

&\\

&=\frac{n_{m,-i}}{ n+ \alpha-1}

\end{aligned}
```

然而，

```math
\frac{\sum_{m=1}^{K} n_{m,-i}}{n+\alpha-1} = \frac{n-1}{n+\alpha-1}
```
因此，
```math
P(z_i=m|z_{-i}) = \frac{n_{m,-i}}{ n+ \alpha-1},existing
```
```math
P(z_i=m|z_{-i}) = \frac{\alpha}{ n+ \alpha-1},new
```

#### 应用

举个例子，假设我们有一组数据$\textbf{x}$，每个数据$x\_i$都是来自于一个参数为$\theta\_i$的高斯分布，这个高斯分布的参数为$\theta\_i$，在贝叶斯统计中，所有的参数都是来自于一个分布，因此我们也假设这个高斯分布的参数$\theta\_i$来自于一个先验分布$H$。比如假设这个$H$是一个标准正态分布$N(0,1)$，从生成过程看，如果有$N$个数据那么我们从这个标准正态分布中抽样$N$次，每次的结果作为高斯分布的参数，继续抽取得到样本。显然这$N$个样本都是不一样的，尽管它们都是在$0$附近。如果，我们希望某些样本的高斯分布是相同的，也就是说依然需要从标准正态分布中抽取$N$个在$0$附近的样本作为高斯分布的参数，但是要求某些样本是相同的。那么，首先我们构造一个以标准正态分布为基分布的狄利克雷过程。有一个参数是集中参数为$\alpha$。首先，第一步，我们从这个DP中抽取一个样本，使用stick-breaking构造法，将得到一个样本如下：
```math
G = \sum_{k=1}^\infty \pi_k \delta_{\phi_k}
```

这里的 $\delta\_{\phi_k}$是一个以$\phi\_k$为中心的狄克拉函数，也就是指示函数，$\phi\_k$是来自于$H$的第$k$个样本，$\pi\_k$是其权重。那么，其实我们就得到了一个关于H的样本的离散情况，尽管每次抽取的$\phi\_k$都不一样，但是，根据DP的性质，在前面抽取的$\phi\_k$的权重应该是比较大的。可以这么理解，我们从标准正态分布中抽取前1000个样本，肯定也是0附近的数据出现的可能性高，也就是这些值越可能在前面抽样得到，也就是其权重可能更高。这样我们就能得到一个离散化的标准正态分布了，那么每个数据所属的分布的参数再从这个离散化的分布中抽取，其结果和从标准正态分布中抽取类似，但是不同数据之间是可能抽取到相同的参数的。也就是说，从DP中抽取一个样本，不仅是抽取了一系列参数（也就是分布），还抽取了每个参数（分布）的权重。

相关博客

最热博客