function [alpha,beta,eta] = ldabayes(d,l,k,emmax,demmax)
% [alpha,beta,eta] = ldabayes(d,l,k,[emmax,demmax])
% Latent Dirichlet Allocation, full Bayesian model.
% Copyright (c) 2004 Daichi Mochihashi, all rights reserved.
% $Id: ldabayes.m,v 1.2 2004/10/26 02:13:26 dmochiha Exp $
% d      : data of documents
% l      : # of lexicon
% 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 < 5
  demmax = 20;
  if nargin < 4
    emmax = 100;
  end
end
n = length(d);
eta = zeros(l,1);
beta = ones(l,k) / l;
alpha = normalize(rand(1,k));
gammas = zeros(n,k);
lik = 0;
plik = lik;
tic;

for j = 1:emmax
  fprintf(1,'iteration %d/%d..\t',j,emmax);
  % vb-estep
  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
  alpha = newton_alpha(gammas);
  eta = newton_eta(repmat(eta,1,k) + betas);
  beta = mnormalize(repmat(eta,1,k) + betas,1);
  % converge?
  lik = lda_lik(d,beta,gammas);
  fprintf(1,'likelihood = %g\t',lik);
  if (j > 1) && converged(lik,plik,1.0e-4)
    fprintf(1,'\nconverged.\n');
    return;
  end
  plik = lik;
  % ETA
  elapsed = toc;
  fprintf(1,'ETA:%s (%d sec/step)\r', ...
	  rtime(elapsed * (emmax / j  - 1)),round(elapsed / j));
end
fprintf(1,'\n');

% $Id: ldabayes.m,v 1.2 2004/10/26 02:13:26 dmochiha Exp $
