diff --git a/SIRC.m b/SIRC.m index 82497f9c7e04271a3f01923e9df9628a0be63c31..89b84728fa67396218e225d977715152fb66757c 100644 --- a/SIRC.m +++ b/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