annotate toolboxes/FullBNT-1.0.7/Kalman/sample_lds.m @ 0:e9a9cd732c1e tip

first hg version after svn
author wolffd
date Tue, 10 Feb 2015 15:05:51 +0000
parents
children
rev   line source
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