MCMC，LDA，文本建模，来点干货(一)

Fri 27 September 2013 | in Study | tags: probabilitybayesmachine_learningstatisticsldamarkovmcmcevernote

很久没写 blog 了，PyMoTW 也没想好下一个对象（不过这次有涉及到一些 matplotlib…），刚好最近啃完一篇关于文本建模的科普「LDA数学八卦」而且啃得痛苦无比，也做了很多的笔记(扫描到 Evernote 上了)，这次就先聊聊它们吧。

先梳理一下做的笔记，大致结构与原文「LDA数学八卦」相仿，感兴趣的话请务必翻翻那篇文章：

1. 讲的是 Gamma 函数的推导和 Beta 函数的引入
2. 写了 Gamma 分布以及 Binomial 分布与 Poisson 分布的关系（略去了从 Binomial 到 Gamma 的推导因为涉及到了一个比较复杂的等式）
3. 从一个例子入手着重讲了 Beta 分布。这部分强烈推荐看原文，魔鬼和骰子的故事讲的很生动
4. Beta-Binomial 共轭，是前面例子的延伸，对理解之后的 Dirichlet-Multinomial 共轭和更后面的 LDA 建模有很大的帮助
5. 推广到 D-M 共轭，同时介绍了这两个分布的参数估计
6. 噔噔噔噔，进入 Markov Chain 部分~讲了 Markov Chain 的平稳分布收敛条件和一些性质，是后文 Markov Chain Monte Carlo 算法的基础
7. MCMC。平稳分布的判定条件(detailed balance condition)、如何人工满足这个条件（α 接受概率的引入）以及以此为基础的 MCMC 采样算法（最简单的 Metropolis-Hastings 采样算法）。这部分难度比之前大了一些（个人感觉）
8. M-H 算法衍生的大名鼎鼎的 Gibbs Sampling，以及从二维到高维的推广
9. 进入文本建模部分。Unigram Model，涉及 D-M 共轭和参数估计
10. PLSA (Probabilistic Latent Semantic Analysis) 简述
11. 比 PLSA 更进一步（或者说更 Bayesian）的 LDA (Latent Dirichlet Allocation)。看了原文不下4，5遍才算是勉强理解。重点在于理解文档生成的方式以及 D-M 共轭是怎么出现的
12. 有了模型之后利用 Gibbs Sampling 进行参数计算和文档统计，这部分数学稍重

可以发现1-5是概率知识的普及，6-8是 MCMC 的介绍，从9开始才是真正的文本建模部分。

之后两张纸是「漫谈 HMM：Forward-Backward Algorithm」以及「漫谈 HMM：Definition」（都来自目前在 MIT 的pluskid 大牛的博文，个人偶像之一）。HMM 经常忘，刚好笔记也不是太长就附在后面了。有兴趣的话可以看看原文再看看笔记，希望有帮助。

结果再唠叨几句。在实验室做事情花的时间也不少了，但是看 paper 始终会觉得这儿不懂那儿不懂，然后翻回之前提到这个算法的 paper 再重复这个过程。我觉得相比起在 paper 中间那么点儿地方去了解一个算法不如翻回最早提到该算法的 paper 或者是别人整理出来的科普，拿出纸笔去算去写，才能有更好的领悟吧。像刚刚结束的组会提到了 wordnet random walk，当时不明白为什么一定会 converge，后来翻了翻笔记发现它就是个典型的 Markov random walk，根据笔记 part 6 里提到的如果马氏链满足平稳分布收敛条件（非周期且各个状态连通），最终一定会收敛到平稳分布去；还提到了 explicit semantic analysis，如果了解了 PLSA 那么对 ESA 的理解就更没问题了。发现刚知道的知识在实际问题中可以如此应用的时候便觉得花这么些功夫还是值得的。

接下来分析几篇相关的文吧，也算强化自己对这些知识的理解。

「Gibbs sampler in various languages (revisited)」: 比较了各个语言的 Gibbs Sampler 的效率（当然 C 完胜，Python 在 PyPy 的帮助下只慢了三倍…）。值得关注的是一个典型的 Gibbs Sampler 是怎么写的（Python 只用了13行），功夫其实是花在对要采样的分布进行分析上的。由于本人微积分都还给老师了所以只好请准备考研的室友把那个分布的 conditional probability 算了一下，了不起的是那么复杂一个积分他居然还算对了…赞。

「如何生成随机数（上）」: Warning: 该文没有(下)…依然来自于刚刚提到的偶像 pluskid。浅显易懂的讲了随机数的生成：如何从 Uniform Distribution 产生的随机数映射到更复杂的分布，以及比较简单的一个采样算法 Inverse Transform Sampling Method，当然也少不了方便直观展示的 python 代码。虽然烂尾了但还是好评。

另外还想谈的一篇文由于要说的太多，下次再讲吧（里面遇到的问题想了好久最后还是通过洗了一个澡才解决的…）。

那么，下回再见。

PS: 顺便把自己觉得写得最好的一页贴出来…

读了 plda (A parallel C++ implementation of fast Gibbs sampling of Latent Dirichlet Allocation) 的源码，不长，对理解 LDA 很有帮助，但还是要做好花时间的准备。其实核心算法并不复杂，见如下的 python 伪代码（基本上就是把上一章的算法用伪代码表示出来，就当复习了吧）：

# init topic labels for every word randomlyfor m in LDACorpus:    for w in m:        # note that same VOCAB has many word POSITIONS        topic[m][w] = random topic `k`        topic_word_distribution[k][w] += 1        doc_topic_distribution[m][k] += 1        # word_topic_distribution[w][k] += 1while interating:    for m in LDACorpus:        for w in m:            # generate topic distribution for w            distribution = {}            for k in Topics:                # calculate P(z_k | W_^w, z_^w)                diff = -1 if topic[m][w] == k else 0                d_factor = doc_topic_distribution[m][k] + diff                t_factor = topic_word_distribution[k][w] + diff                distribution[k] = (d_factor + alpha) * (t_factor + beta) / \                          (sum(topic_word_distribution[k]) + beta * Vocab_size)            # sample from distribution            k = sample(distribution)            topic[m][w] = k            # update model            topic_word_distribution[k][w] += 1            doc_topic_distribution[m][k] += 1    # check convergence and post process

topic_word_distribution 和 doc_topic_distribution 都非常直观，记录每个文档的 topic 分布和每个 topic 的词汇分布，其中topic_word_distribution 作为最后的 model 返回即可。初始化就不用说了，随机给每篇文档每个词（注意：包括重复词，这里可以看作每个『位置』）取一个 topic，并更新两个 model。接下来迭代的过程就是重新遍历每个词，利用 Gibbs Sampling 来估计这个词的 topic 分布然后取样、更新 model（实际上就是 condition 在除开这个词之后的语料计算 P(z_k | w_~i, z_~i)）。distribution[k] 的计算涉及到上篇文章笔记里做的参数估计，也是 LDA 整个部分数学要求最高的地方（「LDA 八卦」作者语）。最后检查 convergence 并简单的平均一下（似乎就是将 model 参数除以 [总迭代 - burn_in] 取个平均， 这部分看得不太仔细，如理解错误敬请提出）。

plda 的代码是我见过的写得最好看的 C++ 代码之一（对不起其实我代码看得少…），感觉上可以直接投入工业强度级别的使用（无责任瞎猜）。核心算法就上面几十行，plda 包装一下再支持个并行最后变成2000行左右的包，可读性和健壮性都挺好，阅读的过程中也体会到一个算法是怎样解构、重现的，哪些地方需要什么样的数据结构又需要什么样的异常处理，对基本上没用 C++ 写过大型工程的我来说有很大的帮助。

好，那再回到上次结尾提到的一篇文章，讲的是 Metropolis-Hastings algorithm （以下简称 MH algorithm）。笔记里的 MH 算法是这样的

而这篇教程里讲的 MH 部分是这样的：

# initial guess for alpha as array.guess = [3.0]# Prepare storing MCMC chain as array of arrays.A = [guess]# define stepsize of MCMC.stepsizes = [0.005]  # array of stepsizesaccepted  = 0.0# Metropolis-Hastings with 10,000 iterations.for n in range(10000):    old_alpha  = A[len(A)-1]  # old parameter value as array    old_loglik = evaluateLogLikelihood(old_alpha, D, N, M_min,                    M_max)    # Suggest new candidate from Gaussian proposal distribution.    new_alpha = numpy.zeros([len(old_alpha)])    for i in range(len(old_alpha)):        # Use stepsize provided for every dimension.        new_alpha[i] = random.gauss(old_alpha[i], stepsizes[i])    new_loglik = evaluateLogLikelihood(new_alpha, D, N, M_min,                    M_max)    # Accept new candidate in Monte-Carlo fashing.    if (new_loglik > old_loglik):        A.append(new_alpha)        accepted = accepted + 1.0  # monitor acceptance    else:        u = random.uniform(0.0,1.0)        if (u < math.exp(new_loglik - old_loglik)):            A.append(new_alpha)            accepted = accepted + 1.0  # monitor acceptance        else:            A.append(old_alpha)print "Acceptance rate = "+str(accepted/10000.0)

在进入问题之前先简单介绍一下，这个 snippet 是为了估计一个参数 α。先用已知的 α（值为2.35）代入一个特定的模型生成1000000个数据，再用这些数据去估计 α（MH 的结果是2.3507正负0.0015，相当准确了）。整个过程中模型（即概率分布）是已知的，未知的只是 α 一个参数而已。

当时看这段 code 最奇怪的是中间偏后接受采样的部分：为什么和 likelihood 有关？那个 if-else 的逻辑是什么意思？如果你觉得对 MH 理解比较深的话可以思考一下，至少我从刚刚的笔记出发推出这个算法是花了点时间的（而且是想半天想不通无奈去洗澡才突然灵感爆发）。

揭晓答案。

new_loglik > old_loglik 意味着

<img src="http://latex.codecogs.com/gif.latex?P(\vec{M}|\alpha_{new})%3EP(\vec{M}|\alpha_{old})" title="P(\vec{M}|\alpha_{new})>P(\vec{M}|\alpha_{old})" >="" <="" p="" style="background-color: transparent; border: 0px; margin: 0px; outline: 0px; padding: 0px; vertical-align: baseline;">

简单的 Bayes' rule:

<img src="http://latex.codecogs.com/gif.latex?\frac{P(\alpha_{new}|\vec{M})P(\vec{M})}{P(\alpha_{new})}%3E\frac{P(\alpha_{old}|\vec{M})P(\vec{M})}{P(\alpha_{old})}" title="\frac{P(\alpha_{new}|\vec{M})P(\vec{M})}{P(\alpha_{new})}>\frac{P(\alpha_{old}|\vec{M})P(\vec{M})}{P(\alpha_{old})}" >="" <="" p="" style="background-color: transparent; border: 0px; margin: 0px; outline: 0px; padding: 0px; vertical-align: baseline;">

假设 alpha 为 Uniform distribution（待会儿试着解释这个假设），那么等价于

<img src="http://latex.codecogs.com/gif.latex?P(\alpha_{new}|\vec{M})%3EP(\alpha_{old}|\vec{M})" title="P(\alpha_{new}|\vec{M})>P(\alpha_{old}|\vec{M})" >="" <="" p="" style="background-color: transparent; border: 0px; margin: 0px; outline: 0px; padding: 0px; vertical-align: baseline;">

回到笔记里的 MH 算法，为防止歧义用 A 表示原先的 α，再 conditionalize 到观察数据上：

假设 transition kernel 是 Gaussian 分布（代码中 new_alpha[i] 的赋值也能看出）

那么转移应该是对称的

此时

因而从均匀分布中取样的 u 小于 A(a_old, a_new) 恒成立，而代码中满足这个 if 就可以直接接受这个 sample！否则再去比较 u 与 A(a_old, a_new)，也就是 math.exp(new_loglik - old_loglik)：同样的步骤可以推导出这两个家伙是等价的。

这段代码实现的 MH 算法搞清楚了，回到之前那个假设：为什么 alpha 是均匀分布呢？答案是：我也不知道...只是做这个假设的话可以满足 proposal distribution（也就是那个 Gaussian），所以这个假设应该是没问题的...吧？（也有可能我把整个过程都理解错了...求指导！）

OK，这个系列的第二篇也算结束了，列一下该系列下一篇要谈的东西和要做的功课吧，备忘。

1. Probabilistic Topic Models, by Mark Steyvers from University of California, Irvine
2. Approximate Inference: Markov Chain Monte Carlo, slides of Probabilistic Graphical Models by Eric Xing
3. 微博上关于LDA和PLSA的讨论

-EOF-

转自：http://blog.edfward.com/2013/10/mcmcldawen-ben-jian-mo-lai-dian-gan-huo-er.html