之前已经见过了变分推断这种求近似解的手段,它似乎为我们对复杂模型提供了一个漂亮的解法,但是我们也发现变分推断求解过程比较复杂,又比较依赖指数家族共轭分布(虽然非指数家族或者非共轭可解,但是手段更加复杂),MC(蒙特卡洛方法)就提供了另一种手段,其中在概率图模型上常用的就是MCMC(蒙特卡洛马尔科夫)方法,用一句话描述MCMC,就是对于在随机向量\(X\)状态空间上的随机过程\({X(1),...,X(t)}\),如果给出适当的概率转移矩阵并作用于\(X\),就可以构造出马尔可夫链使\(X(t)\)的分布收敛到其真实后验分布.因此,这一类算法称为MCMC算法。其中比较常见的就是Gibbs采样算法。

这里有一些有意思的前置知识。我们常说统计学分为频率学派与贝叶斯学派,频率学派认为我们要求的变量其实是个我们所未知的定值,如MLE(极大似然),MAP(极大后验)等方法其实是求了在其定义下的最优值,MLE是使似然函数最大的值,而MAP是使后验函数最大的值,但是贝叶斯学派认为这些值本身就是个分布,此时如果还是用MLE中的概率最大的单一值的话,就是损失掉很多信息,所以取期望是个更好的选项,回忆我们在面对复杂的生成模型时写出的联合分布往往只有参数节点和观察到的节点,希望将其中的变量都给积分积掉(其实本质上就是使用其期望),而且要注意这里的积分其实和随机变量本身是离散的还是连续的没有关系。比如对于观察到的\(\mathcal{X}\)和我们感兴趣的\(p(y|\pi)\),我们往往使用

\[p(y|x) = \int p(y|\pi)p(\pi|\mathcal{X})\]

明确了这一点,我们采用MC方法的合理性就明了了,很多积分是难求甚至是不可求的,而MC方法提供了一种通过采样接近的方法

\[\int f(x)p(x)dx = E_{p(x)}[f(x)]\simeq \frac{1}{n}\sum_{i=1}^{n}f(x_i)\]

即通过多次采样取期望的方式估计积分值。

好啦,这里我们已经有了求解积分的方法,剩下的就是对分布\(p(x)\)采样了,面对一个复杂的生成模型,怎样才能采样出对应的分布,MCMC方法就解决了这个问题,首先重要的一点是马尔科夫链细致平稳分布,即对于满足一定性质的马氏链(我们常见的概率图模型基本都满足这个性质,详细的可以看这篇介绍)具有转移概率矩阵P,进过一定次数的迭代后,会趋于一个稳定分布。所以,对于\(x\)来看,如果我们找到一个合适的概率转移矩阵,那么经过一定次数的迭代,我们就能获得合理的\(p(x)\)生成\(x\)的采样。gibbs采样就是利用这种方法求解。

gibbs采样在实际应用的推导中显得较为简单,它假设对于\(p(x)\)\(p(y)\)两个难求的分布,如果我们可以通过\(p(x)=\int p(x,y)dy\)来求解,即从\(p(x|y=y_{i-1})\)中生成\(x_i\),同理,对于\(y\)也有\(p(y)=\int p(x,y)dx\),即\(p(y|x = x_i)\)生成\(y_i\),则这种方式可以生成一个gibbs序列,在一定的迭代次数之后就会收敛到固定的分布,此时采样就可以求得我们想要的积分。 gibbs采样的形式比较简单,对于要采样的\(Z\),有\(z=<z_ 0,z_1,...,z_k>\),其中\(k>1\),则我们可以

初始化\(z_{(0)}:=<z_1^{(0)},...,z_j^{(0)}>\)
  for \(t = 1\) to T do
    for \(i = 1\) to k do
      \(z_i^{(t+1)}\sim P(Z_i|z_1^{(t+1)},...,z_{i-1}^{(t+1)},z_{i+1}^{(t)},...,z_k^{(t)})\)
    end for
  end for

\[P(Z_i|z_1^{(t+1)},...,z_{i-1}^{(t+1)},z_{i+1}^{(t)},...,z_k^{(t)})=\frac{P(z_1^{(t+1)},...,z_{i-1}^{(t+1)},z_{i+1}^{(t)},...,z_k^{(t)})}{P(Z_i|z_1^{(t+1)},...,z_k^{(t)})}\]

之前已经看过LDA的变分推断的推导,这里用更复杂一点的B-LDA(论文戳这里)为例子,先看生成模型(详细一点的读书笔记可以看这里)。

与LDA相比,针对twitter数据,论文进行了一些新的假设首先是简化LDA中一个文章多个主题的情况,在twitter的短文中一条推文只有一个主题,其次加入背景词,认为twitter中有一些常见的词并不是由主题生成而是由背景词生成,另一方面认为用户的行文本身其实也是由主题决定的,为此还有一个生成行为的多项分布。

具体的生成过程为:

对于\(1\leq t \leq T\),生成\(\phi_t \sim Dir(\beta)\),\(\psi_t \sim Dir(\eta)\) 生成\(\phi' \sim Dir(\beta')\),\(\varphi\sim Dir(\gamma)\) 对于每一个用户从\(u=1,...,U\) 利用\(\theta_u\sim Dir(\alpha)\)生成每个用户的主体分布 对于该用户的第n个推文,\(n=1,...,N_u\) 从该用户的\(\theta_u\)生成该推文主题\(z_{u,n}\) 对于每一个单词\(l=1,...,L_{u,n}\) 生成\(y_{u,n,l} \sim Bernoulli(\varphi)\),来决定词是由背景词还是主题生成 当\(y_{u,n,l}=1\)的时候,从主题中生成单词,即\(w_{u,n,l} \sim \phi_{z_u,n}\),当\(y_{u,n,l}=0\)的时候,从背景词中生成单词,即\(w_{u,n,l}\sim \phi'\)。 生成用户行为\(b_{u,n} \sim \psi_{z_u,n}\)

这里可以看到我们想要采样求解的隐变量决定词生成来源的\(y\)和主题\(z\),首先先写出联合分布有

\[\begin{split}&p(Z,Y,B,W|\alpha, \beta, \beta',\gamma,\eta) \\ =&p(W|Z,Y,\beta, \beta')p(Y|\gamma)p(B|Z,\eta)p(Z|\alpha)\end{split}\]

然后依次展开

对于\(p(Z|\alpha)\)

\[\begin{split}p(Z|\alpha) &= \int p(Z|\theta)p(\theta|\alpha) \\ &= \int \big(\prod_u^U\prod_{t}^T\theta^{n_u^t}_{t,u}\big)\big(\prod_{u}^U\frac{1}{\Delta(\alpha)}\prod_t^T\theta_{t,u}^{\alpha_t-1}\big)d\theta \\ &= \int\big(\prod_{u}^U\frac{1}{\Delta(\alpha)}\prod_t^T\theta_{t,u}^{n_u^t+\alpha_t-1}\big)d\theta \\ &= \int \prod_{u}^U\frac{\Delta(n_u^t+\alpha)}{\Delta(\alpha)}\frac{1}{\Delta(n_u^t+\alpha)}\prod_t^T\theta_{t,u}^{n_u^t+\alpha_t-1} d\theta \\ &= \prod_{u}^U\frac{\Delta(n_u^t+\alpha)}{\Delta(\alpha)} \end{split}\]

这里我们发现由于采用了共轭的先验,所以对于参数很容易就能通过积分消掉,既能减少采样的复杂度,又能简化式子的形式。我们观察这个式子的内容,很容易就能发现,其形式的规律。通过相似的方法还能求解出其他的带超参的。

\[p(B|Z,\eta) = \int p(B|Z,\psi)p(\psi|\eta)d\psi = \prod_{t=1}^T \frac{\Delta(n_b^t+\eta)}{\Delta(\eta)}\]

\[p(Y|\gamma)=\frac{\Delta(n_{(.)}^y+\gamma)}{\Delta(\gamma)}\]

这里的\(n_{(.)}^y\)实际上是同时代表了\(y=1\)\(y=0\)的情况。

现在只还剩\(p(W|Z,Y,\beta,\beta')\),这其实是模型核心的部分,通过指示变量\(y\)指示从背景词分布中生成,还是从主题中生成,根据之前的规律,其实我们期望积分积掉\(\phi\)\(\phi'\)

\[\begin{split}p(W|Z,Y,\beta,\beta') & = \int_\phi \int_{\phi'} p(W|Z,Y,\beta, \beta')p(\phi|\beta)p(\phi'|\beta')d\phi d\phi' \\ &=\int_\phi \int_{\phi'} (\prod_w^W\phi'^{n_{y=0}^w}_w\prod_t^T\phi_t^{n^w_{t,y=1}}) (\prod_t^T\frac{1}{\Delta(\beta)}\prod_w^W\phi_{t,w}^{\beta_{w}-1})(\prod_w^W\phi'^{\beta'_w-1}_w)d\phi d\phi' \\ &= \frac{\Delta(n_{y=0}^w+\beta')}{\Delta(\beta')}\prod_{t=1}^T\frac{\Delta(n_{t,y=1}^w+\beta)}{\Delta(\beta)} \end{split}\]

有了以上几个式子,我们就可以写出来最初的式子为

\[\begin{split}&p(Z,Y,B,W|\alpha, \beta, \beta',\gamma,\eta) \\ &=p(W|Z,Y,\beta, \beta')p(Y|\gamma)p(B|Z,\eta)p(Z|\alpha)\\ &=\frac{\Delta(n_{(.)}^y+\gamma)}{\Delta(\gamma)}\frac{\Delta(n_{y=0}^w+\beta')}{\Delta(\beta')}\prod_{t=1}^T\frac{\Delta(n_{t,y=1}^w+\beta)}{\Delta(\beta)}\\ &\prod_{t=1}^T \frac{\Delta(n_b^t+\eta)}{\Delta(\eta)}\prod_{u}^U\frac{\Delta(n_u^t+\alpha)}{\Delta(\alpha)} \end{split}\]

下一步就是求要采样的条件分布。联合分布中\(B,W\)都是观察到的变量,并且我们在采样的过程中已经积分掉了参数,剩下的就是针对\(Z、Y\)的采样对于每一个Z,由于每一个用户的每一条twitter都有自己的\(Z\),所以在求每一个条件分布的时候,应该对应每一个\(Z_{u,n}\)即第u个用户的第n个推文,如果我们要剥离这条推文,那么与这条推文相关的用户行为即\(B_{u,n}\)也应该被剥离,这里用\(c\)代表\(\{u,n\}\),于是有

\[\begin{split}&p(z_c|Z_{\lnot c}, W, Y, B)\\ &\propto\frac{p(Z,Y,B,W|\alpha, \beta, \beta',\gamma,\eta)}{p(Z_{\lnot c},Y,B_{\lnot c},W|\alpha, \beta, \beta',\gamma,\eta)} \\ &=\frac{\Delta(n_{z,y=1}^w+\beta)}{\Delta(n_{z,y=1,\lnot c}^w+\beta)} \frac{\Delta(n^b_z+\eta)}{\Delta(n^b_{z,\lnot c}+\eta)} \frac{\Delta(n^u_z+\alpha)}{\Delta(n_{z,\lnot c}^u+\alpha)} \end{split}\]

将伽马函数消掉就可以得到最后的迭代式,只要注意到\(\Gamma(x + 1) = x\Gamma(x)\),就能利用各个统计量\(n\)之间的差消掉伽马函数,得到迭代式。

同理,对于指示标签\(y\),考虑到对于每一个词都是由\(y_{u,n,l}\)指示的,所以应该剥离掉\(y_{u,n,l}\),这里用\(d\)表示\({u,n,l}\)。由于剥离了指示标签\(y_d\),其实在生成模型中与其相关的项目也是可以直接看出来的,即由\(y\)指示的部分和生成的部分。

\[\begin{split}&p(y_d|Y_{\lnot d},Z,W,B) \\ &=\frac{\Delta(n_{(.)}^y+\gamma)}{\Delta(n^y_{\lnot d}+\gamma)}\frac{\Delta(n^w_{y=0}+\beta')}{\Delta(n^w_{y=0,\lnot d}+\beta')}\prod_{t=1}^U\frac{\Delta(n^w_{t,y=1}+\beta)}{\Delta(n^w_{t,y=1,\lnot d}+\beta)} \end{split}\]

这里其实包含了\(y=1\)\(y=0\)两种情况,分别会导致对后面两项不同的结果,而且剥离一个\(y_d\)其实指示对统计量造成了1的影响,所以得到的形式很简单。

这样我们就得到了所有需要采样的变量的表达式,然后我们先对未知的部分随机初始化,然后通过采样的参数生成具体的分布,最后多次采样得到值就可以了,有一些要注意的是,首先由于随机初始化要经过一定的迭代次数才能达到稳态,所以最初的阶段称为build-in,这段时间一般对结果采样,另一方面,为了获得结果,可以每次都对变量采样,也可以隔一定的次数采样,由于其实采样的结果是分布而不是值(其实贝叶斯就是这样认为的),所以根据目的不同,其实针对结果也有很多处理的方式。而且,对于部分隐变量已知的情况(比如部分\(z\)已知时),可以在采样中不改变标签来利用已知信息。

这里我们可以发现,相比变分推断这种方法,采样的方法从形式上来看其实更加简单,但是难以直接判断收敛,但是对于分布如果是非指数家族的时候始终很实用的方法




X