From 38db85e80473dd82aab86e981c1d8b16dab6c955 Mon Sep 17 00:00:00 2001
From: "onica.klaudia" <onica.klaudia@hallgato.ppke.hu>
Date: Wed, 2 Jun 2021 21:06:00 +0200
Subject: [PATCH] Estimation of all parameters with least square method
---
SIRC.asv | 64 --------------------------------------------------------
SIRC.m | 17 ++++++++-------
2 files changed, 9 insertions(+), 72 deletions(-)
delete mode 100644 SIRC.asv
diff --git a/SIRC.asv b/SIRC.asv
deleted file mode 100644
index 5c1c37f..0000000
--- a/SIRC.asv
+++ /dev/null
@@ -1,64 +0,0 @@
-%% 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
diff --git a/SIRC.m b/SIRC.m
index 093d9e1..82497f9 100644
--- a/SIRC.m
+++ b/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);
\ No newline at end of file
+theta_SI2 = lsq(X_SI, Y_SI2);
+
+% 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
--
GitLab