% % Example11_24.m % % Requires: SIR_fun2.m, SIR_rhs2.m SIR2.txt, SIR2_integer.txt % clear all global N % % Load data which consists of 141 points of infectious data along with the % associated times. % %load SIR2.txt load SIR2_integer.txt % % Set parameters and initial conditions. These values are close to the flu data % but with larger I0 and smaller variance to prevent negative values of I(t). % S0 = 760; I0 = 12; R0 = 0; N = S0 + I0 + R0; beta = 0.0022; % Estimated values for flu data: beta = 2.3040e-3; gamma = 0.4470; % gamma = 0.4284; beta_mean = beta; gamma_mean = gamma; params = [beta gamma]; h = 1e-16; tf = 14; dt = 0.1; t_data = SIR2_integer(:,1); % I_data = SIR2(:,2); I_data = SIR2_integer(:,2); t_data = 0:dt:tf; Y0 = [S0; I0]; n = length(t_data); p = 2; % % Specify the parameters and compute the solution I(t). When integrating % the system, we take advantage of the fact that R(t) = N - S(t) - I(t). % % Construct synthetic data having the observation variance sigma^2 = 81; % Note that this is significantly less than the observation variance % sigma^2 = 322.8 estimated for the flu data. % q_init = [beta gamma]; ode_options = odeset('RelTol',1e-8); [t,Y] = ode45(@SIR_rhs2,t_data,Y0,ode_options,q_init); Y(:,3) = N - Y(:,1) - Y(:,2); I = Y(:,2); sigma = 9; sigma2 = sigma^2; error = sigma*randn(size(t_data')); % I_data = I + error; % I_data_integer = round(I_data); % SIR_data = [t_data' I_data]; % SIR_data_integer = [t_data' I_data_integer]; % save('SIR2.txt','SIR_data','-ascii') % save('SIR2_integer.txt','SIR_data_integer','-ascii') Res = I_data - I; % % Plot the data and model % figure(1) plot(t,Y(:,1),'k-',t,Y(:,2),'r--',t,Y(:,3),'b-.','linewidth',3) axis([0 tf 0 800]) set(gca,'Fontsize',[22]); xlabel('Time (Days)') ylabel('Responses S(t), I(t), R(t)') legend('Susceptible S(t)','Infectious I(t)','Recovered R(t)','Location','East') %print -depsc SIR2_responses figure(2) plot(t,I,'k-',t,I_data,'bx',t,0*t,'k-','linewidth',3) axis([0 tf 0 350]) set(gca,'Fontsize',[22]); xlabel('Time (Days)') ylabel('Infectious I(t)') legend('Model','Data','Location','NorthEast') %print -depsc SIR2_model_data figure(3) plot(t,Res,'xb',t,0*t,'k-','linewidth',3) axis([0 tf -25 25]); hold on set(gca,'Fontsize',[22]); plot(t,2*sigma*ones(length(t)),'--k','linewidth',3) plot(t,-2*sigma*ones(length(t)),'--k','linewidth',3) hold off xlabel('Time (Days)') ylabel('Residuals') %print -depsc SIR2_residuals % % Estimate the parameters beta and gamma and compute the resulting trajectories I(t) % ode_options = odeset('RelTol',1e-6); q_init = [beta gamma]; modelfun = @(q)SIR_fun2(q,Y0,t_data,I_data); [q_opt,fval] = fminsearch(modelfun,q_init); beta_opt = q_opt(1) gamma_opt = q_opt(2) [t,Y] = ode45(@SIR_rhs2,t_data,Y0,ode_options,q_opt); Y(:,3) = N - Y(:,1) - Y(:,2); I = Y(:,2); % % Use the complex-step derivative approximation to compute the sensitivity matrix. Compute % the observation error estimate and V. % Y0 = [S0; I0]; for i=1:p if i==1 beta_complex = complex(beta_opt,h); q = [beta_complex,gamma]; [t,Ypert] = ode45(@SIR_rhs2,t_data,Y0,ode_options,q); I_beta = imag(Ypert(:,2))/h; Sens(:,i) = I_beta; else gamma_complex = complex(gamma_opt,h); q = [beta,gamma_complex]; [t,Ypert] = ode45(@SIR_rhs2,t_data,Y0,ode_options,q); I_gamma = imag(Ypert(:,2))/h; Sens(:,i) = I_gamma; end end Res = I_data - I; sigma2_est = (1/(n-p))*Res'*Res V = sigma2_est*inv(Sens'*Sens) %% % Construct the 95% confidence intervals. % alpha = 0.05; tval = tinv(1-alpha/2,n-p); int_beta = [beta_opt - sqrt(V(1,1))*tval beta_opt + sqrt(V(1,1))*tval] int_r = [gamma_opt - sqrt(V(2,2))*tval gamma_opt + sqrt(V(2,2))*tval] % %