update LDS

chapter13/LDS/kalmanFilter.m

@@ -1,5 +1,7 @@
-function [mu, V, llh] = kalmanFilter(X, model)
-% Kalman filter 
+function [mu, V, llh] = kalmanFilter(model, X)
+% Kalman filter (forward algorithm for linear dynamic system)
+% NOTE: This is the exact implementation of the Kalman filter algorithm in PRML.
+% However, this algorithm is not practical. It is numerical unstable. 
 % Input:
 %   X: d x n data matrix
 %   model: model structure
@@ -23,7 +25,7 @@
 I = eye(k);
 
 PC = P*C';
-R = (C*PC+S);
+R = C*PC+S;
 K = PC/R;                                        % 13.97
 mu(:,1) = mu0+K*(X(:,1)-C*mu0);                     % 13.94
 V(:,:,1) = (I-K*C)*P;                               % 13.95

chapter13/LDS/kalmanSmoother.m

@@ -1,5 +1,7 @@
-function [nu, U, Ezz, Ezy, llh] = kalmanSmoother(X, model)
+function [nu, U, Ezz, Ezy, llh] = kalmanSmoother(model, X)
 % Kalman smoother (forward-backward algorithm for linear dynamic system)
+% NOTE: This is the exact implementation of the Kalman smoother algorithm in PRML.
+% However, this algorithm is not practical. It is numerical unstable. 
 % Input:
 %   X: d x n data matrix
 %   model: model structure
@@ -28,7 +30,7 @@
 
 % forward
 PC = P0*C';
-R = (C*PC+S);
+R = C*PC+S;
 K = PC/R;
 mu(:,1) = mu0+K*(X(:,1)-C*mu0);
 V(:,:,1) = (I-K*C)*P0;

chapter13/LDS/ldsEm.m

@@ -1,18 +1,33 @@
-function [model, llh] = ldsEm(X, model)
+function [model, llh] = ldsEm(X, init)
 % EM algorithm for parameter estimation of linear dynamic system.
+% NOTE: This is the exact implementation of the EM algorithm in PRML.
+% However, this algorithm is not practical. It is numerical unstable and 
+% there is too much redundant degree of freedom. 
 % Input:
 %   X: d x n data matrix
 %   model: prior model structure
 % Output:
 %   model: trained model structure
 %   llh: loglikelihood
 % Written by Mo Chen ([email protected]).
-tol = 1e-4;
+d = size(X,1);
+if isstruct(init)   % init with a model
+    model = init;
+elseif numel(init) == 1  % random init with latent k
+    k = init;
+    model.A = randn(k,k);
+    model.G = iwishrnd(eye(k),k);
+    model.C = randn(d,k);
+    model.S = iwishrnd(eye(d),d);
+    model.mu0 = randn(k,1);
+    model.P0 = iwishrnd(eye(k),k);
+end
+tol = 1e-2;
 maxIter = 100;
 llh = -inf(1,maxIter);
 for iter = 2:maxIter
 %     E-step
-    [nu, U, Ezz, Ezy, llh(iter)] = kalmanSmoother(X, model);
+    [nu, U, Ezz, Ezy, llh(iter)] = kalmanSmoother(model,X);
     if llh(iter)-llh(iter-1) < tol*abs(llh(iter-1)); break; end   % check likelihood for convergence
 %     M-step 
     model = maximization(X, nu, U, Ezz, Ezy);

demo/ch13/lds_demo.m

@@ -6,9 +6,9 @@
 n = 100;
 
 [X,Z,model] = ldsRnd(d,k,n);
-[mu, V, llh] = kalmanFilter(X, model);
-
-[nu, U, Ezz, Ezy, llh] = kalmanSmoother(X, model);
-[model, llh] = ldsEm(X, model);
-plot(llh);
+[mu, V, llh] = kalmanFilter(model, X);
 
+[nu, U, Ezz, Ezy, llh] = kalmanSmoother(model, X);
+% [model, llh] = ldsEm(X,k);
+% plot(llh);
+%