某师兄在临死(毕业)前,挖了个坑希望给他闲得XX疼的时候想出来的三硬币问题给出EM、Gibbs Sampling、VBEM、loop bp以及VBEP各种解法。现在他已经成为一名出色的sql工程师了。我算是给出这个Expectation Propagation解法,告慰他的在天之灵,少传播一点负能量。

引言

学习变分期望传播,好的材料还是看巨硬的这个网站,上面列了一大堆东西,反正我是没看完,我的印象中这也确实是他们的人搞出来的算法。他们还有Infer.net这种黑科技。

回到正题,EP可以看做一种loopy belief propagation的推广,他和VB不一样的地方在于变分推断中,我们是极小化\(KL(q||p)\),而EP则是直接极小化\(KL(p || Q)\)。正片文章的推导过程参(chao)考(xi) 这篇论文

三硬币

先来回顾一下这个问题,最初的问题是,三枚均匀程度未知硬币A、B和C,先投A如果是正面就投B,如果投反面就接着投C,我们只能看到最后B和C的结果,最终得到一个诸如正反正反正正反的这种序列,然后我们想知道每一个硬币的均匀程度以及每一次投的那个硬币。

这其实是一个很典型的带有隐变量的问题,长得和LDA也很像,然后将这个问题再推广一下就变成A是有K面的色子,B、C等都是有T面的色子,画出下面这个生成模型。

生成模型

生成模型

\(\alpha\)为超参 每一次色子A以\(\pi\)的概率生成\(Z\),然后\(Z\)又会选择剩下的几个塞子中的一个,最后选择\(\theta_{kt}\)的概率来生成展示的面\(w\)。整个模型挺简单的。

求解

Generative Aspect Model

这个模型求解的难处在于Z未知,所以求后验分布的时候需要遍历每一种Z的情况,这显然是不可能的,变分推断中,就是把\(\pi\)\(Z\)之间的关系切断,然后进行近似。其实换个贝叶斯一点的观点,\(Z\)我们观测不到,所以每一种\(Z\)都是出现了的,只不过他都是以一定的概率存在的,那么相关联的概率分布\(\pi\)的每一种结果也都出现了,所以我们干脆就放弃仔细刻画\(Z\)转而直接去搞定\(\pi\)\(Z\)两个变量分布结合以后的表现。也就是放弃topic,转向了aspect

具体过程

其实也不具体,因为打起来太麻烦。

先写出联合分布 \[\begin{split} p(w, \pi|\alpha) &= P(\pi|\alpha)P(w|\pi) \\ &= p(\pi | \alpha) \prod_{n = 1}^N p(w_n | \pi) \\ & = p(\pi| \alpha) \prod_{n = 1}^N(\sum_k \pi_k \theta_{ktn}) \\ & = p(\pi | \alpha) \prod_t(\sum_k \pi_k \theta_{kt})^{n_t} \end{split}\]

这里就把Z放弃,然后直接用\(\pi\)这个\(1 * K\)K的概率向量和\(\theta\)这个\(K * T\)的矩阵的乘积代替。然后对于这个求不出来的\(p(w | \pi)\),与变分推断相似的我们使用分布\(Z(\pi)\)代替。

即使用多个\(z_t(\pi) = s_t \prod_k \pi^{\beta_{tk}}_k\)代替\(p(w|\pi) = \sum_t \pi \theta_{tw}\)

然后整个分布就可以展开成

\[\begin{split} p(w,\pi|\alpha) & \sim p(\pi|\alpha)\prod_t(s_t \prod_k \pi_k^{\beta_{tk}})^{n_t} \\ & = \frac{\Gamma(\sum_k \alpha_k)}{\prod_k \Gamma(\alpha_k)}\prod_t s_t \prod_k \theta_k^{\alpha_k - 1 + \sum_t \beta_{tk}n_t} \end{split}\]

然后提出来\(\gamma_k = \alpha_k + \sum_t \eta_{tk}n_t\)就能得到

\[p(w,\pi|\alpha) \sim q(\theta) = \frac{\Gamma(\sum_k \alpha_k)}{\prod_k \Gamma(\alpha_k)}\prod_t s_t \prod_k \theta_k^{\gamma_k- 1 }\]

得到了没有归一化的\(q(\pi)\)之后,我们尝试归一化

\[q(\pi) = \frac{\frac{\Gamma(\sum_k \alpha_k)}{\prod_k \Gamma(\alpha_k)}\prod_t s_t \prod_k \theta_k^{\gamma_k- 1 }}{\int \frac{\Gamma(\sum_k \alpha_k)}{\prod_k \Gamma(\alpha_k)}\prod_t s_t \prod_k \theta_k^{\gamma_k- 1 } d\pi}\]

利用gamma函数的性质以及狄利克雷分布的性质,可以消掉积分号最后得到\(\frac{\Gamma(\sum_k \gamma_k)}{\prod_k \Gamma(\gamma_k)}\prod_k \theta_k^{\gamma_k- 1 }\)

我们想要利用moment matching的方法优化的话,就要一步一步的走,像gibbs sampling一样,去掉其中的一个,然后接近真实分布,然后再去掉另一个,再走下一步,所以,还需要去掉某一个\(t\)后的形式。 即 \[q^{\backslash t}(\pi) = \frac{\Gamma(\sum_k(\gamma_k - \beta_{tk}))}{\prod_k \Gamma(\gamma_k - \beta_{tk})}\prod_k \theta_k^{\gamma^{\backslash t}_k - 1}\]

其中 \(\gamma^{\backslash t} = \gamma_k - \beta_{tk}\)

然后把这个塞回到之前的\(p(w|\pi)\)里面,然后在归一化就得到完整的后验

\[\hat{q}{\pi} = \frac{\sum_k \gamma_k^{\backslash t}}{\sum_k \theta_{kt} \gamma_k^{\backslash t}}(\sum_k \pi_k \theta_{kt})\frac{\Gamma(\sum_k \gamma_k^{\backslash t})}{\prod_k \Gamma(\gamma_k^{\backslash t})}\prod_k \pi_k^{\gamma_k^{\backslash w } - 1}\]

然后就进入moment maching的阶段了 先求出来近似后验的一阶距为

\[E_{\hat{q}}[\pi_k] = \frac{\sum_k \gamma_k^{\backslash t}}{\sum_k \theta_{kt} \gamma_k^{\backslash t}}\frac{\gamma_k^{\backslash t}}{\sum_k \gamma_k^{\backslash t}}\frac{\theta_{kt} + \sum_k \theta_{kt}\gamma_k^{\backslash t}}{\sum_k \gamma_k^{\backslash t} + 1}\]

相似的计算二阶矩为

\[ E_{\hat{q}}[\pi_k^2] = \frac{\sum_k \gamma_k^{\backslash t}}{\sum_k \theta_{kt} \gamma_k^{\backslash t}}\frac{\gamma_k^{\backslash t}}{\sum_k \gamma_k^{\backslash t}} \frac{\gamma_k^{\backslash t} + 1}{\sum_k \gamma_k^{\backslash t} + 1}\frac{2\theta_{kt} + \sum_k \theta_{kt}\gamma_k^{\backslash t}}{\sum_k \gamma_k^{\backslash t} + 2} \]

由于\(\theta\)服从狄利克雷分布,则其一阶距应该为\(\frac{\gamma_k}{\sum_k \gamma_k}\),二阶矩应该为\(\frac{1}{k} \sum_k \frac{\gamma_k(\gamma_k + 1)}{\sum_k \gamma_k (\gamma_k + 1)}\)。两种形式对应起来,就能真正的参数应该为

\[\gamma_k = \frac{\sum_k (E_{\hat{q}}[\pi_k^2] - E_{\hat{q}}[\pi_k])}{\sum_k (E_{\hat{q}}[\pi_k]^2 - E_{\hat{q}}[\pi^2_k])}E_{\hat{q}[\pi_k]}\]

求得了真实的\(\gamma\)之后,就可以带回狄利克雷分布,由于我们假设了分布\(z\)所以使用分布\(z\)\(\hat{q}\)应该与\(Dir(\gamma)\)相等。

这样就可以解出了近似的参数

\[\beta_{tk} = \gamma_k - \gamma_k^{\backslash w}\] \[s_w = \frac{\sum_k \theta_{kt} \gamma_k^{\backslash t}}{\sum_k \gamma_k^{\backslash t}} \frac{\Gamma(\Gamma(\gamma_k^{\backslash t}))}{\prod_k \Gamma(\gamma_k)}\frac{\prod_k \Gamma(\gamma_k^{\backslash t})}{\Gamma(\sum_k \gamma_k^{\backslash t})}\]

然后不断迭代更新的那个t,直到收敛。

其实从求解过程就能看到,EP依赖于共轭分布的性质,这样才能保证后验分布的形式已知,然后依赖于指数分布家族的性质,分布可以用一阶距、二阶矩准确描述。




X