image thumbnail

Hierarchical Kalman Filter for clinical time series prediction

by

 

It is an implementation of hierarchical (a.k.a. multi-scale) Kalman filter using belief propagation.

timeSeriesDataGeneration(useTrueParameter)
function ts_data = timeSeriesDataGeneration(useTrueParameter)
    ts_data.memo = ...
    {'% x{1}(1) ~ N(mu_0{1}, V_0{1});';...
     '% x{1}(t1) = A{1}*x{1}(t1-1) + q{1}; N(0, V_Q{1})';...
     '% y{1}(t1) = B{1}*x{1}(t1) + w{1}; N(0, V_W{1})';...
     '% x{2}(t2) = A{2}*x{1}(parent) + q{2}; N(0, V_Q{2})';...
     '% y{2}(t2) = B{2}*x{1}(t2) + w{2}; N(0, V_W{2})';};
    
    %% generating fake data.
    ts_data.s = RandStream('mt19937ar','Seed',4);
    RandStream.setGlobalStream(ts_data.s);
    Lv1 = 1;
    Lv2 = 2;
    ts_data.T(Lv1) = 40;  % total number of data points in level 1
    ts_data.T(Lv2) = 200; % total number of data points in level 2
    ts_data.num_Levels = numel(ts_data.T); % total number of data points in level 2
    [ts_data.connectionMap{Lv1}.toChildren, ts_data.connectionMap{Lv2}.toParents] = linspaceInt(1, ts_data.T(2), ts_data.T(1)); % connection map between level 1 and level 2.
    ts_data.dim_x(Lv1) = 2; % dimension of hidden state at level 1
    ts_data.dim_x(Lv2) = 2; % dimension of hidden state at level 2
    ts_data.dim_y(Lv1) = 2;    % dimension of observations at level 1
    ts_data.dim_y(Lv2) = 2;    % dimension of observations at level 2
    for Lvi = 1:ts_data.num_Levels
        ts_data.x{Lvi} = zeros(ts_data.dim_x(Lvi), ts_data.T(Lvi));    
        ts_data.y{Lvi} = zeros(ts_data.dim_y(Lvi), ts_data.T(Lvi));
    end    
    %% for hidden state in level 1
    %%
    ts_data.mu_0_true{Lv1} = [4; 4];       % true parameter
    ts_data.V_0_true{Lv1}  = [1, 0; 0, 1]; % true parameter
    ts_data.mu_0{Lv1} = [2; 2];            % initial parameter
    ts_data.V_0{Lv1}  = [1.2, 0; 0, 1.2];      % initial parameter
    if(useTrueParameter)
        ts_data.mu_0{Lv1} = ts_data.mu_0_true{Lv1};            % initial parameter
        ts_data.V_0{Lv1}  = ts_data.V_0_true{Lv1};      % initial parameter
    end
    ts_data.mu_0_initial{Lv1} = ts_data.mu_0{Lv1};            % initial parameter
    ts_data.V_0_initial{Lv1}  = ts_data.V_0{Lv1};      % initial parameter
    %%    
    ts_data.V_Q_ture{Lv1}  = [0.1,  0; 0, 0.3];      % true parameter
    ts_data.A_ture{Lv1}    = [1, 0.04; 0.02,    1];  % true parameter  
    ts_data.V_Q{Lv1}  = [0.5,  0; 0, 0.5];           % initial parameter
    ts_data.A{Lv1}    = [1.2, 0.1; 0.1,    1.2];       % initial parameter
    if(useTrueParameter)
        ts_data.V_Q{Lv1}  = ts_data.V_Q_ture{Lv1};           % initial parameter
        ts_data.A{Lv1}    = ts_data.A_ture{Lv1};       % initial parameter
    end
    ts_data.V_Q_initial{Lv1}  = ts_data.V_Q{Lv1};           % initial parameter
    ts_data.A_initial{Lv1}    = ts_data.A{Lv1};       % initial parameter
        
    %% for observation in level 1    
    %%
    ts_data.V_W_ture{Lv1}  = [0.05, 0; 0, 0.1];  % true parameter
    ts_data.B_ture{Lv1}    = [ 1, 0.1; 0.2, 1];  % true parameter
    ts_data.V_W{Lv1}  = [0.05, 0; 0, 0.1];       % initial parameter
    ts_data.B{Lv1}    = [ 1, 0.1; 0.2, 1];       % initial parameter
    if(useTrueParameter)
        ts_data.V_W{Lv1}  = ts_data.V_W_ture{Lv1};       % initial parameter
        ts_data.B{Lv1}    = ts_data.B_ture{Lv1};       % initial parameter
    end
    ts_data.V_W_initial{Lv1}  = ts_data.V_W{Lv1};       % initial parameter
    ts_data.B_initial{Lv1}    = ts_data.B{Lv1};       % initial parameter
        
    %% for hidden state in level 2
    %%
    ts_data.V_Q_ture{Lv2} = [0.08,    0;  0,   0.05]; % true parameter
    ts_data.A_ture{Lv2}   = [-0.6, 0.05; 0.01, 0.4];  % true parameter
    ts_data.V_Q{Lv2} = [0.1,    0;  0,   0.1];      % initial parameter
    ts_data.A{Lv2}   = [0.1, 0.1; 0.1, 0.2];       % initial parameter
    if(useTrueParameter)
        ts_data.V_Q{Lv2} = ts_data.V_Q_ture{Lv2};      % initial parameter
        ts_data.A{Lv2}   = ts_data.A_ture{Lv2};       % initial parameter
    end
    ts_data.V_Q_initial{Lv2} = ts_data.V_Q{Lv2};      % initial parameter
    ts_data.A_initial{Lv2}   = ts_data.A{Lv2};       % initial parameter
    %% for observation in level 2
    ts_data.V_W_ture{Lv2}  = [0.07, 0; 0, 0.09];       % true parameter
    ts_data.B_ture{Lv2}    = [ 0.6, 0.3; 0.06, 1.1];   % true parameter
    ts_data.V_W{Lv2}  = [0.1, 0; 0, 0.1];            % initial parameter
    ts_data.B{Lv2}    = [ 0.1, 0; 0.1, 1];        % initial parameter
    if(useTrueParameter)
        ts_data.V_W{Lv2}  = ts_data.V_W_ture{Lv2};            % initial parameter
        ts_data.B{Lv2}    = ts_data.B_ture{Lv2};        % initial parameter
    end
    ts_data.V_W_initial{Lv2}  = ts_data.V_W{Lv2};            % initial parameter
    ts_data.B_initial{Lv2}    = ts_data.B{Lv2};        % initial parameter
    %% build simulation data for level 1
    ts_data.x{Lv1}(:, 1) = mvnrnd(ts_data.mu_0_true{Lv1}, ts_data.V_0_true{Lv1}, 1)';
    ts_data.y{Lv1}(:, 1) = ts_data.B_ture{Lv1} * ts_data.x{Lv1}(:, 1) + mvnrnd(zeros(1, ts_data.dim_y(Lv1)), ts_data.V_W_ture{Lv1}, 1)';
    ts_data.alpha{Lv1} = ts_data.y{Lv1}(:, 1) * ts_data.y{Lv1}(:, 1)';
    for t = 2:ts_data.T(Lv1)
        ts_data.x{Lv1}(:, t) = ts_data.A_ture{Lv1} * ts_data.x{Lv1}(:, t-1) + mvnrnd(zeros(1, ts_data.dim_x(Lv1)), ts_data.V_Q_ture{Lv1}, 1)';
        ts_data.y{Lv1}(:, t) = ts_data.B_ture{Lv1} * ts_data.x{Lv1}(:, t) + mvnrnd(zeros(1, ts_data.dim_y(Lv1)), ts_data.V_W_ture{Lv1}, 1)';
        ts_data.alpha{Lv1} = ts_data.alpha{Lv1} + ts_data.y{Lv1}(:, t) * ts_data.y{Lv1}(:, t)';
    end
    
    % build simulation data for level 2
    ts_data.alpha{Lv2} = 0;
    for i = 1:ts_data.T(Lv2)
        ts_data.x{Lv2}(:, i) = ts_data.A_ture{Lv2} * ts_data.x{Lv1}(:, ts_data.connectionMap{Lv2}.toParents(i)) + mvnrnd(zeros(1, ts_data.dim_x(Lv2)), ts_data.V_Q_ture{Lv2}, 1)';
        ts_data.y{Lv2}(:, i) = ts_data.B_ture{Lv2} * ts_data.x{Lv2}(:, i) + mvnrnd(zeros(1, ts_data.dim_y(Lv2)), ts_data.V_W_ture{Lv2}, 1)';
        ts_data.alpha{Lv2} = ts_data.alpha{Lv2} + ts_data.y{Lv2}(:, i) * ts_data.y{Lv2}(:, i)';
    end
    figure(1); clf;
    subplot(1,2,1)
    plot(ts_data.x{Lv1}(1, :), ts_data.x{Lv1}(2, :), 'b*'); hold on;
    plot(ts_data.y{Lv1}(1, :), ts_data.y{Lv1}(2, :), 'kx'); hold on;
    legend('hidden state', 'observation');
    subplot(1,2,2)
    plot(ts_data.x{Lv2}(1, :), ts_data.x{Lv2}(2, :), 'b*'); hold on;
    plot(ts_data.y{Lv2}(1, :), ts_data.y{Lv2}(2, :), 'kx'); hold on;
    legend('hidden state', 'observation');
    ts_data.Lv1 = Lv1;
    ts_data.Lv2 = Lv2;
    ts_data.newLoglik = 0;
    ts_data.oldLoglik = -1000;
end

Contact us