SIRC.m

%% 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_noise = y + noise;
% 
% figure(2)
% plot(t, y_noise);
% xlabel('Time(days)');
% ylabel('Number of individuals');
% legend('S', 'I', 'R', 'C');


[t_noisy, y_noisy] = ode45(@(t, y) noisy_deriv(y, nu, mu, epsilon, beta, gamma, Gamma, q), [0 time], y0);

% figure(3)
% plot(t_noisy, y_noisy);
% xlabel('Time(days)');
% ylabel('Number of individuals');
% legend('S', 'I', 'R', 'C');

%% Least square estimation

% The matrices of known properties
% X_S = [ones(size(y_noisy, 1)-1, 1), y_noisy(1:end-1, 2), y_noisy(1:end-1, 4), y_noisy(1:end-1, 1)];
% X_I = [y_noisy(1:end-1, 2), y_noisy(1:end-1, 4), y_noisy(1:end-1, 1)];
% X_R = [y_noisy(1:end-1, 2), y_noisy(1:end-1, 4), y_noisy(1:end-1, 3)];
% X_C = [y_noisy(1:end-1, 2), y_noisy(1:end-1, 4)];


% dSdt + dIdt = nu - mu*S - (gamma+mu)*I
X_SI = [ones(size(y_noisy, 1)-1, 1), y_noisy(1:end-1, 1), y_noisy(1:end-1, 2)];
Y_SI = y_noisy(2:end,1) - y_noisy(1:end-1,1) + y_noisy(2:end,2) - y_noisy(1:end-1,2);
Y_SI2 = y_noisy(2:end,1) + y_noisy(2:end,2);
theta_SI = lsq(X_SI, Y_SI);
theta_SI2 = lsq(X_SI, Y_SI2);