Least square estimation algorithm

SIRC.asv

+%% 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

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

lsq.m

+function params = lsq(X, Y)
+    params = (X'*X)\(X'*Y);
+end
+