From 9ad3c71d67bdb3cfb3f87e97bb4d10fe28e8e188 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rivny=C3=A1k=20T=C3=ADmea?= <rivnyak.timea@hallgato.ppke.hu> Date: Mon, 7 Jun 2021 15:34:19 +0000 Subject: [PATCH] SIRC model with recursive LSQ --- SIRC_RLSQ.m | 110 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 110 insertions(+) create mode 100644 SIRC_RLSQ.m diff --git a/SIRC_RLSQ.m b/SIRC_RLSQ.m new file mode 100644 index 0000000..7f33ccb --- /dev/null +++ b/SIRC_RLSQ.m @@ -0,0 +1,110 @@ +clc +%% 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'); + +%% Recursive Least Square + +h = 1.0; +lambda0 = 0.95; % 0.95-0.99 +alpha = 0.7; +N = length(y); + +% dSdt = nu - beta*I*S - epsilon*beta*C*S - mu*S; +y_S = (y(2:end, 1) - y(1:end-1, 1)) / h; +phi_S = zeros(N-1, 4); +phi_S(:,1) = ones(N-1, 1); +phi_S(:,2) = y(1:end-1, 2) .* y(1:end-1, 1); +phi_S(:,3) = y(1:end-1, 4) .* y(1:end-1, 1); +phi_S(:,4) = y(1:end-1, 1); + +[theta_S, p_S] = recursiveLSQ(phi_S, y_S, lambda0, alpha) +figure(4) +plot(1:size(theta_S, 1), theta_S); +legend('nu', '-beta', '-epsilon*beta', '-mu') + +% dIdt = beta*I*S + epsilon*beta*C*S - gamma*I - mu*I; +y_I = (y(2:end, 2) - y(1:end-1, 2)) / h; +phi_I = zeros(N-1, 4); +phi_I(:,1) = y(1:end-1, 2) .* y(1:end-1, 1); +phi_I(:,2) = y(1:end-1, 4) .* y(1:end-1, 1); +phi_I(:,3) = y(1:end-1, 2); +phi_I(:,4) = y(1:end-1, 2); + +[theta_I, p_I] = recursiveLSQ(phi_I, y_I, lambda0, alpha) +figure(5) +plot(1:size(theta_I, 1), theta_I); +legend('beta', 'epsilon*beta', '-gamma', '-mu') + +% dRdt = gamma*(1-q)*I + Gamma*C - mu*R; +y_R = (y(2:end, 3) - y(1:end-1, 3)) / h; +phi_R = zeros(N-1, 3); +phi_R(:,1) = y(1:end-1, 2); +phi_R(:,2) = y(1:end-1, 4); +phi_R(:,3) = y(1:end-1, 3); + +[theta_R, p_R] = recursiveLSQ(phi_R, y_R, lambda0, alpha) +figure(6) +plot(1:size(theta_R, 1), theta_R); +legend('gamma*(1-q)', 'Gamma', '-mu') + +% dCdt = gamma*q*I - Gamma*C - mu*C; +y_C = (y(2:end, 4) - y(1:end-1, 4)) / h; +phi_C = zeros(N-1, 3); +phi_C(:,1) = y(1:end-1, 2); +phi_C(:,2) = y(1:end-1, 4); +phi_C(:,3) = y(1:end-1, 4); + +[theta_C, p_C] = recursiveLSQ(phi_C, y_C, lambda0, alpha) +figure(7) +plot(1:size(theta_C, 1), theta_C); +legend('gamma*q', '-Gamma', '-mu') + -- GitLab