rewrite ldsRnd and lds_demo

chapter13/LDS/ldsRnd.m

@@ -1,4 +1,4 @@
-function [X, Z, model] = ldsRnd(d, k, n)
+function [Z, X] = ldsRnd(model, n)
 % Generate a data sequence from linear dynamic system.
 % Input:
 %   d: dimension of data
@@ -8,25 +8,27 @@
 %   X: d x n data matrix
 %   model: model structure
 % Written by Mo Chen ([email protected]).
-A = randn(k,k);
-G = iwishrnd(eye(k),k);
-C = randn(d,k);
-S = iwishrnd(eye(d),d);
-mu0 = randn(k,1);
-P0 = iwishrnd(eye(k),k);
+mu0 = model.mu0;
+P0 = model.P0;
+A = model.A;
+G = model.G;
+C = model.C;
+S = model.S;
+
+k = size(G,1);
+d = size(S,1);
 
 X = zeros(d,n);
 Z = zeros(k,n);
-Z(:,1) = gaussRnd(mu0,P0);              % 13.80
+Z(:,1) = gaussRnd(mu0,P0);                  % 13.80
 X(:,1) = gaussRnd(C*Z(:,1),S);
 for i = 2:n
-    Z(:,i) = gaussRnd(A*Z(:,i-1),G);           % 13.75, 13.78
-    X(:,i) = gaussRnd(C*Z(:,i),S);      % 13.76, 13.79
+    Z(:,i) = gaussRnd(A*Z(:,i-1),G);        % 13.75, 13.78
+    X(:,i) = gaussRnd(C*Z(:,i),S);          % 13.76, 13.79
 end
-
+model.mu0 = mu0; % prior mean
+model.P0 = P0;  % prior covairance
 model.A = A; % transition matrix 
 model.G = G; % transition covariance
 model.C = C; % emission matrix
 model.S = S;  % emision covariance
-model.mu0 = mu0; % prior mean
-model.P0 = P0;  % prior covairance

demo/ch13/lds_demo.m

@@ -1,14 +1,68 @@
-% demos for LDS in ch13
+close all;
+%% generate data
+clear; 
+d = 2;
+k = 4;
+n = 50;
 
-clear; close all;
-d = 3;
-k = 2;
-n = 100;
+A = [1 0 1 0; 
+     0 1 0 1;
+     0 0 1 0;
+     0 0 0 1]; 
+G = 0.001*eye(k);
 
-[X,Z,model] = ldsRnd(d,k,n);
-[mu, V, llh] = kalmanFilter(model, X);
+C = [1 0 0 0;
+     0 1 0 0];
+S = eye(d);
 
-[nu, U, Ezz, Ezy, llh] = kalmanSmoother(model, X);
-% [model, llh] = ldsEm(X,k);
-% plot(llh);
-% 
+mu0 = [8; 10; 1; 0];
+P0 = eye(k);
+
+model.A = A;
+model.G = G;
+model.C = C;
+model.S = S;
+model.mu0 = mu0;
+model.P0 = P0;
+
+[z,x] = ldsRnd(model, n);
+figure;
+hold on
+plot(x(1,:), x(2,:), 'ro');
+plot(z(1,:), z(2,:), 'b*-');
+legend('observed', 'latent')
+axis equal
+hold off
+
+%% filter
+[mu, V, llh] = kalmanFilter(model, x);
+figure
+hold on
+plot(x(1,:), x(2,:), 'ro');
+plot(mu(1,:), mu(2,:), 'b*-');
+legend('observed', 'filtered')
+axis equal
+hold off
+
+%% smoother
+[nu, U, llh] = kalmanSmoother(model, x);
+figure
+hold on
+plot(x(1,:), x(2,:), 'ro');
+plot(nu(1,:), nu(2,:), 'b*-');
+legend('observed', 'smoothed')
+axis equal
+hold off
+
+%% EM
+[model, llh] = ldsEm(x,model);
+nu = kalmanSmoother(model, x);
+figure
+hold on
+plot(x(1,:), x(2,:), 'ro');
+plot(nu(1,:), nu(2,:), 'b*-');
+legend('observed', 'smoothed with fitted model')
+axis equal
+hold off
+figure;
+plot(llh);