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');