added comment to LDS

tcchriszhao · Jan 27, 2016 · 37a7f94 · 37a7f94
1 parent 62af89e
commit 37a7f94
Show file tree

Hide file tree

Showing 4 changed files with 29 additions and 27 deletions.
diff --git a/chapter13/LDS/demo.m b/chapter13/LDS/demo.m
@@ -2,11 +2,12 @@
 clear; close all;
 d = 3;
 k = 2;
-n = 100;
-
+n = 2000;
+% 
 [X,model] = ldsRnd(d,k,n);
+%% WARNING: The standard kalman filter as descripted in PRML is numerically unstable. Sometime, this demo fails, that is not my implementation problem.
 [mu, V, llh] = kalmanFilter(X, model);
 % % [nu, U, Ezz, Ezy, llh] = kalmanSmoother(X, model);
-[model, llh] = ldsEm(X, model);
+% [model, llh] = ldsEm(X, model);
 % plot(llh);
 
diff --git a/chapter13/LDS/kalmanFilter.m b/chapter13/LDS/kalmanFilter.m
@@ -17,9 +17,9 @@
 
 PC = P*C';
 R = (C*PC+S);
-K = PC/R;
-mu(:,1) = mu0+K*(X(:,1)-C*mu0);
-V(:,:,1) = (I-K*C)*P;
+K = PC/R;                                        % 13.97
+mu(:,1) = mu0+K*(X(:,1)-C*mu0);                     % 13.94
+V(:,:,1) = (I-K*C)*P;                               % 13.95
 llh(1) = logGauss(X(:,1),C*mu0,R);
 for i = 2:n
     [mu(:,i), V(:,:,i), llh(i)] = ...
@@ -28,12 +28,12 @@
 llh = sum(llh);
 
 function [mu, V, llh] = forwardStep(x, mu, V, A, G, C, S, I)
-P = A*V*A'+G;
-PC = P*C';
+P = A*V*A'+G;                                               % 13.88
+PC = P*C';                                                      
 R = C*PC+S;
-K = PC/R;
+K = PC/R;                                                   % 13.92
 Amu = A*mu;
-CAmu = C*Amu;
-mu = Amu+K*(x-CAmu);
-V = (I-K*C)*P;
-llh = logGauss(x,CAmu,R);
+CAmu = C*Amu;                                               
+mu = Amu+K*(x-CAmu);                                        % 13.89
+V = (I-K*C)*P;                                              % 13.90
+llh = logGauss(x,CAmu,R);                                   % 13.91
diff --git a/chapter13/LDS/kalmanSmoother.m b/chapter13/LDS/kalmanSmoother.m
@@ -46,19 +46,20 @@
 end
 
 function [mu, V, Amu, P, llh] = forwardStep(x, mu0, V0, A, G, C, S, I)
-P = A*V0*A'+G;
+P = A*V0*A'+G;                                              % 13.88
 PC = P*C';
 R = C*PC+S;
-K = PC/R;
+K = PC/R;                                                   % 13.92
 Amu = A*mu0;
 CAmu = C*Amu;
-mu = Amu+K*(x-CAmu);
-V = (I-K*C)*P;
-llh = logGauss(x,CAmu,R);
+mu = Amu+K*(x-CAmu);                                        % 13.89
+V = (I-K*C)*P;                                              % 13.90
+llh = logGauss(x,CAmu,R);                                   % 13.91
+
 
 function [nu, U, Ezz, Ezy] = backwardStep(nu0, U0, mu, V, Amu, P, A)
-J = V*A'/P;   % smoother gain matrix
-nu = mu+J*(nu0-Amu);
-U = V+J*(U0-P)*J';
-Ezy = J*U0+nu0*nu';     % E[z_t,z_{t+1}] eq 13.106
-Ezz = U+nu*nu';
+J = V*A'/P;                                                 % 13.102
+nu = mu+J*(nu0-Amu);                                        % 13.100
+U = V+J*(U0-P)*J';                                          % 13.101
+Ezy = J*U0+nu0*nu';                                         % 13.106 
+Ezz = U+nu*nu';                                             % 13.107
diff --git a/chapter13/LDS/ldsEm.m b/chapter13/LDS/ldsEm.m
@@ -23,13 +23,13 @@
 Ezz2 = Ezzn-Ezz(:,:,1);
 Ezy = sum(Ezy,3);
 
-A = Ezy/Ezz1;
+A = Ezy/Ezz1;                                           % 13.113
 EzyA = Ezy*A';
-G = (Ezz2-(EzyA+EzyA')+A*Ezz1*A')/(n-1);
+G = (Ezz2-(EzyA+EzyA')+A*Ezz1*A')/(n-1);                % 13.114
 Xnu = X*nu';
-C = Xnu/Ezzn;
+C = Xnu/Ezzn;                                           % 13.115
 XnuC = Xnu*C';
-S = (X*X'-(XnuC+XnuC')+C*Ezzn*C')/n;
+S = (X*X'-(XnuC+XnuC')+C*Ezzn*C')/n;                    % 13.116
 
 model.A = A;
 model.G = G;