diff --git a/SIRC_RLSQ.m b/SIRC_RLSQ.m index 7f33ccb4c9d992e456a922e11081c34225611e94..d7dd2b4c814cded9c7d10bf83e7ef8c845940b67 100644 --- a/SIRC_RLSQ.m +++ b/SIRC_RLSQ.m @@ -12,7 +12,7 @@ C0 = 0; S0 = N - I0 - R0 - C0; % Initial parameteres nu = 14; -mu = 0.015; +mu = 0.013; epsilon = 1.2; beta = 0.0002; gamma = 0.05; @@ -31,31 +31,11 @@ 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 +%% Recursive Least Square without noise h = 1.0; lambda0 = 0.95; % 0.95-0.99 -alpha = 0.7; +alpha = 0.5; N = length(y); % dSdt = nu - beta*I*S - epsilon*beta*C*S - mu*S; @@ -67,22 +47,21 @@ 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); +figure(2) +plot(theta_S); legend('nu', '-beta', '-epsilon*beta', '-mu') -% dIdt = beta*I*S + epsilon*beta*C*S - gamma*I - mu*I; +% dIdt = beta*I*S + epsilon*beta*C*S - (gamma + mu)*I; y_I = (y(2:end, 2) - y(1:end-1, 2)) / h; -phi_I = zeros(N-1, 4); +phi_I = zeros(N-1, 3); 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') +figure(3) +plot(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; @@ -92,19 +71,18 @@ 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); +figure(4) +plot(theta_R); legend('gamma*(1-q)', 'Gamma', '-mu') -% dCdt = gamma*q*I - Gamma*C - mu*C; +% dCdt = gamma*q*I - (Gamma + mu)*C; y_C = (y(2:end, 4) - y(1:end-1, 4)) / h; -phi_C = zeros(N-1, 3); +phi_C = zeros(N-1, 2); 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') +figure(5) +plot(theta_C); +legend('gamma*q', '-(Gamma+mu)')