added comment

chapter13/HMM/hmmEm.m

@@ -28,9 +28,9 @@
     energy(iter) = sum(log(c(c>0)));
     if energy(iter)-energy(iter-1) < tol*abs(energy(iter-1)); break; end   % check likelihood for convergence
 %     M-step 
-    A = normalize(A.*(alpha(:,1:n-1)*bsxfun(@times,beta(:,2:n).*M(:,2:n),1./c(2:end))'),2);
-    s = gamma(:,1);
-    M = bsxfun(@times,gamma*X',1./sum(gamma,2))*X;
+    A = normalize(A.*(alpha(:,1:n-1)*bsxfun(@times,beta(:,2:n).*M(:,2:n),1./c(2:end))'),2);      % 13.19
+    s = gamma(:,1);                                                                             % 13.18
+    M = bsxfun(@times,gamma*X',1./sum(gamma,2))*X;                                          
 end
 energy = energy(2:iter);
 model.A = A;

chapter13/HMM/hmmEm.m~

@@ -0,0 +1,40 @@
+function [model, energy] = hmmEm(x, init)
+% EM algorithm to fit the parameters of HMM model (a.k.a Baum-Welch algorithm)
+%   x: 1 x n sequence of observations
+%   init: model or k
+% Written by Mo Chen ([email protected]).
+n = size(x,2);
+d = max(x);
+X = sparse(x,1:n,1,d,n);
+
+if isstruct(init)   % init with a model
+    A = init.A;
+    E = init.E;
+    s = init.s;
+elseif numel(init) == 1  % random init with latent k
+    k = init;
+    A = normalize(rand(k,k),2);
+    E = normalize(rand(k,d),2);
+    s = normalize(rand(k,1),1);
+end
+M = E*X;
+
+tol = 1e-4;
+maxIter = 100;
+energy = -inf(1,maxIter);
+for iter = 2:maxIter
+%     E-step
+    [gamma,alpha,beta,c] = hmmSmoother_(M,A,s);
+    energy(iter) = sum(log(c(c>0)));
+    if energy(iter)-energy(iter-1) < tol*abs(energy(iter-1)); break; end   % check likelihood for convergence
+%     M-step 
+    A = normalize(A.*(alpha(:,1:n-1)*bsxfun(@times,beta(:,2:n).*M(:,2:n),1./c(2:end))'),2);      
+    s = gamma(:,1);
+    M = bsxfun(@times,gamma*X',1./sum(gamma,2))*X;
+end
+energy = energy(2:iter);
+model.A = A;
+model.E = E;
+model.s = s;
+
+

chapter13/HMM/hmmFilter_.m

@@ -9,6 +9,6 @@
 alpha = zeros(K,T);
 [alpha(:,1),energy(1)] = normalize(s.*M(:,1),1);
 for t = 2:T
-    [alpha(:,t),energy(t)] = normalize((At*alpha(:,t-1)).*M(:,t),1);
+    [alpha(:,t),energy(t)] = normalize((At*alpha(:,t-1)).*M(:,t),1);    % 13.59
 end
 energy = sum(log(energy(energy>0)));

chapter13/HMM/hmmSmoother_.m

@@ -7,11 +7,11 @@
 alpha = zeros(K,T);
 [alpha(:,1),c(1)] = normalize(s.*M(:,1),1);
 for t = 2:T
-    [alpha(:,t),c(t)] = normalize((At*alpha(:,t-1)).*M(:,t),1);
+    [alpha(:,t),c(t)] = normalize((At*alpha(:,t-1)).*M(:,t),1);  % 13.59
 end
 beta = ones(K,T);
 for t = T-1:-1:1
-    beta(:,t) = A*(beta(:,t+1).*M(:,t+1))/c(t+1);
+    beta(:,t) = A*(beta(:,t+1).*M(:,t+1))/c(t+1);   % 13.62
 end
 gamma = alpha.*beta;
 

chapter13/HMM/hmmViterbi_.m

@@ -12,7 +12,7 @@
 Z(:,1) = 1:k;
 v = log(s(:))+M(:,1);
 for t = 2:n
-    [v,idx] = max(bsxfun(@plus,A,v),[],1);
+    [v,idx] = max(bsxfun(@plus,A,v),[],1);    % 13.68
     v = v(:)+M(:,t);
     Z = Z(idx,:);
     Z(:,t) = 1:k;