sample_lds.m
上传用户:mozhenmi
上传日期:2008-02-18
资源大小:13k
文件大小:2k
源码类别:

其他小程序

开发平台:

Matlab

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