如何理解狄利克雷过程(Dirichlet Process)

标签:#DirichletProcess##狄利克雷过程##非参贝叶斯# 时间:2018/01/04 20:10:37 作者:小木

[TOC]

假设我们有N个数据X=\{x_1, \cdots, x_N\},是来自于K个高斯分布的样本,那么所有的样本的联合概率分布是:

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

在这里,模型的参数是

\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表示,也就是

K \propto \log N

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

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

\theta_i \sim H(\theta)

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

G \sim DP(\alpha,H)

这里的H就是刚才那个H,称为基分布,另一个参数\alpha是一个正值,它的作用是决定G的离散程度,即GH多么的接近。实际上,当\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等。每个值上的竖线代表这个值的权重,从样本的抽取结果来看就是该样本出现的次数。


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

以下图为例(来自英文版维基百科,所以说把维基百科封禁真是个傻逼决策):

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

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

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

(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)表示Ga_1区域的权重总和,H(a_1)表示Ha_1区域的权重总和。也就是说每个DP的样本下的某种划分都是服从一个Dirichlet分布的。

下面,我们回忆一下Dirichlet分布的一个性质,假设:

P(X_1,\cdots,X_K) \sim DIR(\alpha_1,\cdots,\alpha_K)

那么,有:

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

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中,有:

E[G(a_k)] = H(a_k)

推导很简单,下面求和之后等于\alpha,和上面的一起划掉之后就剩下了这个了。我们发现\alpha并没有在均值中出现,也就是说这个\alpha对均值没有任何影响。接下来我们看一下方差:

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

\theta_1 \sim H

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

\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)如下:

\theta_2 \sim H

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

抽样

\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}

\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}

然而,

\frac{\sum_{m=1}^{K} n_{m,-i}}{n+\alpha-1} = \frac{n-1}{n+\alpha-1}

因此,

P(z_i=m|z_{-i}) = \frac{n_{m,-i}}{ n+ \alpha-1},existing

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构造法,将得到一个样本如下:

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中抽取一个样本,不仅是抽取了一系列参数(也就是分布),还抽取了每个参数(分布)的权重。

支付宝扫码打赏

如果文章对您有帮助,欢迎打赏鼓励作者