function [z,x] = sample_lds(A, C, U, V, init_state, T )
% SAMPLE_LDS Simulate a run of a stochastic linear dynamical system.
% [z,x] = sample_lds( A, C, U, V, init_state, T )
% 
%   z(t+1) =  A*z(t) + w(t),  w ~ N(0, U),  z(0) = init_state
%   x(t)   =  C*z(t) + v(t),  v ~ N(0, V)
%
% Input:
% A(:,:) - the transition matrix
% C(:,:) - the observation matrix
% U(:,:) - the transition covariance
% V(:,:) - the observation covariance
% init_state(:) - the initial mean
% T - the num. time steps to run for
%
% Output:
% z(:,t)    - the hidden state vector at time t.
% x(:,t)    - the observation vector at time t.

[os ss] = size(C);

state_noise_samples = gsamp(zeros(length(U),1), U, T)';
obs_noise_samples   = gsamp(zeros(length(V),1), V, T)';

z = zeros(ss, T);
x = zeros(os, T);

z(:,1) = init_state(:);
x(:,1) = C*z(:,1) + obs_noise_samples(:,1);

for t=2:T
  z(:,t) = A*z(:,t-1) + state_noise_samples(:,t);
  x(:,t) = C*z(:,t)   + obs_noise_samples(:,t);
end