PLSA的Python实现

概率潜在语义分析是在LSA基础上提出的分析“文档-主题-词”之间的关系。

转自概率语言模型及其变形系列(1)-PLSA及EM算法 - Coding for Dreams - 博客频道 - CSDN.NET http://blog.csdn.net/yangliuy/article/details/8330640

1.首先定义文档集合D和词集合W及共现频率矩阵N，Z代表隐含的主题。每个主题在所有词项上服从Multinomial （多项式）分布，每个文档在所有主题上服从Multinomial 分布。该方法使用概率模型模拟潜在的语义空间，将文档和词映射到同一语义空间，把一个词汇在不同的潜在语义变量上的概率分布理解为词汇的不同意思，从而成功地解决“一词多义”和“一义多词”的问题。文档和词同现的联合概率为：

PLSA的概率生成图模型如下：

(1) 以的概率选中文档；

(2) 以的概率选中主题；

(3) 以的概率产生一个单词。

 

而和分布对应了两组Multinomial 分布，需要估计这两组分布的参数。下面给出用EM算法估计PLSA参数的详细推导过程。

EM算法的步骤是：

(1)E步骤：求隐含变量Given当前估计的参数条件下的后验概率。

(2)M步骤：最大化Complete data对数似然函数的期望，此时我们使用E步骤里计算的隐含变量的后验概率，得到新的参数值。

两步迭代进行直到收敛。

当原始数据的似然函数很复杂时，我们通过增加一些隐含变量来增强我们的数据，得到“complete data”,而“complete data”的似然函数更加简单，方便求极大值。于是，原始的数据就成了“incomplete data”。我们将会看到，我们可以通过最大化“complete data”似然函数的期望来最大化"incomplete data"的似然函数，以便得到求似然函数最大值更为简单的计算途径。

针对我们PLSA参数估计问题，在E步骤中，直接使用贝叶斯公式计算隐含变量在当前参数取值条件下的后验概率，有

在这个步骤中，我们假定所有的和都是已知的，因为初始时随机赋值，后面迭代的过程中取前一轮M步骤中得到的参数值。

在M步骤中，我们最大化Complete data对数似然函数的期望。在PLSA中，Incomplete data 是观察到的，隐含变量是主题，那么complete data就是三元组，其期望是

注意这里是已知的，取的是前面E步骤里面的估计值。下面我们来最大化期望，这又是一个多元函数求极值的问题，可以用拉格朗日乘数法。拉格朗日乘数法可以把条件极值问题转化为无条件极值问题，在PLSA中目标函数就是,约束条件是

由此我们可以写出拉格朗日函数

这是一个关于和的函数，分别对其求偏导数，我们可以得到

注意这里进行过方程两边同时乘以和的变形，联立上面4组方程，我们就可以解出M步骤中通过最大化期望估计出的新的参数值

解方程组的关键在于先求出,其实只需要做一个加和运算就可以把的系数都化成1，后面就好计算了。

然后使用更新后的参数值，我们又进入E步骤，计算隐含变量 Given当前估计的参数条件下的后验概率。如此不断迭代，直到满足终止条件。

注意到我们在M步骤中还是使用对Complete Data的MLE，那么如果我们想加入一些先验知识进入我们的模型，我们可以在M步骤中使用MAP估计。正如文本语言模型的参数估计-最大似然估计、MAP及贝叶斯估计中投硬币的二项分布实验中我们加入“硬币一般是两面均匀的”这个先验一样。而由此计算出的参数的估计值会在分子分母中多出关于先验参数的preduo counts,其他步骤都是一样的。具体可以参考Mei Qiaozhu 的Notes。

PLSA的实现网上有很多实现code。

# -*- coding: utf-8 -*-

import math
import operator
import random
import gzip
import sys
import marshal

def cos_sim(p, q):
    sum0 = sum(map(lambda x:x*x, p))
    sum1 = sum(map(lambda x:x*x, q))
    sum2 = sum(map(lambda x:x[0]*x[1], zip(p, q)))
    return sum2/(sum0**0.5)/(sum1**0.5)

def _rand_mat(sizex, sizey):
    ret = []
    for i in xrange(sizex):
        ret.append([])
        for _ in xrange(sizey):
            ret[-1].append(random.random())
        norm = sum(ret[-1])
        for j in xrange(sizey):
            ret[-1][j] /= norm
    return ret

class Plsa:

    def __init__(self, corpus, topics=2):
        self.topics = topics
        self.corpus = corpus
        self.docs = len(corpus)
        self.each = map(sum, map(lambda x:x.values(), corpus))
        self.words = max(reduce(operator.add, map(lambda x:x.keys(), corpus)))+1
        self.likelihood = 0
        self.zw = _rand_mat(self.topics, self.words)
        self.dz = _rand_mat(self.docs, self.topics)
        self.dw_z = None
        self.p_dw = []
        self.beta = 0.8

    def save(self, fname, iszip=True):
        d = {}
        for k, v in self.__dict__.items():
            if hasattr(v, '__dict__'):
                d[k] = v.__dict__
            else:
                d[k] = v
        if sys.version_info[0] == 3:
            fname = fname + '.3'
        if not iszip:
            marshal.dump(d, open(fname, 'wb'))
        else:
            f = gzip.open(fname, 'wb')
            f.write(marshal.dumps(d))
            f.close()

    def load(self, fname, iszip=True):
        if sys.version_info[0] == 3:
            fname = fname + '.3'
        if not iszip:
            d = marshal.load(open(fname, 'rb'))
        else:
            try:
                f = gzip.open(fname, 'rb')
                d = marshal.loads(f.read())
            except IOError:
                f = open(fname, 'rb')
                d = marshal.loads(f.read())
            f.close()
        for k, v in d.items():
            if hasattr(self.__dict__[k], '__dict__'):
                self.__dict__[k].__dict__ = v
            else:
                self.__dict__[k] = v

    def _cal_p_dw(self):
        self.p_dw = []
        for d in xrange(self.docs):
            self.p_dw.append({})
            for w in self.corpus[d]:
                tmp = 0
                for _ in range(self.corpus[d][w]):
                    for z in xrange(self.topics):
                        tmp += (self.zw[z][w]*self.dz[d][z])**self.beta
                self.p_dw[-1][w] = tmp

    def _e_step(self):
        self._cal_p_dw()
        self.dw_z = []
        for d in xrange(self.docs):
            self.dw_z.append({})
            for w in self.corpus[d]:
                self.dw_z[-1][w] = []
                for z in xrange(self.topics):
                    self.dw_z[-1][w].append(((self.zw[z][w]*self.dz[d][z])**self.beta)/self.p_dw[d][w])

    def _m_step(self):
        for z in xrange(self.topics):
            self.zw[z] = [0]*self.words
            for d in xrange(self.docs):
                for w in self.corpus[d]:
                    self.zw[z][w] += self.corpus[d][w]*self.dw_z[d][w][z]
            norm = sum(self.zw[z])
            for w in xrange(self.words):
                self.zw[z][w] /= norm
        for d in xrange(self.docs):
            self.dz[d] = [0]*self.topics
            for z in xrange(self.topics):
                for w in self.corpus[d]:
                    self.dz[d][z] += self.corpus[d][w]*self.dw_z[d][w][z]
            for z in xrange(self.topics):
                self.dz[d][z] /= self.each[d]

    def _cal_likelihood(self):
        self.likelihood = 0
        for d in xrange(self.docs):
            for w in self.corpus[d]:
                self.likelihood += self.corpus[d][w]*math.log(self.p_dw[d][w])

    def train(self, max_iter=100):
        cur = 0
        for i in xrange(max_iter):
            print '%d iter' % i
            self._e_step()
            self._m_step()
            self._cal_likelihood()
            print 'likelihood %f ' % self.likelihood
            if cur != 0 and abs((self.likelihood-cur)/cur) < 1e-8:
                break
            cur = self.likelihood

    def inference(self, doc, max_iter=100):
        doc = dict(filter(lambda x:x[0]<self.words, doc.items()))
        words = sum(doc.values())
        ret = []
        for i in xrange(self.topics):
            ret.append(random.random())
        norm = sum(ret)
        for i in xrange(self.topics):
            ret[i] /= norm
        tmp = 0
        for _ in xrange(max_iter):
            p_dw = {}
            for w in doc:
                p_dw[w] = 0
                for _ in range(doc[w]):
                    for z in xrange(self.topics):
                        p_dw[w] += (ret[z]*self.zw[z][w])**self.beta
            # e setp
            dw_z = {}
            for w in doc:
                dw_z[w] = []
                for z in xrange(self.topics):
                    dw_z[w].append(((self.zw[z][w]*ret[z])**self.beta)/p_dw[w])
            # m step
            ret = [0]*self.topics
            for z in xrange(self.topics):
                for w in doc:
                    ret[z] += doc[w]*dw_z[w][z]
            for z in xrange(self.topics):
                ret[z] /= words
            # cal likelihood
            likelihood = 0
            for w in doc:
                likelihood += doc[w]*math.log(p_dw[w])
            if tmp != 0 and abs((likelihood-tmp)/tmp) < 1e-8:
                break
            tmp = likelihood
        return ret

    def post_prob_sim(self, docd, q):
        sim = 0
        for w in docd:
            tmp = 0
            for z in xrange(self.topics):
                tmp += self.zw[z][w]*q[z]
            sim += docd[w]*math.log(tmp)
        return sim

######### unittest #################################
import unittest

class TestPlsa(unittest.TestCase):

    def test_train(self):
        corpus = [{0:2,3:5},{0:5,2:1},{1:2,4:5}]
        p = Plsa(corpus)
        p.train()
        self.assertTrue(cos_sim(p.dz[0], p.dz[1])>cos_sim(p.dz[0], p.dz[2]))
        self.assertTrue(p.post_prob_sim(p.corpus[0], p.dz[1])>p.post_prob_sim(p.corpus[0], p.dz[2]))

    def test_inference(self):
        corpus = [{0:2,3:5},{0:5,2:1},{1:2,4:5}]
        p = Plsa(corpus)
        p.train()
        z = p.inference({0:4, 6:7})
        self.assertTrue(abs(cos_sim(p.dz[0], p.dz[1])-cos_sim(p.dz[0], z))<1e-8)

if __name__ == '__main__':
    unittest.main()