Skip to content
Snippets Groups Projects
Commit 9ad3c71d authored by Rivnyák Tímea's avatar Rivnyák Tímea
Browse files

SIRC model with recursive LSQ

parent 664b3848
No related branches found
No related tags found
No related merge requests found
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')
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment