wolffd@0: % Make a point move in the 2D plane
wolffd@0: % State = (x y xdot ydot). We only observe (x y).
wolffd@0: % Generate data from this process, and try to learn the dynamics back.
wolffd@0: 
wolffd@0: % X(t+1) = F X(t) + noise(Q)
wolffd@0: % Y(t) = H X(t) + noise(R)
wolffd@0: 
wolffd@0: ss = 4; % state size
wolffd@0: os = 2; % observation size
wolffd@0: F = [1 0 1 0; 0 1 0 1; 0 0 1 0; 0 0 0 1]; 
wolffd@0: H = [1 0 0 0; 0 1 0 0];
wolffd@0: Q = 0.1*eye(ss);
wolffd@0: R = 1*eye(os);
wolffd@0: initx = [10 10 1 0]';
wolffd@0: initV = 10*eye(ss);
wolffd@0: 
wolffd@0: seed = 1;
wolffd@0: rand('state', seed);
wolffd@0: randn('state', seed);
wolffd@0: T = 100;
wolffd@0: [x,y] = sample_lds(F, H, Q, R, initx, T);
wolffd@0: 
wolffd@0: % Initializing the params to sensible values is crucial.
wolffd@0: % Here, we use the true values for everything except F and H,
wolffd@0: % which we initialize randomly (bad idea!)
wolffd@0: % Lack of identifiability means the learned params. are often far from the true ones.
wolffd@0: % All that EM guarantees is that the likelihood will increase.
wolffd@0: F1 = randn(ss,ss);
wolffd@0: H1 = randn(os,ss);
wolffd@0: Q1 = Q;
wolffd@0: R1 = R;
wolffd@0: initx1 = initx;
wolffd@0: initV1 = initV;
wolffd@0: max_iter = 10;
wolffd@0: [F2, H2, Q2, R2, initx2, initV2, LL] =  learn_kalman(y, F1, H1, Q1, R1, initx1, initV1, max_iter);
wolffd@0: