function fit_Bombay_plague_data
global times deaths N preds
% uses data crudely pulled off the Kermack/McKendrick figure...
% quite inaccurate and likely missing some points...

times=[0.92540765,1.7937949,3.0680332,4.0520263,6.0180006,7.060356,9.138358,...
    11.033225,12.066189,13.103177,14.068723,14.751514,15.098108,16.089928,...
17.133066,18.274591,19.07981,21.024542,22.026312,23.132172,23.956396,24.950676,...
25.997614,26.98351,27.388912,28.085121,28.548662,29.070456,29.938843];

deaths=[12.695271,16.181248,15.079244,20.871433,...
53.17172,54.362595,125.79737,293.89972,...
391.76483,448.19812,643.88617,769.35913,...
778.5795,703.8099,696.9446,869.6216,...
927.19684,582.00507,404.8067,348.45563,...
210.38066,110.29155,64.29619,50.523354,...
50.538906,37.905846,35.621853,29.88743,33.37341];

N=29;   % 29 data points

initial_guess=[890 0.2 3.4];

[params SS_min]=fminsearch(@error_sum_of_squares,initial_guess);

fprintf('a: %g, b: %g, c: %g\n',params(1),params(2),params(3));
fprintf('SS_min: %g\n',SS_min);

% call SS function one last time, so that preds vector contains the best
% fitting values
error_sum_of_squares(params)
preds_best_fitting=preds;

% curious... how does this stack up to Kermack & McKendrick curve?
params=[890 0.2 3.4];
error_sum_of_squares(params)
preds_Kermack_McKendrick=preds;

figure(1)
plot(times,deaths,'x')
hold on
plot(times,preds_best_fitting,'-r')
plot(times,preds_Kermack_McKendrick,'--k')
legend('Data','Best fitting curve','Kermack/McKendrick curve')
hold off

end

function SS=error_sum_of_squares(params)
global times deaths N preds

% fit a*sech^2(b*t-c)
% params=[a b c]

SS=0;
for j=1:N
   preds(j)=params(1)* (sech(params(2)*times(j)-params(3)))^2;
   SS=SS+(preds(j)-deaths(j))^2;
end

end