Skip to content
Snippets Groups Projects
Select Git revision
  • c026d2c82e81e3d6642131a93b7bcce99229d7d5
  • master default protected
2 results

SIRC_noisyRLSQ.m

Blame
  • SIRC_noisyRLSQ.m 2.31 KiB
    clc
    %% System initalization
    % Population
    N = 1000;
    % Initial number of infected individuals
    I0 = 10;
    % Initial number of recovered individuals
    R0 = 0;
    % Initial number of carrier individuals
    C0 = 0;
    % Everyone else is susceptible
    S0 = N - I0 - R0 - C0;
    % Initial parameteres
    nu = 14;
    mu = 0.013;
    epsilon = 1.2;
    beta = 0.0002;
    gamma = 0.05;
    Gamma = 0.03;
    q = 0.8;
    
    time = 200;
    %t = linspace(0, time, time);
    
    y0 = [S0, I0, R0, C0];
    [t, y] = ode45(@(t, y) deriv(y, nu, mu, epsilon, beta, gamma, Gamma, q), [0 time], y0);
    
    
    %% Add noise
    
    [t_noisy, y_noisy] = ode45(@(t, y) noisy_deriv(y, nu, mu, epsilon, beta, gamma, Gamma, q), [0 time], y0);
    
    figure(1)
    plot(t_noisy, y_noisy);
    xlabel('Time(days)');
    ylabel('Number of individuals');
    legend('S', 'I', 'R', 'C');
    
    %% Recursive Least Square with ar additive noise
    
    h = 1.0;
    lambda0 = 0.99; % 0.95-0.99
    alpha = 0.8;
    N = length(y_noisy);
    
    % dSdt = nu - beta*I*S - epsilon*beta*C*S - mu*S;
    y_S = (y_noisy(2:end, 1) - y_noisy(1:end-1, 1)) / h;
    phi_S = zeros(N-1, 4);
    phi_S(:,1) = ones(N-1, 1);
    phi_S(:,2) = y_noisy(1:end-1, 2) .* y_noisy(1:end-1, 1);
    phi_S(:,3) = y_noisy(1:end-1, 4) .* y_noisy(1:end-1, 1);
    phi_S(:,4) = y_noisy(1:end-1, 1);
    
    [theta_S, p_S] = recursiveLSQ(phi_S, y_S, lambda0, alpha)
    figure(2)
    plot(theta_S);
    legend('nu', '-beta', '-epsilon*beta', '-mu')
    
    % dIdt = beta*I*S + epsilon*beta*C*S - (gamma + mu)*I;
    y_I = (y_noisy(2:end, 2) - y_noisy(1:end-1, 2)) / h;
    phi_I = zeros(N-1, 3);
    phi_I(:,1) = y_noisy(1:end-1, 2) .* y_noisy(1:end-1, 1);
    phi_I(:,2) = y_noisy(1:end-1, 4) .* y_noisy(1:end-1, 1);
    phi_I(:,3) = y_noisy(1:end-1, 2);
    
    [theta_I, p_I] = recursiveLSQ(phi_I, y_I, lambda0, alpha)
    figure(3)
    plot(theta_I);
    legend('beta', 'epsilon*beta', '-(gamma+mu)')
    
    % dRdt = gamma*(1-q)*I + Gamma*C - mu*R;
    y_R = (y_noisy(2:end, 3) - y_noisy(1:end-1, 3)) / h;
    phi_R = zeros(N-1, 3);
    phi_R(:,1) = y_noisy(1:end-1, 2);
    phi_R(:,2) = y_noisy(1:end-1, 4);
    phi_R(:,3) = y_noisy(1:end-1, 3);
    
    [theta_R, p_R] = recursiveLSQ(phi_R, y_R, lambda0, alpha)
    figure(4)
    plot(theta_R);
    legend('gamma*(1-q)', 'Gamma', '-mu')
    
    % dCdt = gamma*q*I - (Gamma + mu)*C;
    y_C = (y_noisy(2:end, 4) - y_noisy(1:end-1, 4)) / h;
    phi_C = zeros(N-1, 2);
    phi_C(:,1) = y_noisy(1:end-1, 2);
    phi_C(:,2) = y_noisy(1:end-1, 4);
    
    [theta_C, p_C] = recursiveLSQ(phi_C, y_C, lambda0, alpha)
    figure(5)
    plot(theta_C);
    legend('gamma*q', '-(Gamma+mu)')