edfward

MCMC,LDA,文本建模 (二)

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

# init topic labels for every word randomly
for 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] += 1

while 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_distributiondoc_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 stepsizes
accepted  = 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 意味着

简单的 Bayes’ rule:

假设 alpha 为 Uniform distribution(待会儿试着解释这个假设),那么等价于

回到笔记里的 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的讨论