From 8e6a6a1f0ae48896c96ab8e006e7e70b910f9ae5 Mon Sep 17 00:00:00 2001
From: "onica.klaudia" <onica.klaudia@hallgato.ppke.hu>
Date: Thu, 3 Jun 2021 18:15:19 +0200
Subject: [PATCH] Least square estimation of equation dIdt
---
SIRC.m | 26 ++++++++++++++++----------
1 file changed, 16 insertions(+), 10 deletions(-)
diff --git a/SIRC.m b/SIRC.m
index 82497f9..89b8472 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
--
GitLab