diff --git a/SIRC.m b/SIRC.m index ab5e2e2dd52a289c42f9950666f2b0d17c19dbc7..d8200c20469ceb0ed9f8233910aa4bcb8f76d156 100644 --- a/SIRC.m +++ b/SIRC.m @@ -53,21 +53,21 @@ legend('S', 'I', 'R', 'C'); %% Least square estimation of the noise-free modell % dSdt + dIdt = nu - mu*S - (gamma+mu)*I -X_SI = [ones(size(y, 1)-1, 1), y(1:end-1, 1), y(1:end-1, 2)]; -Y_SI = y(2:end,1) - y(1:end-1,1) + y(2:end,2) - y(1:end-1,2); +X_SI1 = [ones(size(y, 1)-1, 1), y(1:end-1, 1), y(1:end-1, 2)]; +Y_SI1 = y(2:end,1) - y(1:end-1,1) + y(2:end,2) - y(1:end-1,2); %Y_SI2 = y_noisy(2:end,1) + y_noisy(2:end,2); -theta_SI = abs(lsq(X_SI, Y_SI)); +theta_SI = abs(lsq(X_SI1, Y_SI1)); %theta_SI2 = lsq(X_SI, Y_SI2); % dIdt = beta*I*S + epsilon*beta*C*S - (gamma+mu)*I -X_I = [y(1:end-1,1).*y(1:end-1,2), y(1:end-1,1).*y(1:end-1,4), y(1:end-1,2)]; -Y_I = y(2:end, 2) - y(1:end-1, 2); -theta_I = abs(lsq(X_I, Y_I)); +X_I1 = [y(1:end-1,1).*y(1:end-1,2), y(1:end-1,1).*y(1:end-1,4), y(1:end-1,2)]; +Y_I1 = y(2:end, 2) - y(1:end-1, 2); +theta_I = abs(lsq(X_I1, Y_I1)); % dCdt = gamma*q*I - (Gamma-mu)*C -X_C = [y(1:end-1, 2), y(1:end-1, 4)]; -Y_C = y(2:end, 4) - y(1:end-1, 4); -theta_C = abs(lsq(X_C, Y_C)); +X_C1 = [y(1:end-1, 2), y(1:end-1, 4)]; +Y_C1 = y(2:end, 4) - y(1:end-1, 4); +theta_C = abs(lsq(X_C1, Y_C1)); [nu_lsq1, mu_lsq1, beta_lsq1, epsilon_lsq1, gamma_lsq1, q_lsq1, Gamma_lsq1] = deal(theta_SI(1), ... theta_SI(2), theta_I(1), theta_I(2)/theta_I(1), theta_SI(3)-theta_SI(2), ... @@ -101,20 +101,60 @@ legend('S', 'I', 'R', 'C'); %% Least square estimation after adding noise to the modell % 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); -theta_SI = abs(lsq(X_SI, Y_SI)); +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); +theta_SI = abs(lsq(X_SI2, Y_SI2)); % dIdt = beta*I*S + epsilon*beta*C*S - (gamma+mu)*I -X_I = [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_I = y_noisy(2:end, 2) - y_noisy(1:end-1, 2); -theta_I = abs(lsq(X_I, Y_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); +theta_I = abs(lsq(X_I2, Y_I2)); % dCdt = gamma*q*I - (Gamma-mu)*C -X_C = [y_noisy(1:end-1, 2), y_noisy(1:end-1, 4)]; -Y_C = y_noisy(2:end, 4) - y_noisy(1:end-1, 4); -theta_C = abs(lsq(X_C, Y_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); +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)); \ No newline at end of file + theta_C(1)/(theta_SI(3)-theta_SI(2)), theta_C(2)+theta_SI(2)); + +%% Instrumental variable method + +% 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)]; + +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)); + +[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), ... + theta_C(1)/(theta_SI(3)-theta_SI(2)), theta_C(2)+theta_SI(2)); + +y_iv4 = zeros(size(y, 1), 4); +y_iv4(1,:) = [S0, I0, R0, C0]; +for i = 2:length(y_iv4) + % dSdt = nu - (beta*I + epsilon*beta*C)*S - mu*S + y_iv4(i,1) = y_iv4(i-1,1) + nu_iv4 - (beta_iv4*y_iv4(i-1,2) + ... + epsilon_iv4*beta_iv4*y_iv4(i-1,4))*y_iv4(i-1,1) - mu_iv4*y_iv4(i-1,1); + % dIdt = (beta*I + epsilon*beta*C)*S - gamma*I - mu*I + y_iv4(i,2) = y_iv4(i-1,2) + (beta_iv4*y_iv4(i-1,2) + ... + epsilon_iv4*beta_iv4*y_iv4(i-1,4))*y_iv4(i-1,1) - ... + gamma_iv4*y_iv4(i-1,2) - mu_iv4*y_iv4(i-1,2); + % dRdt = gamma*(1-q)*I + Gamma*C - mu*R + y_iv4(i,3) = y_iv4(i-1,3) + gamma_iv4*(1-q_lsq1)*y_iv4(i-1,2) + ... + Gamma_iv4*y_iv4(i-1,4) - mu_iv4*y_iv4(i-1,3); + % dCdt = gamma*q*I - Gamma*C - mu*C + y_iv4(i,4) = y_iv4(i-1,4) + gamma_iv4*q_iv4*y_iv4(i-1,2) - ... + Gamma_iv4*y_iv4(i-1,4) - mu_iv4*y_iv4(i-1,3); +end + +figure(5) +plot(t, y_iv4); +title('Instrumental variable method (noise-free modell)'); +xlabel('Time(days)'); +ylabel('Number of individuals'); +legend('S', 'I', 'R', 'C'); \ No newline at end of file