function [hmm,loglikes] = esthmm(X, hmm, epsilon, maxittr, minittr, dispfreq)
% Starting with an initial hmm, update its parameters so as to maximize
% the likelihood of the observed sequences given in X. X can either be a
% cell array of row vectors X{i}, each representing an output sequence (of
% possibly different length), or X can be a matrix with each row of X,
% X(i,:), representing an output sequence (in which case all output
% sequences must be of the same length). As a special case,  X might be a 
% single row vector representing one output sequence.

  if nargin<3, epsilon = 1e-4; end;
  if nargin<4, maxittr = 100; end;
  if nargin<5, minittr = 2; end;
  if nargin<6, dispfreq = 1; end;
  if isnumeric(X), X=num2cell(X,2); end;

  loglikes = []; 
  cont = 1; nittr=0; lastll = -inf;
  while (cont),
    [hmm,ll] = hmmEMstep(X,hmm);
    loglikes = [loglikes;ll];
    cont = (((nittr<minittr) | (ll>lastll+epsilon)) ...
	    & nittr<maxittr) ;
    lastll = ll;
    nittr = nittr+1;
    if (mod(nittr,dispfreq)==0), disp(ll); end;
  end;


  % ------------------------------

function [hmm,ll] = hmmEMstep(X,hmm)
  % one EM iteration
    
  [ns,no] = size(hmm.Po);
  ll = 0;
 
  Ns0 = zeros(1,ns);  % soft counts for the starting state: 
  % Ns0(i) = expected number of times initial state was i
  Nss = zeros(ns,ns); % soft counts for state transitions
  % Nss(i,j) = expected number of times state i was visited, followed by
  % state j
  Nsx = zeros(ns,no); % soft counts for state/observation pairs
  % Nsx(i,j) = expected number of times state i was visited emitted
  % output j

  % E-step
  for r=1:length(X),			% iterate over all output sequences
    x = X{r};
    [Ps, Pss, logPx] = hmmposteriors(x,hmm);
    % Strictly speaking, the E step is just hmmposteriors. But we don't want
    % to keep around the complete posterior probabilities for each sequence
    % (which are of size ns*ns*length(x)), until we finish the "E"xpectation
    % on all sequences. Instead, we right-away get from the posteriors
    % what we really need-- the softcounts  of starting states, state
    % transitions and state/observation pairs.
    
    % Add the softcounts from this output sequence to the totals
    Ns0 = Ns0 + Ps(:,1)';		% Nso(i) = expected  number of 
                                        % times the initial state was i
    Nsx = Nsx + Ps*op(x','==',1:no);	% Nsx(i,j) = expected number
                                        % of times we were in state i AND
                                        % emitted output j
    Nss = Nss + sum(Pss,3);		% Nss(i,j) = expected number
                                        % of times we were in state i and
                                        % then transitioned to j
    ll = ll + logPx;
  end;
  % complete the "M"aximization step-- update the parameters to the maximum
  % likelihood setting according to the expected state occupations
  hmm.Pt = op(Nss,'./',sum(Nss,2)) ;
  hmm.Po = op(Nsx,'./',sum(Nsx,2));
  hmm.p0 = Ns0/sum(Ns0);

