-
Notifications
You must be signed in to change notification settings - Fork 4
/
Copy pathSimulateLDS.m
93 lines (86 loc) · 2.39 KB
/
SimulateLDS.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
function [X, Y, Y_est] = SimulateLDS(LDS, T, nTrial)
% Simulate a Linear Dynamical System
%
%
% INPUT:
%
% LDS is a Linear Dynamical System model with fields:
%
% A,C,Q,R,x0,V0
%
% Corresponding to the Model:
%
% x(t+1) = A*x(t) + q(t)
% y(t) = C*x(t) + r(t)
%
% cov([q,r])=[Q 0; 0 R]
%
% x0 initial state estimate as column vector - (nx)x1
% (or a cell array with a vector for each "run")
% V0 cov of the initial state estimate - (nx)x(nx)
%
%
% T Number of time steps
% for one experiment:
% T is scalar
% for one or more experiments:
% T is vector of length E containing scalar trials counts
%
%
% OUTPUT
%
% X State
% Y Output
%
% For one experiment:
% X - (nx)xT
% Y - (ny)xT
% For more than one experiment:
% X,Y are cell arrays of length E containing
% X - (nx)xT(e)
% Y - (ny)xT(e)
% Copyright (C) 2005 Philip N. Sabes
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
%
%
% Modified by Ziqiang Wei
%
% 07/22/2014
%%%%%%%%%%% distribute LDS elements
FN = fieldnames(LDS);
for i=1:length(FN),
eval(sprintf('%s = LDS.%s;',FN{i},FN{i}));
end
%%%%%%%%%%% Dimensions
E = nTrial;
ny = size(C,1);
nx = size(A,1); % k= dim of states
sqQ = sqrtm(Q);
sqR = sqrtm(R);
sqV = sqrtm(V0);
X = zeros(nx, T, nTrial);
Y = zeros(ny, T, nTrial);
Y_est = zeros(ny, T, nTrial);
x = LDS.x0(:,ones(nTrial,1)) + sqV*randn(nx,nTrial);
for t = 1:T
y = C*x + sqR*randn(ny,nTrial); % output eq
y_est = C*x;
X(:,t,:) = x;
Y(:,t,:) = y;
Y_est(:,t,:) = y_est;
x = A*x + sqQ*randn(nx,nTrial);
end