LSQ for additive Gaussian noise modell + IV4

SIRC.m

@@ -32,25 +32,26 @@ y0 = [S0, I0, R0, C0];
 
 
 %% 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_noisy1 = y + noise;
+
+figure(2)
+plot(t, y_noisy1);
+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);
+[t_noisy, y_noisy2] = ode45(@(t, y) noisy_deriv(y, nu, mu, epsilon, beta, gamma, Gamma, q), [0 time], y0);
 
 figure(3)
-plot(t_noisy, y_noisy);
+plot(t_noisy, y_noisy2);
 xlabel('Time(days)');
 ylabel('Number of individuals');
 legend('S', 'I', 'R', 'C');
 
-%% Least square estimation of the noise-free modell
+%% Least square estimation
+%%% 1. For the noise-free modell
 
 % dSdt + dIdt = nu - mu*S - (gamma+mu)*I
 X_SI1 = [ones(size(y, 1)-1, 1), y(1:end-1, 1), y(1:end-1, 2)];
@@ -98,37 +99,85 @@ xlabel('Time(days)');
 ylabel('Number of individuals');
 legend('S', 'I', 'R', 'C');
 
-%% Least square estimation after adding noise to the modell
+%%% Adding noise to the output
 
 % dSdt + dIdt = nu - mu*S - (gamma+mu)*I
-X_SI2 = [ones(size(y_noisy, 1)-1, 1), y_noisy(1:end-1, 1), y_noisy(1:end-1, 2)];
-Y_SI2 = y_noisy(2:end,1) - y_noisy(1:end-1,1) + y_noisy(2:end,2) - y_noisy(1:end-1,2);
+X_SI2 = [ones(size(y, 1)-1, 1), y_noisy1(1:end-1, 1), y_noisy1(1:end-1, 2)];
+Y_SI2 = y_noisy1(2:end,1) - y_noisy1(1:end-1,1) + y_noisy1(2:end,2) - y_noisy1(1:end-1,2);
+%Y_SI2 = y_noisy(2:end,1) + y_noisy(2:end,2);
 theta_SI = abs(lsq(X_SI2, Y_SI2));
+%theta_SI2 = lsq(X_SI, Y_SI2);
 
 % dIdt = beta*I*S + epsilon*beta*C*S - (gamma+mu)*I
-X_I2 = [y_noisy(1:end-1,1).*y_noisy(1:end-1,2), y_noisy(1:end-1,1).*y_noisy(1:end-1,4), y_noisy(1:end-1,2)];
-Y_I2 = y_noisy(2:end, 2) - y_noisy(1:end-1, 2);
+X_I2 = [y_noisy1(1:end-1,1).*y_noisy1(1:end-1,2), y_noisy1(1:end-1,1).*y_noisy1(1:end-1,4), y_noisy1(1:end-1,2)];
+Y_I2 = y_noisy1(2:end, 2) - y_noisy1(1:end-1, 2);
 theta_I = abs(lsq(X_I2, Y_I2));
 
 % dCdt = gamma*q*I - (Gamma-mu)*C
-X_C2 = [y_noisy(1:end-1, 2), y_noisy(1:end-1, 4)];
-Y_C2 = y_noisy(2:end, 4) - y_noisy(1:end-1, 4);
+X_C2 = [y_noisy1(1:end-1, 2), y_noisy1(1:end-1, 4)];
+Y_C2 = y_noisy1(2:end, 4) - y_noisy1(1:end-1, 4);
 theta_C = abs(lsq(X_C2, Y_C2));
 
 [nu_lsq2, mu_lsq2, beta_lsq2, epsilon_lsq2, gamma_lsq2, q_lsq2, Gamma_lsq2] = deal(theta_SI(1), ...
     theta_SI(2), theta_I(1), theta_I(2)/theta_I(1), theta_SI(3)-theta_SI(2), ...
     theta_C(1)/(theta_SI(3)-theta_SI(2)), theta_C(2)+theta_SI(2));
 
-%% Instrumental variable method
+y_lsq2 = zeros(size(y_noisy1, 1), 4);
+y_lsq2(1,:) = [S0, I0, R0, C0];
+for i = 2:length(y_lsq2)
+    % dSdt = nu - (beta*I + epsilon*beta*C)*S - mu*S
+    y_lsq2(i,1) = y_lsq2(i-1,1) + nu_lsq2 - (beta_lsq2*y_lsq2(i-1,2) + ...
+        epsilon_lsq2*beta_lsq2*y_lsq2(i-1,4))*y_lsq2(i-1,1) - mu_lsq2*y_lsq2(i-1,1);
+    % dIdt = (beta*I + epsilon*beta*C)*S - gamma*I - mu*I
+    y_lsq2(i,2) = y_lsq2(i-1,2) + (beta_lsq2*y_lsq2(i-1,2) + ...
+        epsilon_lsq2*beta_lsq2*y_lsq2(i-1,4))*y_lsq2(i-1,1) - ...
+        gamma_lsq2*y_lsq2(i-1,2) - mu_lsq2*y_lsq2(i-1,2);
+    % dRdt = gamma*(1-q)*I + Gamma*C - mu*R
+    y_lsq2(i,3) = y_lsq2(i-1,3) + gamma_lsq2*(1-q_lsq2)*y_lsq2(i-1,2) + ...
+        Gamma_lsq2*y_lsq2(i-1,4) - mu_lsq2*y_lsq2(i-1,3);
+    % dCdt = gamma*q*I - Gamma*C - mu*C
+    y_lsq2(i,4) = y_lsq2(i-1,4) + gamma_lsq2*q_lsq2*y_lsq2(i-1,2) - ...
+        Gamma_lsq2*y_lsq2(i-1,4) - mu_lsq2*y_lsq2(i-1,3);
+end
+
+figure(5)
+plot(t, y_lsq2);
+title('Least-square estimation of the modell with additive Gaussian noise');
+xlabel('Time(days)');
+ylabel('Number of individuals');
+legend('S', 'I', 'R', 'C');
+
+%%% Adding noise to the original equations
+
+% dSdt + dIdt = nu - mu*S - (gamma+mu)*I
+X_SI3 = [ones(size(y_noisy2, 1)-1, 1), y_noisy2(1:end-1, 1), y_noisy2(1:end-1, 2)];
+Y_SI3 = y_noisy2(2:end,1) - y_noisy2(1:end-1,1) + y_noisy2(2:end,2) - y_noisy2(1:end-1,2);
+theta_SI = abs(lsq(X_SI3, Y_SI3));
+
+% dIdt = beta*I*S + epsilon*beta*C*S - (gamma+mu)*I
+X_I3 = [y_noisy2(1:end-1,1).*y_noisy2(1:end-1,2), y_noisy2(1:end-1,1).*y_noisy2(1:end-1,4), y_noisy2(1:end-1,2)];
+Y_I3 = y_noisy2(2:end, 2) - y_noisy2(1:end-1, 2);
+theta_I = abs(lsq(X_I3, Y_I3));
+
+% dCdt = gamma*q*I - (Gamma-mu)*C
+X_C3 = [y_noisy2(1:end-1, 2), y_noisy2(1:end-1, 4)];
+Y_C3 = y_noisy2(2:end, 4) - y_noisy2(1:end-1, 4);
+theta_C = abs(lsq(X_C3, Y_C3));
+
+[nu_lsq3, mu_lsq3, beta_lsq3, epsilon_lsq3, gamma_lsq3, q_lsq3, Gamma_lsq3] = deal(theta_SI(1), ...
+    theta_SI(2), theta_I(1), theta_I(2)/theta_I(1), theta_SI(3)-theta_SI(2), ...
+    theta_C(1)/(theta_SI(3)-theta_SI(2)), theta_C(2)+theta_SI(2));
+
+%% Instrumental variable method for the modell with additive Gaussian noise
 
 % Matrices composed from the instrumental variables
-KSZI_SI = [ones(size(y, 1)-1, 1), y_lsq1(1:end-1, 1), y_lsq1(1:end-1, 2)];
-KSZI_I = [y_lsq1(1:end-1,1).*y_lsq1(1:end-1,2), y_lsq1(1:end-1,1).*y_lsq1(1:end-1,4), y_lsq1(1:end-1,2)];
-KSZI_C = [y_lsq1(1:end-1, 2), y_lsq1(1:end-1, 4)];
+KSZI_SI = [ones(size(y, 1)-1, 1), y_lsq2(1:end-1, 1), y_lsq2(1:end-1, 2)];
+KSZI_I = [y_lsq2(1:end-1,1).*y_lsq2(1:end-1,2), y_lsq2(1:end-1,1).*y_lsq2(1:end-1,4), y_lsq2(1:end-1,2)];
+KSZI_C = [y_lsq2(1:end-1, 2), y_lsq2(1:end-1, 4)];
 
-theta_SI = abs(iv4(X_SI1, Y_SI1, KSZI_SI));
-theta_I = abs(iv4(X_I1, Y_I1, KSZI_I));
-theta_C = abs(iv4(X_C1, Y_C1, KSZI_C));
+theta_SI = abs(iv4(X_SI2, Y_SI2, KSZI_SI));
+theta_I = abs(iv4(X_I2, Y_I2, KSZI_I));
+theta_C = abs(iv4(X_C2, Y_C2, KSZI_C));
 
 [nu_iv4, mu_iv4, beta_iv4, epsilon_iv4, gamma_iv4, q_iv4, Gamma_iv4] = deal(theta_SI(1), ...
     theta_SI(2), theta_I(1), theta_I(2)/theta_I(1), theta_SI(3)-theta_SI(2), ...
@@ -152,9 +201,9 @@ for i = 2:length(y_iv4)
         Gamma_iv4*y_iv4(i-1,4) - mu_iv4*y_iv4(i-1,3);
 end
 
-figure(5)
+figure(6)
 plot(t, y_iv4);
-title('Instrumental variable method (noise-free modell)');
+title('Instrumental variable method for the SIRC modell with additive Gaussian noise');
 xlabel('Time(days)');
 ylabel('Number of individuals');
 legend('S', 'I', 'R', 'C');
\ No newline at end of file