%SIRLINEAR1: 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 and 
%the infected/infective equation from the first-order 
%SIR system. 
%
%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 
dIdt = (I(4:end)-I(2:end-2))/2;
Y = dIdt./I(3:end-1);
X = S(3:end-1);
plot(X,Y,'o');
title('Plot of $\frac{1}{I} \frac{dI}{dt}$ vs. $S$','Interpreter','latex','FontSize',16)
xlabel('$S$','Interpreter','latex','FontSize',14)
ylabel('$\frac{1}{I} \frac{dI}{dt}$','Interpreter','latex','FontSize',14)
pause
c = polyfit(X,Y,1);
a = c(1)
b = -c(2)
plot(X,Y,'o',X,c(1)*X+c(2)) 
title('Plot of $\frac{1}{I} \frac{dI}{dt}$ vs. $S$, with regression line', ...
    'Interpreter','latex','FontSize',16)
xlabel('$S$','Interpreter','latex','FontSize',14)
ylabel('$\frac{1}{I} \frac{dI}{dt}$','Interpreter','latex','FontSize',14)
%
pause
hold off
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)
