%SIRLINEAR2: MATLAB script M-file in which SIR parameter
%values a and b are approximated by a linear fit based
%on the central difference derivative approximation. The
%susceptible and infected/infective equations are both 
%used. 
%
%DATA FROM THE BRITISH MEDICAL JOURNAL, MAR. 4, 1978.
S=[762 740 650 400 250 120 80 50 20 18 15 13 10]; 
I=[1 20 80 220 300 260 240 190 120 80 20 5 2]; 
days=[0 3 4 5 6 7 8 9 10 11 12 13 14]; 
%CENTRAL DIFFERENCE FIT 
dSdt = (S(4:end)-S(2:end-2))/2;
YS = dSdt./S(3:end-1);
XS = I(3:end-1);
dIdt = (I(4:end)-I(2:end-2))/2;
YI = dIdt./I(3:end-1);
XI = S(3:end-1);
%Define a column of ones and a column of zeros
col1 = ones(size(YS))';
col0 = zeros(size(YS))';
%Weights for the dependent variables
w1 = std(YS);
w2 = std(YI);
%Linear fit with a design matrix
Y = [YS'/w1; YI'/w2];
F1 = [-XS' col0];
F2 = [XI' col1];
F = [F1/w1; F2/w2];
p = F\Y;
a = p(1) 
b = - p(2)
%Standard deviation
q = length(Y)-length(p);
s = sqrt(norm(Y-F*p)^2/q);
s1 = s*w1
s2 = s*w2
%
pause
sirrhs = @(t,y) [-a*y(1)*y(2);a*y(1)*y(2)-b*y(2)];
[t,y]=ode45(sirrhs,[0 14],[762;1]);
subplot(2,1,1)
plot(days,S,'o',t,y(:,1))
title('Susceptible Population, model and data','Interpreter','latex','FontSize',16)
subplot(2,1,2)
plot(days,I,'o',t,y(:,2))
title('Infected Population, model and data','Interpreter','latex','FontSize',16)
xlabel('Days','Interpreter','latex','FontSize',14)
%
%To compute the standard deviation
[t,y]=ode45(sirrhs,days,[762;1]);
lE = norm(S'-y(:,1))^2+norm(I'-y(:,2))^2;
q = 2*(length(S)-1)-2;
ls = sqrt(lE/q)