LDA源碼分析（matlab版）

LDA全稱是Latent Dirichlet Allocation。關於LDA的理論知識，可以參見如下參考維基百科。這裡具體講解一下LDA的源碼分析（matlab）
代碼原作者：Daichi Mochihashi
源碼：http://download.csdn.net/detail/nuptboyzhb/5305145
一.LDA源碼在matlab環境下的執行
1.環境配置
將matlab的工作目錄切換到代碼所在目錄
2.調用主函數
>> [alpha,beta] =ldamain(‘train’,20);%訓練資料檔案train
分類數20
二.訓練資料train的資料格式
如：
<feature_id>:<count>  特徵的標號：對應的個數
對於文檔而言，特徵id表徵的是某個單詞，數目則表示單詞出現的次數
train中的每一行表示一個文檔，如下：
1:1 2:4 5:2
1:2 3:3 5:1 6:1 7:1
2:4 5:1 7:1
注意：LDA中train的資料格式與SVM中的差異，在SVM中，訓練資料的格式如下：
與svm的訓練資料格式，相似但不同
<label> <index1>:<value1> <index2>:<value2> ... 
也就是說，SVM中的訓練資料，都有“標籤”。SVM是一個監督學習的過程。而LDA是非監督學習。
三.代碼中的變數意義

n 整數，表示文檔數L  整數，表示單詞數beta 二維數組，行代表單詞，列代表主題，矩陣單元代表某主題產生某詞的機率alpha 數組，對應dirichlet分布的參數k 整數，代表主題數，這個是由使用者佈建的值gamma 一維數組，變分推理中後驗dirichlet分布的參數gammas 充分統計量一維數組，形式同gamma，用於在m-step估計alpha的值q 二維數組，行代表文檔裡的單詞，列代表主題，矩陣單元代表文檔中某主題產生某詞的機率betas 充分統計量,二維數組，形式同q，該變數用於在e-step統計資訊，供m-step估計beta使用。

四.LDA的源碼總覽

雖然LDA的源碼有很多.m檔案，但是主要的檔案時lda.m vbem.m and newton_alpha.m三個檔案。
五.核心程式碼分析

lda.m

function [alpha,beta] = lda(d,k,emmax,demmax)% Latent Dirichlet Allocation, standard model.% d      : data of documents% k      : # of classes to assume% emmax  : # of maximum VB-EM iteration (default 100)% demmax : # of maximum VB-EM iteration for a document (default 20)if nargin < 4  demmax = 20;  if nargin < 3    emmax = 100;  endendn = length(d);l = features(d);%初始化beta = mnormalize(rand(l,k),1);alpha = normalize(fliplr(sort(rand(1,k))));%對應dirichlet分布的參數 gammas = zeros(n,k);ppl = 0;pppl = ppl;tic;fprintf(1,'number of documents      = %d\n', n);fprintf(1,'number of words          = %d\n', l);fprintf(1,'number of latent classes = %d\n', k);for j = 1:emmax  fprintf(1,'iteration %d/%d..\t',j,emmax);  % vb-estep 輸入alpha和beta計算gammas  betas = zeros(l,k);  for i = 1:n % 對每個文檔進行計算    [gamma,q] = vbem(d{i},beta,alpha,demmax);    gammas(i,:) = gamma; %儲存每個文檔的值    betas = accum_beta(betas,q,d{i});  end  % vb-mstep 最大化似然函數(gammas)，求解alpha和beta  alpha = newton_alpha(gammas);  beta = mnormalize(betas,1);  % converge?  ppl = lda_ppl(d,beta,gammas);  fprintf(1,'PPL = %g\t',ppl);  if (j > 1) && converged(ppl,pppl,1.0e-4)    if (j < 5)      fprintf(1,'\n');                 % 迭代次數過少try again!      [alpha,beta] = lda(d,k,emmax,demmax); return;    end    fprintf(1,'\nconverged.\n');    return;  end  pppl = ppl;  % ETA  elapsed = toc;  fprintf(1,'ETA:%s (%d sec/step)\r', ...      rtime(elapsed * (emmax / j  - 1)),round(elapsed / j));endfprintf(1,'\n');

未經允許，文章不得用於商業目的！！