diff --git a/SIRC.asv b/SIRC.asv new file mode 100644 index 0000000000000000000000000000000000000000..5c1c37f860a6324d5242928f7cc152b8126c49d8 --- /dev/null +++ b/SIRC.asv @@ -0,0 +1,64 @@ +%% 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 7443c2425e015d086bde4bde50b80018fc3e6333..093d9e17fce12012d69fdeae67c482e1e2fb4ef7 100644 --- a/SIRC.m +++ b/SIRC.m @@ -24,28 +24,44 @@ time = 200; 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'); +% 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'); +% 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'); +% 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)]; + + +% 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 diff --git a/lsq.m b/lsq.m new file mode 100644 index 0000000000000000000000000000000000000000..1383b5ca98c45973679b42b454a0bdcf80007906 --- /dev/null +++ b/lsq.m @@ -0,0 +1,4 @@ +function params = lsq(X, Y) + params = (X'*X)\(X'*Y); +end +