Least square estimation of equation dIdt

SIRC.m

@@ -44,25 +44,31 @@ y0 = [S0, I0, R0, C0];
 
 [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');
+figure(3)
+plot(t_noisy, y_noisy);
+xlabel('Time(days)');
+ylabel('Number of individuals');
+legend('S', 'I', 'R', 'C');
 
 %% Least square estimation
 
 % 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);
+%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);
+
+% 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);
 
-% dCdt = gamma*q*I - (Gamma-mu)*C;
+% 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
+[nu_lsq, mu_lsq, beta_lsq, epsilon_lsq, gamma_lsq, q_lsq, Gamma_lsq] = 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