wolffd@0
|
1 function [x,y] = sample_lds(F, H, Q, R, init_state, T, models, G, u)
|
wolffd@0
|
2 % SAMPLE_LDS Simulate a run of a (switching) stochastic linear dynamical system.
|
wolffd@0
|
3 % [x,y] = switching_lds_draw(F, H, Q, R, init_state, models, G, u)
|
wolffd@0
|
4 %
|
wolffd@0
|
5 % x(t+1) = F*x(t) + G*u(t) + w(t), w ~ N(0, Q), x(0) = init_state
|
wolffd@0
|
6 % y(t) = H*x(t) + v(t), v ~ N(0, R)
|
wolffd@0
|
7 %
|
wolffd@0
|
8 % Input:
|
wolffd@0
|
9 % F(:,:,i) - the transition matrix for the i'th model
|
wolffd@0
|
10 % H(:,:,i) - the observation matrix for the i'th model
|
wolffd@0
|
11 % Q(:,:,i) - the transition covariance for the i'th model
|
wolffd@0
|
12 % R(:,:,i) - the observation covariance for the i'th model
|
wolffd@0
|
13 % init_state(:,i) - the initial mean for the i'th model
|
wolffd@0
|
14 % T - the num. time steps to run for
|
wolffd@0
|
15 %
|
wolffd@0
|
16 % Optional inputs:
|
wolffd@0
|
17 % models(t) - which model to use at time t. Default = ones(1,T)
|
wolffd@0
|
18 % G(:,:,i) - the input matrix for the i'th model. Default = 0.
|
wolffd@0
|
19 % u(:,t) - the input vector at time t. Default = zeros(1,T)
|
wolffd@0
|
20 %
|
wolffd@0
|
21 % Output:
|
wolffd@0
|
22 % x(:,t) - the hidden state vector at time t.
|
wolffd@0
|
23 % y(:,t) - the observation vector at time t.
|
wolffd@0
|
24
|
wolffd@0
|
25
|
wolffd@0
|
26 if ~iscell(F)
|
wolffd@0
|
27 F = num2cell(F, [1 2]);
|
wolffd@0
|
28 H = num2cell(H, [1 2]);
|
wolffd@0
|
29 Q = num2cell(Q, [1 2]);
|
wolffd@0
|
30 R = num2cell(R, [1 2]);
|
wolffd@0
|
31 end
|
wolffd@0
|
32
|
wolffd@0
|
33 M = length(F);
|
wolffd@0
|
34 %T = length(models);
|
wolffd@0
|
35
|
wolffd@0
|
36 if nargin < 7,
|
wolffd@0
|
37 models = ones(1,T);
|
wolffd@0
|
38 end
|
wolffd@0
|
39 if nargin < 8,
|
wolffd@0
|
40 G = num2cell(repmat(0, [1 1 M]));
|
wolffd@0
|
41 u = zeros(1,T);
|
wolffd@0
|
42 end
|
wolffd@0
|
43
|
wolffd@0
|
44 [os ss] = size(H{1});
|
wolffd@0
|
45 state_noise_samples = cell(1,M);
|
wolffd@0
|
46 obs_noise_samples = cell(1,M);
|
wolffd@0
|
47 for i=1:M
|
wolffd@0
|
48 state_noise_samples{i} = sample_gaussian(zeros(length(Q{i}),1), Q{i}, T)';
|
wolffd@0
|
49 obs_noise_samples{i} = sample_gaussian(zeros(length(R{i}),1), R{i}, T)';
|
wolffd@0
|
50 end
|
wolffd@0
|
51
|
wolffd@0
|
52 x = zeros(ss, T);
|
wolffd@0
|
53 y = zeros(os, T);
|
wolffd@0
|
54
|
wolffd@0
|
55 m = models(1);
|
wolffd@0
|
56 x(:,1) = init_state(:,m);
|
wolffd@0
|
57 y(:,1) = H{m}*x(:,1) + obs_noise_samples{m}(:,1);
|
wolffd@0
|
58
|
wolffd@0
|
59 for t=2:T
|
wolffd@0
|
60 m = models(t);
|
wolffd@0
|
61 x(:,t) = F{m}*x(:,t-1) + G{m}*u(:,t-1) + state_noise_samples{m}(:,t);
|
wolffd@0
|
62 y(:,t) = H{m}*x(:,t) + obs_noise_samples{m}(:,t);
|
wolffd@0
|
63 end
|
wolffd@0
|
64
|
wolffd@0
|
65
|