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

SIRC.m

Blame
  • SIRC.m 7.73 KiB
    %% 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.015;
    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);
    
    % figure(1)
    % plot(t, y);
    % xlabel('Time(days)');
    % ylabel('Number of individuals');
    % legend('S', 'I', 'R', 'C');
    
    
    %% Add noise
    noise = wgn(size(t,1), 4, 0, 42);
    y_noisy1 = y + noise;
    
    figure(2)
    plot(t, y_noisy1);
    xlabel('Time(days)');
    ylabel('Number of individuals');
    legend('S', 'I', 'R', 'C');
    
    
    [t_noisy, y_noisy2] = ode45(@(t, y) noisy_deriv(y, nu, mu, epsilon, beta, gamma, Gamma, q), [0 time], y0);
    
    figure(3)
    plot(t_noisy, y_noisy2);
    xlabel('Time(days)');
    ylabel('Number of individuals');
    legend('S', 'I', 'R', 'C');
    
    %% Least square estimation
    %%% 1. For the noise-free modell
    
    % dSdt + dIdt = nu - mu*S - (gamma+mu)*I
    X_SI1 = [ones(size(y, 1)-1, 1), y(1:end-1, 1), y(1:end-1, 2)];
    Y_SI1 = y(2:end,1) - y(1:end-1,1) + y(2:end,2) - y(1:end-1,2);
    %Y_SI2 = y_noisy(2:end,1) + y_noisy(2:end,2);
    theta_SI = abs(lsq(X_SI1, Y_SI1));
    %theta_SI2 = lsq(X_SI, Y_SI2);
    
    % dIdt = beta*I*S + epsilon*beta*C*S - (gamma+mu)*I
    X_I1 = [y(1:end-1,1).*y(1:end-1,2), y(1:end-1,1).*y(1:end-1,4), y(1:end-1,2)];
    Y_I1 = y(2:end, 2) - y(1:end-1, 2);
    theta_I = abs(lsq(X_I1, Y_I1));
    
    % dCdt = gamma*q*I - (Gamma-mu)*C
    X_C1 = [y(1:end-1, 2), y(1:end-1, 4)];
    Y_C1 = y(2:end, 4) - y(1:end-1, 4);
    theta_C = abs(lsq(X_C1, Y_C1));
    
    [nu_lsq1, mu_lsq1, beta_lsq1, epsilon_lsq1, gamma_lsq1, q_lsq1, Gamma_lsq1] = deal(theta_SI(1), ...
        theta_SI(2), theta_I(1), theta_I(2)/theta_I(1), theta_SI(3)-theta_SI(2), ...
        theta_C(1)/(theta_SI(3)-theta_SI(2)), theta_C(2)+theta_SI(2));
    
    y_lsq1 = zeros(size(y, 1), 4);
    y_lsq1(1,:) = [S0, I0, R0, C0];
    for i = 2:length(y_lsq1)
        % dSdt = nu - (beta*I + epsilon*beta*C)*S - mu*S
        y_lsq1(i,1) = y_lsq1(i-1,1) + nu_lsq1 - (beta_lsq1*y_lsq1(i-1,2) + ...
            epsilon_lsq1*beta_lsq1*y_lsq1(i-1,4))*y_lsq1(i-1,1) - mu_lsq1*y_lsq1(i-1,1);
        % dIdt = (beta*I + epsilon*beta*C)*S - gamma*I - mu*I
        y_lsq1(i,2) = y_lsq1(i-1,2) + (beta_lsq1*y_lsq1(i-1,2) + ...
            epsilon_lsq1*beta_lsq1*y_lsq1(i-1,4))*y_lsq1(i-1,1) - ...
            gamma_lsq1*y_lsq1(i-1,2) - mu_lsq1*y_lsq1(i-1,2);
        % dRdt = gamma*(1-q)*I + Gamma*C - mu*R
        y_lsq1(i,3) = y_lsq1(i-1,3) + gamma_lsq1*(1-q_lsq1)*y_lsq1(i-1,2) + ...
            Gamma_lsq1*y_lsq1(i-1,4) - mu_lsq1*y_lsq1(i-1,3);
        % dCdt = gamma*q*I - Gamma*C - mu*C
        y_lsq1(i,4) = y_lsq1(i-1,4) + gamma_lsq1*q_lsq1*y_lsq1(i-1,2) - ...
            Gamma_lsq1*y_lsq1(i-1,4) - mu_lsq1*y_lsq1(i-1,3);
    end
    
    figure(4)
    plot(t, y_lsq1);
    title('Least-square estimation of the noise-free modell');
    xlabel('Time(days)');
    ylabel('Number of individuals');
    legend('S', 'I', 'R', 'C');
    
    %%% Adding noise to the output
    
    % dSdt + dIdt = nu - mu*S - (gamma+mu)*I
    X_SI2 = [ones(size(y, 1)-1, 1), y_noisy1(1:end-1, 1), y_noisy1(1:end-1, 2)];
    Y_SI2 = y_noisy1(2:end,1) - y_noisy1(1:end-1,1) + y_noisy1(2:end,2) - y_noisy1(1:end-1,2);
    %Y_SI2 = y_noisy(2:end,1) + y_noisy(2:end,2);
    theta_SI = abs(lsq(X_SI2, Y_SI2));
    %theta_SI2 = lsq(X_SI, Y_SI2);
    
    % dIdt = beta*I*S + epsilon*beta*C*S - (gamma+mu)*I
    X_I2 = [y_noisy1(1:end-1,1).*y_noisy1(1:end-1,2), y_noisy1(1:end-1,1).*y_noisy1(1:end-1,4), y_noisy1(1:end-1,2)];
    Y_I2 = y_noisy1(2:end, 2) - y_noisy1(1:end-1, 2);
    theta_I = abs(lsq(X_I2, Y_I2));
    
    % dCdt = gamma*q*I - (Gamma-mu)*C
    X_C2 = [y_noisy1(1:end-1, 2), y_noisy1(1:end-1, 4)];
    Y_C2 = y_noisy1(2:end, 4) - y_noisy1(1:end-1, 4);
    theta_C = abs(lsq(X_C2, Y_C2));
    
    [nu_lsq2, mu_lsq2, beta_lsq2, epsilon_lsq2, gamma_lsq2, q_lsq2, Gamma_lsq2] = deal(theta_SI(1), ...
        theta_SI(2), theta_I(1), theta_I(2)/theta_I(1), theta_SI(3)-theta_SI(2), ...
        theta_C(1)/(theta_SI(3)-theta_SI(2)), theta_C(2)+theta_SI(2));
    
    y_lsq2 = zeros(size(y_noisy1, 1), 4);
    y_lsq2(1,:) = [S0, I0, R0, C0];
    for i = 2:length(y_lsq2)
        % dSdt = nu - (beta*I + epsilon*beta*C)*S - mu*S
        y_lsq2(i,1) = y_lsq2(i-1,1) + nu_lsq2 - (beta_lsq2*y_lsq2(i-1,2) + ...
            epsilon_lsq2*beta_lsq2*y_lsq2(i-1,4))*y_lsq2(i-1,1) - mu_lsq2*y_lsq2(i-1,1);
        % dIdt = (beta*I + epsilon*beta*C)*S - gamma*I - mu*I
        y_lsq2(i,2) = y_lsq2(i-1,2) + (beta_lsq2*y_lsq2(i-1,2) + ...
            epsilon_lsq2*beta_lsq2*y_lsq2(i-1,4))*y_lsq2(i-1,1) - ...
            gamma_lsq2*y_lsq2(i-1,2) - mu_lsq2*y_lsq2(i-1,2);
        % dRdt = gamma*(1-q)*I + Gamma*C - mu*R
        y_lsq2(i,3) = y_lsq2(i-1,3) + gamma_lsq2*(1-q_lsq2)*y_lsq2(i-1,2) + ...
            Gamma_lsq2*y_lsq2(i-1,4) - mu_lsq2*y_lsq2(i-1,3);
        % dCdt = gamma*q*I - Gamma*C - mu*C
        y_lsq2(i,4) = y_lsq2(i-1,4) + gamma_lsq2*q_lsq2*y_lsq2(i-1,2) - ...
            Gamma_lsq2*y_lsq2(i-1,4) - mu_lsq2*y_lsq2(i-1,3);
    end
    
    figure(5)
    plot(t, y_lsq2);
    title('Least-square estimation of the modell with additive Gaussian noise');
    xlabel('Time(days)');
    ylabel('Number of individuals');
    legend('S', 'I', 'R', 'C');
    
    %%% Adding noise to the original equations
    
    % dSdt + dIdt = nu - mu*S - (gamma+mu)*I
    X_SI3 = [ones(size(y_noisy2, 1)-1, 1), y_noisy2(1:end-1, 1), y_noisy2(1:end-1, 2)];
    Y_SI3 = y_noisy2(2:end,1) - y_noisy2(1:end-1,1) + y_noisy2(2:end,2) - y_noisy2(1:end-1,2);
    theta_SI = abs(lsq(X_SI3, Y_SI3));
    
    % dIdt = beta*I*S + epsilon*beta*C*S - (gamma+mu)*I
    X_I3 = [y_noisy2(1:end-1,1).*y_noisy2(1:end-1,2), y_noisy2(1:end-1,1).*y_noisy2(1:end-1,4), y_noisy2(1:end-1,2)];
    Y_I3 = y_noisy2(2:end, 2) - y_noisy2(1:end-1, 2);
    theta_I = abs(lsq(X_I3, Y_I3));
    
    % dCdt = gamma*q*I - (Gamma-mu)*C
    X_C3 = [y_noisy2(1:end-1, 2), y_noisy2(1:end-1, 4)];
    Y_C3 = y_noisy2(2:end, 4) - y_noisy2(1:end-1, 4);
    theta_C = abs(lsq(X_C3, Y_C3));
    
    [nu_lsq3, mu_lsq3, beta_lsq3, epsilon_lsq3, gamma_lsq3, q_lsq3, Gamma_lsq3] = deal(theta_SI(1), ...
        theta_SI(2), theta_I(1), theta_I(2)/theta_I(1), theta_SI(3)-theta_SI(2), ...
        theta_C(1)/(theta_SI(3)-theta_SI(2)), theta_C(2)+theta_SI(2));
    
    %% Instrumental variable method for the modell with additive Gaussian noise
    
    % Matrices composed from the instrumental variables
    KSZI_SI = [ones(size(y, 1)-1, 1), y_lsq2(1:end-1, 1), y_lsq2(1:end-1, 2)];
    KSZI_I = [y_lsq2(1:end-1,1).*y_lsq2(1:end-1,2), y_lsq2(1:end-1,1).*y_lsq2(1:end-1,4), y_lsq2(1:end-1,2)];
    KSZI_C = [y_lsq2(1:end-1, 2), y_lsq2(1:end-1, 4)];
    
    theta_SI = abs(iv4(X_SI2, Y_SI2, KSZI_SI));
    theta_I = abs(iv4(X_I2, Y_I2, KSZI_I));
    theta_C = abs(iv4(X_C2, Y_C2, KSZI_C));
    
    [nu_iv4, mu_iv4, beta_iv4, epsilon_iv4, gamma_iv4, q_iv4, Gamma_iv4] = deal(theta_SI(1), ...
        theta_SI(2), theta_I(1), theta_I(2)/theta_I(1), theta_SI(3)-theta_SI(2), ...
        theta_C(1)/(theta_SI(3)-theta_SI(2)), theta_C(2)+theta_SI(2));
    
    y_iv4 = zeros(size(y, 1), 4);
    y_iv4(1,:) = [S0, I0, R0, C0];
    for i = 2:length(y_iv4)
        % dSdt = nu - (beta*I + epsilon*beta*C)*S - mu*S
        y_iv4(i,1) = y_iv4(i-1,1) + nu_iv4 - (beta_iv4*y_iv4(i-1,2) + ...
            epsilon_iv4*beta_iv4*y_iv4(i-1,4))*y_iv4(i-1,1) - mu_iv4*y_iv4(i-1,1);
        % dIdt = (beta*I + epsilon*beta*C)*S - gamma*I - mu*I
        y_iv4(i,2) = y_iv4(i-1,2) + (beta_iv4*y_iv4(i-1,2) + ...
            epsilon_iv4*beta_iv4*y_iv4(i-1,4))*y_iv4(i-1,1) - ...
            gamma_iv4*y_iv4(i-1,2) - mu_iv4*y_iv4(i-1,2);
        % dRdt = gamma*(1-q)*I + Gamma*C - mu*R
        y_iv4(i,3) = y_iv4(i-1,3) + gamma_iv4*(1-q_lsq1)*y_iv4(i-1,2) + ...
            Gamma_iv4*y_iv4(i-1,4) - mu_iv4*y_iv4(i-1,3);
        % dCdt = gamma*q*I - Gamma*C - mu*C
        y_iv4(i,4) = y_iv4(i-1,4) + gamma_iv4*q_iv4*y_iv4(i-1,2) - ...
            Gamma_iv4*y_iv4(i-1,4) - mu_iv4*y_iv4(i-1,3);
    end
    
    figure(6)
    plot(t, y_iv4);
    title('Instrumental variable method for the SIRC modell with additive Gaussian noise');
    xlabel('Time(days)');
    ylabel('Number of individuals');
    legend('S', 'I', 'R', 'C');