Estimation of all parameters with least square method

SIRC.asv

-%% 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)];
-% 
-% theta_s = lsq(X_S, y_noisy(2:end,1));
-
-% dSdt + dIdt = nu - mu*S - (gamma+mu)*I
-X_SI = [ones(size(y_noisy, 1)-1, 1), 
\ No newline at end of file

SIRC.m

@@ -52,16 +52,17 @@ y0 = [S0, I0, R0, C0];
 
 %% 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);
+theta_SI2 = lsq(X_SI, Y_SI2);
\ No newline at end of file
+
+% dCdt = gamma*q*I - (Gamma-mu)*C;
+X_C = [y_noisy(1:end-1, 2), y_noisy(1:end-1, 4)];
+Y_C = y_noisy(2:end, 4) - y_noisy(1:end-1, 4);
+theta_C = lsq(X_C, Y_C);
+
+[nu_lsq, mu_lsq, gamma_lsq, q_lsq, Gamma_lsq] = deal(theta_SI(1), theta_SI(2), ...
+    theta_SI(3)-theta_SI(2), theta_C(1)/(theta_SI(3)-theta_SI(2)), theta_C(2)+theta_SI(2));
\ No newline at end of file