为什么能联想到这些词，因为在我们的认知下，这些词出现的频率比较高，换句话说，我们接触的这些词比较多，比较频繁，所以我们能在确定主题之后很快的从这些词中选择一个

这些是前话，现在来开始说一下工程上怎么实现吧，我自己是用python来写的，在这里跟大家分享一下

首先是LDA模型的类定义：

class LDAModel:
alpha = float      #超参数alpha
beta = float       #超参数beta
D = int            #文档数目
K = int            #主题个数
W = int            #词的个数
NumberOfIterations = int   #迭代次数
SaveStep = int     #存储的步数
Dictionary = object  #整个语料的词典
Z = object         # D * doc.size()大小的矩阵，Z[i][j]表示第i文档的第j个词背分配的主题
W = object         # D * doc.size()大小的矩阵， W[i][j]表示第i文档的第j个词
IDListSet = object  # D * doc.size()大小的矩阵， IDListSet[i][j]表示第i篇文档的第j个词在词典中的编号
nw = object        # W * K 大小的矩阵， nw[w][z]表示词w被分配到主题z的次数
nd = object        # D * K 大小的矩阵，nd[d][z]文档d中被分配为主题z的词的个数
nwsum = object     # K * 1 大小的向量，nwsum[z]表示主题z中包含的词的个数
ndsum = object     # D * 1 大小的向量，ndsum[d]表示文档d中包含的词的个数
theta = object     # D * K 大小的矩阵，p(z|d) = theta[d][z]
phi = object       # K * V 大小的矩阵，p(w|z) = phi[z][w]

具体我就不说这些成员的意思，注释上都有，

首先是构造函数，这里要说明的是，在工程上，alpha一般取50/K, beta一般取0.01,吉布斯抽样的迭代次数一般为1000次、

    def __init__(self, alpha, beta, NumberOfIterations, SaveStep, K):self.alpha = alphaself.beta = betaself.NumberOfIterations = NumberOfIterationsself.SaveStep = SaveStepself.K = K#初始化大小为K * 1的向量，初始值为0self.nwsum = ListUtil.Initial(self.K)

有一些列表工具类的方法我已经打包了，先列出来

#ListUtil.py

import string
def Normalize(list, smoother=0.0):"""对向量list进行归一化处理，得到每个元素出现的概率:param list: 向量:param smoother: 平滑值，缺省值为0; 为了防止0概率的出现"""sum = Sum(list)K = len(list)newlist = []if sum > 0:newlist = [float((item + smoother) / (sum + K * smoother)) for item in list]return newlistdef Sum(list):"""计算list中所有元素的和"""res = 0for item in list:res += itemreturn resdef Initial(size, data=0):"""生成一个大小为size, 所有元素都为data的列表:param size: 列表大小:param data: 列表元素"""list = []for i in xrange(size):list.append(data)return listdef InitialMat(M, N, data=0):"""初始化大小为M * N的矩阵，所有元素初始化为data:param M::param N::param data: 矩阵元素"""mat = []for i in xrange(M):row = Initial(N, data)mat.append(row)return matdef InitialEmptyMat(rows):"""初始化一个空的matrix:param rows:"""mat = []for i in xrange(rows):tmp = []   #代表每一个文档包含的词，初始化为空mat.append(tmp)return matdef toString(list):"""将list中的元素拼接成字符串方便用作文件操作:param list: 列表元素"""listStr = ""count = 0for ele in list:if type(ele) == int:eleStr = str(ele)elif type(ele) == float:#浮点数转换为字符串，保留8位小数eleStr = str("%.10f"%ele)elif type(ele) == str or type(ele) == unicode:eleStr = eleif count != len(list) - 1:eleStr += " "count += 1listStr += eleStrlistStr += "\n"return listStrdef StringToFloatList(SS):"""string 转换为float:param SS: 从文件中读取的字符串"""res = [string.atof(item) for item in SS.split(" ")]return resdef AssignList(LL):"""将LL中的值拷贝到另一个list中:param LL: 字符串"""newLL = []for ele in LL:newLL.append(ele)return newLLdef FindMax(LL):"""返回列表LL中最大的元素"""LL.sort()return LL[len(LL) - 1]

好，接着我们刚才的定义

现在可以开始初始化模型了，代码如下：

def ModelInit(self, filename):"""读取文档，文本预处理，构造词典，构造语料库"""Docs = LoadData.LoadDataFromFile(os.getcwd() + "/" + filename)self.D = len(Docs)print "Load ", self.D, " docs from the file"#读取停用词表StopWordList = LoadData.LoadStopWords()#对输入文本进行预处理：去标点符号，去停用词，词干化，然后每篇文档生成一个词的列表WordListSet = [Preprocess.PreprocessText(doc, StopWordList) for doc in Docs if type(doc) != unicode]#通过词表集构造词典self.Dictionary = Preprocess.ConstructDictionary(WordListSet)self.W = len(self.Dictionary)print "Total number of words is: ", self.Wprint "Begin to save the dictionary..."self.SaveDictionary()print "Done!!"#IDListSet 大小 D * doc.size()print "Begin to map the word to ID"self.IDListSet = []for wdl in WordListSet:IdList = Preprocess.Word2Id(wdl, self.Dictionary)self.IDListSet.append(IdList)print "Done!!"#ndsum[d] 文档d中包含的词的个数self.ndsum = ListUtil.Initial(self.D)#初始化一个 D * K的矩阵self.theta = ListUtil.InitialMat(self.D, self.K, 0.0)self.phi = ListUtil.InitialMat(self.K, self.W, 0.0)#nd[d][z] 文档d中被分配给主题z的词数self.nd = ListUtil.InitialMat(self.D, self.K, 0)#nw[w][z] 主题z中包含的词w的个数self.nw = ListUtil.InitialMat(self.W, self.K, 0)#Z[d][w] 文档d的第w个词的主题self.Z = []print "Begin to initialize the LDA model..."#初始化计数向量和计数矩阵self.RandomAssignTopic()print "Topic assignment done!!"

首先是从文件中读取文档，LoadData同样是我定义的工具类；然后用辅助类Preprocess去完成文本的预处理（包括去标点符号，去停用词，词干化，构造词典等等），初始化完成之后，再为每个词赋一个初始的topic，你可能要问，LDA中文档主题的选择应该要服从狄利克雷先验分布，但是为什么可以随机赋值，其实这要从马尔可夫链开始说起了，简单的来说，马尔可夫链就是对象的一系列状态的集合，并且对象的当前状态仅仅跟它的上一个状态有关，来看一个具体的例子

这个图反应的是一个人的收入跟它父母收入高低的联系，1代表下层，2代表中层，3代表上层；上面矩阵的第一行就表示，如果父代是上层阶级，那么他的子代有0.65的可能性仍然是上层阶级，有0.28的概率称为中层阶级，有0.07的概率成为下层阶级；上面的矩阵我们就定义为传递概率矩阵P，如果我们假设当代人属于这3个阶级的概率分别是x(0，1) x(0，2) x(0，3)，那么他们的下一代的阶级分布x(1) = x(0) * P，同样，再下一代的阶级分布x(2) = x(1) * P，其他的以此类推，下面来看一个例子：

我们假设初始代的阶级分布x(0) = [0.21 0.68 0.11]，那么通过结算可以得到以下结果

我们再换一个初始的阶级分布，假设现在x(0) = [0.75 0.15 0.1]，计算结果如下：

通过这两个例子，我们发现无论初始的状态怎么样，最后的阶级分布都能够收敛，并且收敛到同一个分布，这就是马尔可夫链的厉害之处，虽然我到现在仍然无法理解这其中的奥秘，但是把他应用到工程中还是可以的

马尔可夫把这个结论总结成了一个定理（不知道怎么打公式，直接上图了）：

这里的P就是我们上面所说的传递概率矩阵，最后的pi就是收敛之后的阶层分布，非周期性马氏链我也不太清楚是什么一回事儿，反正你要记住，我们日常生活中接触的大都是非抽泣性马氏链

正是因为马氏链这种收敛性，所以我们才能够在工程上为文档中的词随机分派主题，因为不管怎样，它最后都能够熟练到一个平稳分布，只是说收敛的快慢问题

现在来贴一下这两个工具类

#Preprocess.py
import string
import nltk
from gensim import corporadef PreprocessText(text, StopWordList):"""预处理一篇文本：剔除标点符号，词干化，去停用词:param text: 传入的文本，类型为字符串:param StopWordList: 停用词表"""WordList = DelPunctuation(text)StemmeredWordList = Stemmer(WordList)FilteredWordList = FilterStopWords(StemmeredWordList, StopWordList)return <span style="font-family: arial, 宋体, sans-serif;">FilteredWordList</span>
def DelPunctuation(text):"""剔除文本中的标点符号:param text:需要剔除标点符号的文本，类型为字符串return:返回文本中的词的序列"""delset = string.punctuation#将标点符号转换为空格newText = text.encode('utf8').translate(None, delset)#文本中的词的列表WordList = [word for word in newText.split(" ") if word != '' and word != ' ']return WordListdef FilterStopWords(WordList, StopWordList):"""返回去停用词后的词表:param WordList::param StopWordList:"""FilteredWordList = filter(lambda x: x.lower() not in StopWordList, WordList)return FilteredWordListdef Stemmer(WordList):"""对文档的词表进行词干化:param WordList:"""stemmer = nltk.LancasterStemmer()StemmeredWordList = [stemmer.stem(w) for w in WordList]return StemmeredWordListdef ConstructDictionary(WordListSet):"""根据输入文档集texts构造词典:rtype : object:param WordListSet: 文档集对应的词表，WordListSet[i]表示第i篇文档中的词"""print "Begin to construct the dictionary"res = corpora.Dictionary(WordListSet)print "Total number of words is: ", len(res)return resdef Word2Id(WordList, Dictionary):"""将词表转换为词典dictionary中的ID:param WordList:"""IDList = []for word in WordList:#遍历字典查找目标项for k, v in Dictionary.items():if v == word:IDList.append(k)return IDList

在文本与处理时，用到了nltk这个强大的自然语言处理的库，程序中使用其中的LancasterStemmer()进行词干化；然后也用到了gensim库，在这个类中，主要是用corpora来构造训练文档集的词典

在贴一下LoadData的代码

#LoadData.py
import os
import stringdef LoadDataFromFile(path):""":param path:短文本存放路径"""#转换为绝对路径fp = open(path, 'r')Docs = []for line in fp:#去掉结尾换行符ll = line.strip('\n').strip('\r')Docs.append(ll)fp.close()print "Done, load ", len(Docs), " docs from the file"return Docsdef LoadStopWords():"""从指定路径读取停用词表return:停用词列表"""path = os.getcwd()path += "/StopWords.txt"fp = open(path, 'r')#获取停用词列表StopWordsList = [line.strip('\n') for line in fp]fp.close()return StopWordsListdef LoadDictionary():"""从指定路径加载训练词典"""path = os.getcwd() + "/dictionary.txt"fp = open(path, 'r')Dictionary = dict()for line in fp:elements = line.strip('\n').split(" ")#词的idk = string.atoi(elements[0])#词本身v = elements[1]Dictionary[k] = vfp.close()return Dictionary

这个类我就不多解释了，学过python的小伙伴应该都能看懂，只是涉及文件操作的路径名你们可以自己diy，我用的是我自己电脑上的文件名

继续解析LDAModel类

初始化模型之后，我们就要开始训练模型参数了（后面的是重点哟）：

刚才忘记贴随机分派主题的代码了：

def RandomAssignTopic(self):"""随机为文档中的词分配主题更新计数向量ndsum, nwsum, 计数矩阵nd, nw的值"""for d in xrange(self.D):DocSize = len(self.IDListSet[d])row = ListUtil.Initial(DocSize)self.Z.append(row)for w in xrange(DocSize):#从主题编号0-K-1中随机抽取一个topic = Sample.UniSample(self.K)#获取词的IDwid = self.IDListSet[d][w]self.Z[d][w] = topic#被分派给topic的词w的数目自增1self.nw[wid][topic] += 1#文档d中被分配给主题topic的词的个数self.nd[d][topic] += 1#主题topic中包含的总的词数self.nwsum[topic] += 1self.ndsum[d] = DocSize

lda的训练过程主要就是吉布斯抽样的过程，具体的来说，吉布斯抽样就是将抽样的一个词w从当前的分布中抽出，然后通过抽出这个词之后的主题分布theta和词的分布phi，来计算这个词被分派到其他主题的概率，先上代码

    def sampling(self, d, w):
"""
Gibbs Sampling为当前词重新分配主题
:param d: 文档编号
:param w: 词在文档中的编号
"""
topic = self.Z[d][w]
#对应位置上的词的ID
wid = self.IDListSet[d][w]
self.nw[wid][topic] -= 1
self.nd[d][topic] -= 1
self.nwsum[topic] -= 1
self.ndsum[d] -= 1
#p为马尔可夫链传递概率，p[z]表示当前词被分配到主题z的概率
p = self.ComputeTransProb(d, w)
#从多项分布中抽取新的主题
newtopic = Sample.MultSample(p)
self.nw[wid][newtopic] += 1
self.nd[d][newtopic] += 1
self.nwsum[newtopic] += 1
self.ndsum[d] += 1
return newtopic
def ComputeTransProb(self, d, w):
"""
对第d篇文档的第w个词
计算Gibbs Sampling过程中的传递概率
:param d: 文档编号
:param w: 词在文档中的编号
"""
#用于平滑
Wbeta = self.W * self.beta
Kalpha = self.K * self.alpha
#第d篇文档，第w个词对应的id
wid = self.IDListSet[d][w]
p = ListUtil.Initial(self.K, 0.0)
for k in xrange(self.K):
#p[k] = p(w|k)*p(k|d)   k为主题
p[k] = (float(self.nw[wid][k]) + self.beta) / (float(self.nwsum[k]) + Wbeta) * (float(self.nd[d][k]) + self.alpha) / (float(self.ndsum[d]) + Kalpha)
return p

其实上面的计算法则就是p(z(i)=k|z', w, alpha,beta) = p(z'(i)=k|d)*p(w|z'(i)=k)，z'就是代表将当前词w剔除之后的主题分布，z(i)对应当前词w的主题，这里就跟文章开头的生成模型的原理呼应上了，我们先以一定的概率选择主题（p(topic|doc)），然后在从主题包含的词中抽取相应的词（p(word|topic)），吉布斯抽样也是沿着doc->topic->word这样的方向进行的，给一张图大家可能更好理解，theta[m][k]表示第m篇文档，第k个主题出现的概率；phi[k][t]表示主题k中，词t出现的概率

计算传递概率之后，我们再从多项分布中抽取相应的主题（这个函数对应MultSample函数，下面为整个Sample的代码）

def UniSample(K):"""产生从O到K－1的整数:param K: 主题个数"""return RandomNumber.RandInt(0, K - 1)def MultSample(ProbList):"""从多项分布ProbList中采样, ProbList表示剔除当前词之后的主题分布:param ProbList: 多项分布"""size = len(ProbList)for i in xrange(1, size):ProbList[i] += ProbList[i - 1]#随机产生一个［0，1）的小数u = RandomNumber.RandFloat()res = 0for k in xrange(size):if ProbList[k] >= u * ProbList[size - 1]:#抽样结果res = kbreak#res为抽样后的主题编号return res

其实吉布斯抽样的目的就是为乐得到在 图模型中的theta和phi，那要怎么样计算theta和phi呢？其实很简单的，吉布斯抽样中，doc-topic, topic-word矩阵的计数是变化的，当抽样收敛之后，我们就得到了最后的计数，通过这些计数来计算频率就好了

    def ComputTheta(self):"""计算p(z|d)矩阵size:D * Kp(z|d) = theta[d][z]"""for d in xrange(self.D):for k in xrange(self.K):self.theta[d][k] = (float(self.nd[d][k]) + self.alpha) / (float(self.ndsum[d]) + self.K * self.alpha)def ComputePhi(self):"""计算p(w|z)size:K * Wp(w|z) = phi[z][w]"""for k in xrange(self.K):for w in xrange(self.W):self.phi[k][w] = (self.nw[w][k] + self.beta) / (self.nwsum[k] + self.W * self.beta)

为了防止0概率的出现，我们分别用alpha和beta做了平滑

首先来看phi[k][w]，我们用频率计数nw[w][k]（主题w中包含的词k的数目）来计算，不幸的是，这个值有可能为0，在一连串的乘式中，0的出现会使其他项毫无意义，所以我们要避免这种情况，怎么避免呢，我们可以假设主题k中事先已经存在了词典中的所有词，然后，再用我们得到的收敛之后的 主题－词计数去更新里面的内容，这样，就可以保证不会出现0概率（因为所有词至少出现一次，同时，主题k中的词数nwsum[k]增加了W），那么phi[k][w]=(nw[w][k] + 1) / (nwsum[k] + W)，这就是鼎鼎大名的 “拉普拉斯平滑”，不过工程上，1对于概率的影响还是太大了（比如1/2和2/3），所以，我们吧1换成了更小的beta（一般是0.01），这样对于概率的影响就变得很小了，也就得到了上面的公式，theta[d][k]也同理

整个参数训练的过程如下

    def estimate(self):"""LDA参数估计"""for i in xrange(1, self.NumberOfIterations + 1):for d in xrange(self.D):for w in xrange(len(self.IDListSet[d])):newtopic = self.sampling(d, w)#为当前词分派新主题self.Z[d][w] = newtopicif i % self.SaveStep == 0:#计算当前的迭代结果self.ComputTheta()self.ComputePhi()self.SaveTempRes(i)

LDA模型的参数训练部分就讲到这里，下一篇跟大家分享一下LDA的参数推导

希望看到的亲能帮我指出文中的错误，文笔不好，大家多见谅