PLSA模型的再理解以及原始碼分析

   

    之前寫過一篇PLSA的博文，其中的收穫就是知道PLSA是LSA在概率層面的擴充套件，知道了PLSA是一種主題模型，知道了PLSA中的引數估計使用的是EM演算法。當時我就認為，這樣子經典好用的演算法，我是會回頭再來理解它的，這樣子才會有更加深刻的心得。所以有了這篇PLSA模型的再理解。

1. 兩種思路解PLSA模型

    參考了很多資料，發現大體上有兩種解決PLSA模型的思路。下面我們大致說一下它們的思路。

思路一：EM演算法中需要更新兩個概率

    PLSA模型的示意圖如下：

其中包括的概率有：

也就是說我們知道的是P(d,w)，需要求解的是P(z,d)和P(w,z)。

    我們根據P(d,w)構建似然函式，而且我們需要最大化這個似然函式，如下：

需要明確的是，這個似然函式中的自變數是P(z,d)和P(w,z),不巧的是，這兩個變數在Log函式的內部是相加的關係，我們知道這為對似然函式的求導帶來了很多麻煩。所以我們使用的是EM演算法。

   EM演算法的基本實現思想是：

(1)E步驟：求隱含變數Given當前估計的引數條件下的後驗概率。

(2)M步驟：最大化Complete data對數似然函式的期望，此時我們使用E步驟裡計算的隱含變數的後驗概率，得到新的引數值。

兩步迭代進行直到收斂。

    具體而言，

    在PLSA模型中的E步驟，我們直接用貝葉斯公式計算隱含變數zk在當前引數取值條件下的後驗概率，記住是後驗概率，如下：

這裡的都是已知的，因為我們一開始會對它們設定一個初始值，這樣後面的迭代就可以依據於此。使用上面的公式我們計算出了隱含變數的後驗概率，然後我們進入M步驟。

    M步驟中，我們需要計算complete data的對數似然函式的期望，而且要最大化這個期望。那麼首先我們要先計算出這個對數似然期望，如下：

（整篇演算法就是這個無法理解，為什麼期望是這樣子的公式）

注意，這個時候的P(z|d,w)是已知的了，因為在E步驟中已經計算過了，這時候我們需要最大化這個期望，其中包含的自變數有P(z,d)和P(w,z)，所以我們需要分別對它們求偏導，令偏導為0：

解方程可得如下的結果：

這個就是我們在M步驟中最大化complete data對似然函式的期望，得到的兩個引數的更新公式。

    不斷重複EM兩個步驟，可以得到最終的兩個引數P(z,d)和P(w,z)的值。

思路二：EM演算法中需要更新3個引數

    這個思路和上面思路的區別就在於對P(d,w)的展開公式使用的不同，思路二使用的是3個概率來展開的，如下：

這樣子我們後面的EM演算法的大致思路都是相同的，就是表達形式變化了，最後得到的EM步驟的更新公式也變化了。當然，思路二得到的是3個引數的更新公式。如下：

你會發現，其實多了一個引數是P(z)，引數P（d|z）變化了（之前是P(z|d)），然後P(w|z)是不變的，計算公式也相同。

2. PLSA演算法的應用效果

    對於資料X，如下：

我們使用思路二的PLSA演算法，更新許多次之後的概率為：

iteration 27:
p(d/z)=
0.0000 0.0769 0.2931 0.2934 0.0000 0.0733 0.1467 0.1166 0.0000
0.1952 0.3223 0.0002 0.0000 0.1952 0.0000 0.0000 0.0918 0.1952

p(w/z)=
0.0733 0.0733 0.0000 0.0823 0.2878 0.0000 0.0000 0.1467 0.0000 0.2200 0.1166 0.0000
0.0651 0.0652 0.1302 0.1222 0.0049 0.1302 0.1302 0.0000 0.1302 0.0000 0.0918 0.1302

p(z)=
0.4702 0.5298

迭代了27次，得到了這三個概率。那麼這三個概率有什麼用呢？

   我們知道LSA模型是利用SVD把原來的資料矩陣X降維，使用兩個詞之間的語義關係被挖掘了出來。兩個原始空間中不相關的卻語義上相關的詞，在LSA模型中是比較相關的，這個是LSA直接的作用。那麼PLSA有什麼作用呢？

   給定一個文件d，我們可以將其分類到一些主題詞類別下。

    PLSA演算法可以通過訓練樣本的學習得到三個概率，而對於一個測試樣本，其中P(w|z)概率是不變的，但是P(z)和P(d|z)都是需要重新更新的，我們也可以使用上面的EM演算法，假如測試樣本d的資料，我們可以得到新學習的P(z)和P(d|z)引數。這樣我們就可以計算：

為什麼要計算P(z|d)呢？因為給定了一個測試樣本d，要判斷它是屬於那些主題的，我們就需要計算P(z|d)，就是給定d，其在主題z下成立的概率是多少，不就是要計算嗎。這樣我們就可以計算文件d在所有主題下的概率了。

    這樣既可以把一個測試文件劃歸到一個主題下，也可以給訓練文件打上主題的標記，因為我們也是可以計算訓練文件它們的的。如果從這個應用思路來說，思路一說似乎更加直接，因為其直接計算出來了。

    如果從上面這個應用的角度來說，LSA和PLSA的目的就是從一群文件集中找到潛在的語義因子latent factors。由於提取到的主題詞比文件中的詞的數量要少很多，而且我們在學習的過程中不需要知道文件的型別資訊，所以說LSA和PLSA是無監督的特徵降維方法。

    

3. PLSA原始碼分析

    知道PLSA演算法的理論和作用，我們還要會用，下面提供一個MATLAB原始碼，人家寫的挺簡單的，我加了一些分析，幫助理解。

main函式：

% function demo
% 
close all; clear all;
load data.mat;
k=2;
[pz pdz pwz pzdw]=plsa(X,k);  % 返回的是4種概率
imagesc(pz); colormap(gray); colorbar; 
title('Probability of topics pz'); 
xlabel('topic z_i'); ylabel('probability');

figure;
imagesc(pdz); colormap(gray); colorbar;
title('Category probability of document pdz'); 
xlabel('document d_i (pdz_{(i,j)} is the prob. that d_i belongs to z_j)'); ylabel('topic z_j');

figure;
imagesc(pwz); colormap(gray); colorbar;
title('Category probability of word pwz'); 
xlabel('word w_i (pwz_{(i,j)} is the prob. that w_i belongs to z_j)'); ylabel('topic z_j');

PLSA子函式：

function [pz pdz pwz pzdw]=plsa(X,k)
% function [pz pdz pwz pzdw]=plsa(X,k)
% Return pLSA probability matrix p of m*n matrix X
% X is the document-word co-occurence matrix
% k is the number of the topics--z
% document--collums,word--rows

err=0.0001;
x = X; 
[m n]=size(x);  % m個詞，n個文件
pz=rand(1,k);  % 隨機初始化pz, k個主題詞
pz2=pz;
pz2(1)=pz2(1)+2*err;

% 初始化的都是隨機值
pdz=rand(k,n);  % 二維
pwz=rand(k,m);  % 二維
pzdw=rand(m,n,k);     % 三維   %initialize
h=0;

deno=zeros(1,k);        %denominator of p(d/z) and p(w/z)
denopzdw=zeros(m,n);       %denominator of p(z/d,w)
numepdz=zeros(k,n);        %numerator of p(d/z)
numepwz=zeros(k,m);        % numerator of p(w/z)
R=sum(sum(x));  %  所有的文件的詞數總和

for ki=1:k
    for i=1:m
        for j=1:n
            % P(w|z)和P(d|z)的分母，是每個zk一個值，所以是個向量
            deno(ki)=deno(ki)+x(i,j)*pzdw(i,j,ki); % 疊加，計算的是P(w|z)和P(d|z)的分母
        end
    end
end
% p(d/z)
for ki=1:k
    for j=1:n
        for i=1:m
            % P(d|z)的分子，每個zk和dj一個值，所以是一個矩陣
            numepdz(ki,j)=numepdz(ki,j)+x(i,j)*pzdw(i,j,ki); % 疊加，計算的是P(d|z)的分子
        end
        pdz(ki,j)=numepdz(ki,j)/deno(ki);  % 分子除以分母，一個矩陣
    end
end
% disp(pdz);

% p(w/z)
for ki=1:k
    for i=1:m
        for j=1:n
            % 計算的是P(w|z)的分子，每個zk和wi一個值，一個矩陣
            numepwz(ki,i)=numepwz(ki,i)+x(i,j)*pzdw(i,j,ki);
        end
        % 分子除以分母，矩陣中的每個元素
        pwz(ki,i)=numepwz(ki,i)/deno(ki);
    end
end
% disp(pwz);

% p(z)
for ki=1:k
    pz(ki)=deno(ki)./R;  % 即為前面的分母
end
%denominator of p(z/d,w)
for i=1:m
    for j=1:n
        for ki=1:k
            % 計算的是P(z|d,w)的分母，對zk疊加，每個wi和dj一個值，一個矩陣
            denopzdw(i,j)=denopzdw(i,j)+pz(ki)*pdz(ki,j)*pwz(ki,i);
        end
    end
end
% p(z/d,w)
for i=1:m
     for j=1:n
         for ki=1:k
             % 計算的是P(z|d,w)，每個zk和wi和dj一個值，一個三維矩陣
             pzdw(i,j,ki)=pz(ki)*pdz(ki,j)*pwz(ki,i)/denopzdw(i,j);
         end
     end
end
% fprintf('p(z/d,w)=\n');
% disp(pzdw)

%  iteration 
iteration=0;
fprintf('iteration:\n');
while abs(pz2(1)-pz(1))>err || abs(pz2(2)-pz(2))>err
    iteration=iteration+1;  % 迭代次數
    deno=zeros(1,k);        %denominator of p(d/z) and p(w/z)
    denopzdw=zeros(m,n);       %denominator of p(z/d,w)
    numepdz=zeros(k,n);        %numerator of p(d/z)
    numepwz=zeros(k,m);        % numerator of p(w/z)
    fprintf('iteration %d:\n',iteration);
    for ki=1:k
        for i=1:m
            for j=1:n
                deno(ki)=deno(ki)+x(i,j)*pzdw(i,j,ki); % 繼續迭代，deno回到0，因為需要疊加，不過pzdw就是上次的新值
            end
        end
    end
% p(d/z)
    for ki=1:k
        for j=1:n
            for i=1:m
                numepdz(ki,j)=numepdz(ki,j)+x(i,j)*pzdw(i,j,ki);  % 計算公式同上，只不過此時的pzdw是更新過後的值
            end
            pdz(ki,j)=numepdz(ki,j)/deno(ki);
        end
    end
    fprintf('p(d/z)=\n');
    disp(pdz)
    % p(w/z)
    for ki=1:k
        for i=1:m
            for j=1:n
                numepwz(ki,i)=numepwz(ki,i)+x(i,j)*pzdw(i,j,ki);
            end
            pwz(ki,i)=numepwz(ki,i)/deno(ki);
        end
    end
    fprintf('p(w/z)=\n');
    disp(pwz)

    % p(z)
    pz=pz2;
    for ki=1:k
        pz2(ki)=deno(ki)./R;
    end
    fprintf('p(z)=\n');
    disp(pz2)
    %denominator of p(z/d,w)
    for i=1:m
        for j=1:n
            for ki=1:k
                denopzdw(i,j)=denopzdw(i,j)+pz2(ki)*pdz(ki,j)*pwz(ki,i);
            end
        end
    end
    % p(z/d,w)
    for i=1:m
         for j=1:n
             for ki=1:k
                 pzdw(i,j,ki)=pz2(ki)*pdz(ki,j)*pwz(ki,i)/(denopzdw(i,j)+eps);
             end
         end
    end
end  %end while

return;

參考：

http://www.cnblogs.com/ywl925/p/3552815.html

http://blog.csdn.net/huangxy10/article/details/8091478

http://blog.csdn.net/huangxy10/article/details/8091558

LSA和PLSA 唐克坦

http://blog.csdn.net/puqutogether/article/details/41720073