function eta = newton_eta(lambdas,maxiter,ini_eta)
% eta = newton_eta(lambdas,[maxiter])
% Newton-Raphson iteration of LDA lexical Dirichlet prior.
% lambdas : matrix of Dirichlet posteriors (V * K)
% maxiter : # of maximum iteration of Newton-Raphson (default 50)
% $Id: newton_eta.m,v 1.3 2004/10/26 02:13:02 dmochiha Exp $
[V,K] = size(lambdas);
if ~(K > 1)
  eta = lambdas(:,1);
  return;
end
if nargin < 3
  ini_eta = mean(lambdas,2) * 1.0e-6;
  if nargin < 2
    maxiter = 50;
  end
end

p = sum(psi(lambdas),2) - sum(psi(sum(lambdas,1)));
l = 0;
eta = ini_eta;
peta = zeros(V,1);

for t = 1:maxiter
  l = l + 1;
  eta0 = sum(eta);
  peta0 = psi(eta0);
  c = sum((K * (psi(eta) - psi(eta0)) - p) ./ psi(1,eta)') / ...
      (1 / psi(1,eta0) - sum(1 ./ psi(1,eta)));
  
  eta = eta - (c - p + K * (psi(eta) - psi(eta0))) ./ (K * psi(1,eta)');
  if any(eta < 0)
    eta = newton_eta(lambdas,maxiter,ini_eta / 10); % try again!
    return;
  end
  
  if (l > 1) && converged(eta,peta,1.0e-4)
    break;
  end
  peta = eta;
end
