Least square estimation of the noise-free system

SIRC.m

@@ -50,25 +50,71 @@ xlabel('Time(days)');
 ylabel('Number of individuals');
 legend('S', 'I', 'R', 'C');
 
-%% Least square estimation
+%% Least square estimation of the noise-free modell
 
 % 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);
+X_SI = [ones(size(y, 1)-1, 1), y(1:end-1, 1), y(1:end-1, 2)];
+Y_SI = 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 = lsq(X_SI, Y_SI);
+theta_SI = abs(lsq(X_SI, Y_SI));
 %theta_SI2 = lsq(X_SI, Y_SI2);
 
+% dIdt = beta*I*S + epsilon*beta*C*S - (gamma+mu)*I
+X_I = [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_I = y(2:end, 2) - y(1:end-1, 2);
+theta_I = abs(lsq(X_I, Y_I));
+
+% dCdt = gamma*q*I - (Gamma-mu)*C
+X_C = [y(1:end-1, 2), y(1:end-1, 4)];
+Y_C = y(2:end, 4) - y(1:end-1, 4);
+theta_C = abs(lsq(X_C, Y_C));
+
+[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');
+
+%% Least square estimation after adding noise to the modell
+
+% 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);
+theta_SI = abs(lsq(X_SI, Y_SI));
+
 % dIdt = beta*I*S + epsilon*beta*C*S - (gamma+mu)*I
 X_I = [y_noisy(1:end-1,1).*y_noisy(1:end-1,2), y_noisy(1:end-1,1).*y_noisy(1:end-1,4), y_noisy(1:end-1,2)];
 Y_I = y_noisy(2:end, 2) - y_noisy(1:end-1, 2);
-theta_I = lsq(X_I, Y_I);
+theta_I = abs(lsq(X_I, Y_I));
 
 % 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);
+theta_C = abs(lsq(X_C, Y_C));
 
-[nu_lsq, mu_lsq, beta_lsq, epsilon_lsq, gamma_lsq, q_lsq, Gamma_lsq] = deal(theta_SI(1), ...
+[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));
\ No newline at end of file