Dirichlet process (三)

今天我们讨论一下Dirichlet process在Bayesian nonparametric里的具体应用,也即它的Posterior的形式。由于电子计算机还未大力发展,在有效sampling手段仍然缺乏的70年代,Bayesian的研究一大块是注重于构造conjugate的prior。所谓conjugate也即要求Posterior和prior能在同一个family里。因为一般的prior都是属于比较容易sample的分布,比如Gaussian啊,或者Gamma之类的,如果Posterior也有同样的形式,而改变的仅仅是参数,那么我们就可以快速求出posterior的均值方差,甚至是从里面sample。Ferguson在构造Dirichlet process的时候也是本着同样的考虑。他选择Dirichlet distribution作为每一个m partition上的prior,一方面因为这是一个自然的选择,另一方面也是因为Dirichlet distribution是multinomial distribution的conjuage prior。也即假如我们的model是一个multinomial distribution

x\sim multi(p_1, p_2, \cdots, p_m)

那么我们如果在p_1, p_2, \cdots, p_m上加一个Dirichlet priorDirichlet(a_1, a_2, \cdots, a_m),其Posterior的形式 (通过简单计算) 就是

p|x \sim Dirichlet(a_1 + \delta_1(x),\cdots, a_m + \delta_m(x))

其中\delta_i(x)取1,如果x = i,否则取0。意思就是x落在哪个category里,Posterior里那个category的参数上就加1。类似的,如果我们有n个数据,那么它的posterior就应该是

p|x \sim Dirichlet(a_1 + \sum_{j = 1}^n\delta_1(x_j),\cdots, a_m + \sum_{j = 1}^n \delta_m(x_j))

对于Dirichlet process,我们也有同样的结论。下面我们从两个不同的角度来论证Dirichlet process的posterior。我们考虑Dirichlet process的infinite sum representation,

DP(x) = \sum_{i = 1}^\infty p_i\delta_{x_i}(x)

根据上一篇的结论,从这里sample出来的分布函数都是离散分布,所以它所model的数据的分布是这样的\mathbb{P}(x = x_i) = p_i,其Likelihood为,

L(x) = p_{x_1}^{\delta_{x_1}(x)}\cdot p_{x_2}^{\delta_{x_2}(x)}\cdots = \prod_{i = 1}^\infty p_{x_i}^{\delta_{x_i}(x)}

而Dirichlet process的density有类似的形式。对于参数为\alpha(\cdot)的DP,其p_i, i = 1,2,\cdots由一个intensity为\alpha的Gamma process决定,而x_i由一个\alpha(\cdot)/\alpha的概率分布决定,所以其density (如果存在的话)应该有如下的形式,

\mathcal{P} = p_{x_1}^{\alpha(x_1)}\cdot p_{x_2}^{\alpha(x_2)}\cdots = \prod_{i=1}^\infty p_{x_i}^{\alpha(x_i)}

那么根据Bayes rule,我们知道posterior就应该是

DP|x \propto L(x)\cdot \mathcal{P} = \prod_{i=1}^\infty p_{x_i}^{\alpha(x_i) + \delta_{x_i}(x)} = \prod_{i=1}^\infty p_{x_i}^{\alpha(x_i) + \delta_{x}(x_i)}

也即,posterior是一个参数为\alpha(\cdot) + \delta_x(\cdot)的Dirichlet process。当然更严格的证明应该类似于Ferguson论文里,考虑任何一个m-partition,证明其posterior在这个m partition上服从一个Dirichlet distribution,且其参数为\alpha(A_i) + \delta_x(A_i)。更一般的,如果我们观测到n个数据点,那么我们的posterior就应该是参数为\alpha(\cdot) + \sum_{i = 1}^n \delta_{x_i}(\cdot)的 Dirichlet process。

我们还可以从另一个角度看这个问题。在上一篇里我们提到了Blackwell和McQueen证明了DP和Polya urn scheme之间的等价性。注意Polya urn scheme是这样定义的:你有一个基准测度\alpha(x),然后第一次从里面sample一个点x_1,然后更新测度为\alpha(x) + \delta_{x_1}(x)。下一次从这个新测度里sample x_2并更新测度为\alpha(x) + \delta_{x_1}(x) + \delta_{x_2}(x),依次类推。这个不断重复的过程定义出来的概率测度就是Polya urn scheme,同时也是DP。这种构造其实提供了一种很简单的计算DP posterior的方法。

DP(x)|x_1,\cdots, x_n = Polya(x)| \mbox{ The first n points sampled are } x_1, x_2, \cdots, x_n

由于Polya urn的自相似性,从第n+1步开始的过程仍然是Polya urn,只是基础测度变成了\alpha(x) + \sum_{i = 1}^n \delta_{x_i}(x)。这也就意味着posterior应该是参数为\alpha(\cdot) + \sum_{i = 1}^n\delta_{x_i}(\cdot)的Dirichlet process。

至此我们就已经得到了在观测到数据点后DP posterior的形式。我们来看一看它跟frequentist之间的联系。很明显,DP的posterior就是frequentist估计概率分布时用的empirical distribution。区别就是DP在empirical distribution前面多了一个\alpha(\cdot)测度,相当于研究者对于数据分布的initial guess。正因为这两者的相似性,对于一维的数据,DP posterior可以直接用来构造很多Nonparametric testing,详细的可以参考Ferguson的paper。同时,Frequentist还有另一种估计概率分布的方法,那就是kernel density estimation (KDE)。每个kernel是一个Gaussian或者其他的分布,这些Gaussian的mean就在某一个观测到的数据点上,还有一个bandwidth参数需要调节。KDE对比于empirical distribution的好处是,它是smooth的,而且是连续随机变量的分布。DP也可以类比定义出类似的Bayesian kernel density来,具体就留给读者思考了。

在这三篇blog里,我们简单介绍了最原始的Dirichlet process,讨论了motivation以及一些性质并推导了Posterior。以后可能还会简单的写一篇blog介绍一些DP的衍生物,比如Chinese restaurant process和Indian buffet process。

思考:
1. 利用第一个思路严格证明DP的posterior也是DP而且其参数为\alpha(\cdot) + \delta_{x}(\cdot),其中x是观测到的数据点。

2. 写出与kernel density estimation对应的Dirichlet process model连续分布函数的方法。

Dirichlet process (二)

这一篇当中我们就聊一聊Ferguson老爷子是怎么构造Dirichlet process来作为所有概率分布上的prior的。就像我们之前谈到,由于概率分布函数的空间太大了,我们没有办法直接定义概率。一个很直观的解决办法就是考虑里面一个较小的集合,而这个集合在原空间里又是在某种意义下稠密的。想象一下在《实变函数》里,当我们处理一般的可测函数比较困难的时候,我们就先考虑简单函数或者阶梯函数来近似,而这些近似的函数构成的集合往往都比较小且可数,技术上容易处理。类似的想法也可以应用到这里。假如我把空间partition成一个个小块,原来的分布函数都可以算出在这些小块上的概率。假如我用这些块状的分布函数去近似原来空间里的分布函数,那么近似的结果就会随着分块的细化而越来越好 (比如每个block上的概率值的差别)。这个观察意味着我们可能只需要考虑限制在这些小块上的分布函数就够了。而这些块状的函数的好处是他们的构成的空间很小。

假设我们分割空间使得

\mathcal{X} = A_1\cup A_2\cup \cdots \cup A_m.

A_i\cap A_j = \emptyset.

那原来空间里的任何一个概率分布函数都可以被map成一个m维的向量

\mathbb{P} \rightarrow (\mathbb{P}(A_1),\cdots, \mathbb{P}(A_m)).

这就相当于强行把所有概率函数给嵌入到了m维的欧式空间里。好了,提问!假如我们能够利用partition把所有可能的概率函数嵌入到m维的欧式空间里,那么原来概率函数空间上的prior在这个m维的欧式空间里最合理的形态应该是什么呢?

答案就是Dirichlet distribution.

因为嵌入的是在m个小块上的概率,所以这些向量的每个分量都必须是正的,同时因为他们是概率,所以和必须是1。也就是说,假如我们把原来空间里的概率函数都强行嵌入到一个m维的欧式空间里 (通过分割空间成m块),那么原空间上的Prior对应到这个嵌入空间上的Prior就应该是一个Dirichlet distribution

(\mathbb{P}(A_1),\cdots, \mathbb{P}(A_m)) \sim Dirichlet(\alpha(A_1), \cdots, \alpha(A_m)).

其中\alpha(\cdot)是用户指定的参数函数,来定义每一块上的Dirichlet的参数值。这一切看起来都非常自然,但是我们现在依然面临着两个问题。第一,我们希望让m趋向于无穷,这样我们就可以无限细分空间,使得我们的block function可以任意接近原来空间里的任何一个概率函数。第二就是我们需要在所有的这些1维,2维,…,m维,… 概率函数上定义测度。上面的例子是在给定m和partition的情况下定义出来的概率函数。但是即使在m固定的情况下,我们也有无穷多种partition的可能性,而每一个不同的partition就对应了一个不同的概率函数。我们把所有这些概率函数 m = 1, 2, 3,\cdots 记为集合M。(虽然M包含如此之多的概率函数,它依然远小于我们最早提到的所有概率函数构成的空间)。我们需要在M以及其\sigma-域上建立概率测度。因为这个概率测度如此特殊 (建立在概率函数上, distribution over distribution),它也被称为random measure (而在Ferguson的paper里被称为random probability).

其结果就是Dirichlet process。具体的构造过程其实非常简单,也就是我们前面的步骤,我们只需要在每一个有限Partition (注意一个Partition就是一个概率函数) 上利用\alpha(\cdot) 定义一个Dirichlet distribution就可以了。随着m趋向于无穷,我们就定义了一个point process。文章里简单verify了一下Kolmogorov consistency条件,也即这些Marginally定义的Dirichlet distribution之间不会有矛盾 (我把某个Partition里的一个小块和旁边的合并得到的distribution,应该和我之前直接在这种新partition下定义的distribution一致)。我们主要来看看Dirichlet process的一些直观的性质。



1. 首先我们先来看看有限步的情况,也即这个Process的第m步对应的随机变量是啥。根据定义,我们知道这个随机变量的support应该是所有的m分块的概率函数。这个随机变量应该有两部分组成,一部分是每一块上的概率,另一个是具体的分块。很容易就可以注意到,这两部分其实是互相独立的。所以如果我们把它写出来,其形式大概就是

D_m(x) = p_1 \delta_{A_1}(x) + p_2 \delta_{A_2}(x) + \cdots + p_m \delta_{A_m}(x),

其中\delta_{A_i}是一个指示函数,表示一个数值是否落在集合里,而A_1, \cdots, A_m是一个组随机partiiton。另一方面,p_1,\cdots, p_m是各个集合上的概率,根据我们的定义应该服从Dirichlet distribution。



2. 相同的idea可以被应用到m趋于无穷时的情况,Dirichlet process应该也有相似的表达。当空间被无限细分之后,\delta_{A_i}就自然从一个集合的指示函数退化到了一个单点的delta函数,所以我们有

DP(x) = p_1\delta_{x_1}(x) + p_2\delta_{x_2}(x) + \cdots + p_m\delta_{x_m}(x) + \cdots

其中x_i就是空间里随机抽取的无穷多个点 (有点类似于point process,但是以概率1有无穷多个)。困难之处是现在p_1, \cdots, p_m, \cdots的分布。这是一个无穷维的向量,所以现在不能用Dirichlet distribution来刻画了,为了得到Dirichlet process在这组向量上诱导的分布,我们需要考虑Dirichlet distribution的一个等价定义。假如我们有m个独立变量v_i = Gamma(a_i, 1), i = 1,2,\cdots, m,那么我们就有

\bigg(\frac{v_1}{\sum v_i}, \frac{v_2}{\sum v_i}, \cdots, \frac{v_m}{\sum v_i}\bigg) \sim Dirichlet(a_1, a_2, \cdots, a_m).

这个等价定义提示我们也许可以通过定义一个Gamma process来得到p_1, p_2, \cdots上的概率。举例来说,假如我们的这个Gamma process的每个元素都是一个独立的Gamma(a_i, 1),且我们有\sum_{i = 1}^\infty a_i = \alpha,那p_i就可以简单的定义为\frac{Gamma(a_i, 1)}{\sum_{i = 1}^\infty Gamma(a_i, 1)}。当然这只是一个直观解释,实际构造的Gamma process要更复杂一些。注意Gamma process并不要求每一个增量都是Gamma variable,只要他们的和是就可以了。在Ferguson原文里用了一个很复杂的counting process来构造,具体如下。我们构造第一个变量J_1使得

\mathbb{P}(J_1\leq x_1) = \exp\{N(x_1)\},

其中

N(x) = \exp\{-\alpha\int_x^\infty e^{-y}y^{-1}dy\}.

然后我们就可以sequentially构造J_i使其服从类似的分布,

\mathbb{P}(J_i\leq x_i|J_1 = x_1, \cdots, J_{i - 1} = x_{i - 1}) = \exp\{N(x_i) - N(x_{i - 1})\},

其中x_1 > x_2 >\cdots > x_{i - 1} > x_i。可以看出来这样构造的序列符合Dirichlet process的self consistency的性质,同时Ferguson还证明了\sum_{i = 1}^\infty J_i是一个带独立增量的Gamma process,所以最后利用这个构造我们就可以得到p_i的分布了

p_i = J_i/(\sum_{i = 1}^\infty J_i) \label{eq:1}

另外这个构造还有一个好处就是这些概率p_i一定是单调下降的!但是问题就是sampling比较复杂。我们后面可能会在讲一讲Dirichlet process的另一个构造性证明,也即stick-breaking process (Sethuraman, 1994)。在stick-breaking process的构造里这些p_i不需要单调下降,所以sampling简单很多。同时不同于Gamma process的构造,stick-breaking process是直接sample比值 J_i/(\sum_{i = 1}^\infty J_i),比起Gamma process的构造也是要简单一些。



3. 我们在2当中定义的其实就是Dirichlet process的infinite sum representation,

DP(x) = \sum_{i = 1}^\infty p_i\delta_{x_i}(x)

而且也有了p_i的分布了,似乎已经大功告成。不过等一下,好像我们忘了在一开始我们还有一个用户指定的参数测度\alpha(\cdot)用来指定m个block的Dirichlet分布的参数值的。然而现在p_i里面只有一个常数\alpha是用来把\alpha(\cdot) normalize成概率测度的,那剩下的\alpha(\cdot)/\alpha怎么加回来呢?
注意啦,当我们从DP里sample了一个概率分布后,我们计算它在各个block上的概率值的时候实际是看的有多少x_i落在那个block里 (p_i没法影响这些block概率,因为它的assignment是和location独立的,是exchangeable的),所以很显然我们应该把这个用户测度\alpha(\cdot)/\alpha加在x_i上!加上x_i是independent and identically distributed变量,且服从\alpha(\cdot)/\alpha分布,我们就完成了完整的Ferguson的DP infinite sum representation的定义。这个形式比起之前抽象的描述在实际当中更有用处。在以后的使用中,人们把\alpha称为concentration parameter,而把\alpha(\cdot)/\alpha称为base measure。



4. 从DP的infinite sum representation里我们可以看出,DP是定义在所有discrete distribution上的。因为每个realization都是一些delta函数的无穷和。所以DP很重要的性质就是,从DP里sample出来的概率分布函数以概率1是离散分布,虽然它可以无限近似任何连续分布。另一方面,在上面对DP的构造中,很自然可以看出DP和Polya urn之间的联系,因为是不断在抽取一些discrete atom。如果我们把x_i都看成是一种颜色的小球,那DP就是在一个urn里不断抽取球。一开始我们的urn里有\alpha(\mathcal{X})那么多球。然后我们抽取第一个小球x_1,我们观察到了这个颜色,然后我们把同样颜色的另一个球x_1加回去,这样球的总数就变成了u = \alpha(\mathcal{X}) + \delta_{x_1},然后我们继续进行从这个新的urn里抽取小球,并且每次放回另一个一样颜色的小球,最后我们的分布就变成了

u = \alpha(\mathcal{X}) + \delta_{x_1} + \delta_{x_2} + \cdots

注意这些x_i有一些可能是一样的。在同一年Blackwell和McQueen的一篇文章里,他们讨论了这种Polya urn scheme和DP之间的关系,并且证明了他们的等价性。这种Polya urn形式的DP和后来的Chinese resuturant process之间有着很大的联系 (每个颜色的小球就相当于Chinese restaurant里的一个个table),它们都有着这种rich getting richer的性质,也即早些被抽取的小球以后被抽取的概率也会比较大。



这次就先介绍到这里吧,内容似乎十分多的样子(扶额…)。下一次的讨论会集中在DP的posterior上,以及和传统的frquentist的nonparametric test之间的联系,相信会是一次比较简短愉快的讨论。

思考:
1. 证明infinite sum representation形式的DP可以以正的概率落在任何一个概率分布Q的任何可测partition的\epsilon包里,只要这个概率分布和\alpha(\cdot)是绝对连续的。即证明对任意概率函数以及任意partition A_1, A_2, \cdots, A_m

\mathbb{P}(\max_{A_i}|DP(A_i) - Q(A_i)|\leq \epsilon)>0.

2. 证明DP并不能以total variation distance无限接近任何概率函数,也即存在概率分布Q,使得对任意\epsilon > 0

\mathbb{P}(\max_{A\subseteq\mathcal{X}}|DP(A) - Q(A)|\leq \epsilon) = 0

3. 证明 \bigg(\frac{Gamma(a_1, 1)}{Gamma(\alpha, 1)}, \cdots, \frac{Gamma(a_m, 1)}{Gamma(\alpha, 1)}\bigg)\sim Dirichlet(a_1, a_2, \cdots, a_m),其中\sum_{i = 1}^m a_i = \alpha

4. 证明Ferguson定义的Gamma process得到的概率值p_i一定是单调下降的。

Dirichlet process (一)

这是最早定义Dirichlet process的文章,发表在1973年的Annals of Statistics上。在当时计算机还没有普及的年代,Bayes的方法基本局限在conjugate的框架之下,这也是Ferguson大神在文章开头就指出的定义Prior需要考虑的条件之一:posterior要足够好算(要有close form)。结果即使在如今计算机十分普及的情况下,人们反而更关注计算速度,于是conjugate的prior使用仍最为广泛。因此当时的Ferguson的考虑结果到如今依然有意义。

Ferguson在文章的开头就指出了,non-parametric Bayes最大的问题就在于没有合适的prior。本质上来说统计的所有估计问题,都是在估计随机变量X的分布函数。在参数模型里我们会假设一个具体的分布,比如

X\sim N(u,\sigma^2).

而正态分布可以被均值和方差唯一决定,于是我们所要做的就是估计\mu, \sigma这两个参数。但是使用参数模型的风险是我们选的模型可能不对,比容把正态模型用到非正态上。非参数模型里不会具体假设一个参数模型,这样就避免了选错模型的问题,但是却带来了更多的问题。在参数模型里,以正态模型为例,假设所有正态分布构成的集合是A,我们默认待估计的分布函数是在A中的,我们只要找出A中那一个最好就可以了(从数据来判断哪个最好)。而正态分布可以用两个参数唯一确定。这样我们就把A中的每个元素和二维欧式空间的点一一对应了起来

N(\mu,\sigma^2)\leftrightarrow (\mu,\sigma),

—–就这个意义来说集合A其实不是很大,应该说很小,比一个二维欧式空间还小。但是在非参数模型里,我们不再默认待估计的分布函数是局限在正态里,现在集合A的元素是任何分布函数—取值在0,1之间的单调右连续函数(更准确地说是说在给定概率空间上的所有测度函数,或者是所有随机变量)。现在集合A有多大呢?至少欧式空间是放不下A的。那么集合A太大带来的问题是什么呢?

为了说明这个问题我们引入一些泛函(functional)的观点。我们把likelihood看成是定义在集合A上的一个函数,因为给定集合A中的任何一个元素(分布函数),我们都可以写出对应的likelihood,并算出相应的值。所以likelihood是把A中元素映到R上的一个映射。由于在这种情况下自变量是函数,而非传统的数值,因此这种函数被称为泛函。对于frequentist来说,他们所用的估计方法是maximum likelihood,即找到A中的一个元素,它刚好最大化likelihood,这个元素就是MLE。但是现在集合A不能被嵌入欧式空间当中,这就意味着原来定义的导数,极限的概念都不能用了,从而给这个最优化问题带来了困难。对于最优化而言,集合A最大的问题是没有“距离”的概念(更多的时候还需要“内积”的概念),所以第一步是要在A上引入距离和内积。但是即使有了距离,欧式空间上原来简洁的导数法则也不再适用了,所以我们还是要尽量把A往欧式空间上靠。幸运的是,如果距离和内积选的比较好,那么我们的集合A就具有了“可分”(separable)的性质。“可分”的意思就是说,我们可以在集合A里面找到一个非常小(可数无穷多个元素)的集合A’,使得其在A中稠密(dense,对于A中的任何一点都可以找到A’中的任何一点,使其任意接近)。如此A’的基(basis)也就成了A的基,从而所有A中的元素都可以唯一表示成集合A’的基的线性和的形式。因为A’只有可数无穷多个,这组基也只可能是可数无穷多。这样来看的话,虽然A不能嵌入欧式空间中,但是A比欧式空间也大不了多少,相当于是无穷维的欧式空间罢了,而且还是线性空间。因此欧式空间中的结果都可以自然推广到A上了。

这样讲比较抽象,我们举一个具体的例子。现在假设我们现在想估计变量X的分布函数,简化版的Weierstrass Theorem告诉我们任何连续密度函数都可以用不同参数的Gaussian density的加权和去无限逼近(scale mixture of Gaussian)。这个定理证明很简单,可以参看http://www.math.sc.edu/~schep/weierstrass.pdf,把最上面等式的积分离散化我们就得到了,

f(x) = \lim_{h\rightarrow 0} \sum_{k=-\infty}^\infty \frac{1}{\sqrt{2\pi}h}e^{-\frac{(x-kh)^2}{2h^2}}\cdot hf(kh).

这里我们的A’就是一小部分Gaussian分布(N(kh,h), k\in\mathcal{Z},可数多个),其中hf(kh), k\in\mathcal{Z}就是需要估计的权重,而A就是所有有连续密度函数的分布函数。我们可以看出来它其中的元素都可以被表示成A’中元素的线性和。于是通过选择A’作为基,A中的任何一个元素f都对应于一组坐标(\cdots, hf(-2h), hf(-h), hf(0), hf(h),\cdots)

对于Bayes来说,情况就略有不同。因为Bayes不依赖于最优化,所以A上有没有“距离”并不重要。但是Bayes需要一个在集合A上的prior P,也就是集合A上的一个测度(叫测度不叫分布函数是为了和集合A的元素区分开来,不然很容易混淆,虽然测度本身和分布函数没有太大区别)作为我们对于所有可能的X的分布函数的先验知识。而这就要求A上必须要有拓扑,也就是要能定义“开集”。尽管这个测度P被称为random measure(因为它是一个过程,是建立在分布函数上的分布),但是除了support是A这一点外,P和一般的prior并没有太大的区别。回到我们之前说的normal model,在那个例子里A的元素是所有的正态分布N(\mu,\sigma^2),在\mu,\sigma上加的任何prior也可以被看成是random measure — 只不过一般的random measure的测度都定义在离散点集上,比如经典的Poisson random measure。而在Normal model里的概率可以认为是定义在单点上。离散点集之间的关系比二维欧式空间里的点之间的关系要复杂很多 (比如距离就很难定义),这也就是为什么random measure要复杂的多。所以对于non-parametric来说,唯一的问题就在于A的结构太复杂,不能被直接嵌入欧式空间中来利用欧式空间上已有的各种测度。因此想要在集合A上建立一个概率空间,并找到对应的\sigma-field就会遇到麻烦—A上没有拓扑,没有定义开集,就更无从谈起Borel集了。而这也就是Ferguson这篇paper的主要目的,如何在这么大的A集合上建立一个足够好的测度空间。足够好就意味着这个测度要能覆盖到整个集合A,其次就是posterior要足够好算,要conjugate。前者要求定义这个测度的时候,要顾及到全部分布函数,后者要求定义的方法必须和有限的情形联系起来,是有限情形的推广(这样才能保证posterior好算)。在下一篇里我就会介绍Ferguson是如何构造这个测度P的。

思考:1.证明非参问题里的集合A无法被嵌入有限维欧式空间中?也即证明对于任意n, 总能找到一个从A到[0,1]^n上一个Lebesgue测度不为0的集合的满射。

2.证明任何连续函数都可以被可列个Gaussian density的加权和逼近(即证明文中的离散化公式,然后对于\forall\epsilon>0, 总存在一个h_0,使得当h<h_0|f(x)-\sum N(x;kh,h)\cdot hf(kh)|<\epsilon, x\in R