% SIR Example Code Identifiability Practical
% Marisa Eisenberg 7-29-16 (marisae@umich.edu)
clear

%% Load Angola Cholera Data

CholeraData;

%% Setup

% A note---I often find it's convenient to make small functions that take
% your parameters/data/model output and return (for example) the initial 
% conditions or measurement equation for your model. 
% This is because the initial conditions and the measurement equation both 
% depend on the parameter k. So we'll make two functions:
%
%   x0fcn: a function that takes in the current parameter values and
%   returns the initial conditions, i.e. I(0) = data(1)*k, S(0) = 1- I(0),
%   R(0) = 0. We'll assume the parameter order is: beta, gamma, k.

x0fcn = @(params) [1-params(3)*data(1); params(3)*data(1); 0];

%   yfcn: a function that takes in the current model output and parameter
%   values, and returns the measured cases, y = I/k. Note that k =
%   params(3), and I = x(:,2).

yfcn = @(x,params) x(:,2)/params(3);

% Note the "params" and "x" here are just dummy variables--they will be
% filled once we call these functions! However, "data(1)" is actually
% referencing the CholeraData we just loaded. Be careful--if you change the
% data later, this will not be reflected unless you reload the x0 function!
% You can change this if you want, e.g. to make x0fcn = @(params,data) etc.


% Starting parameter values:
params = [0.4; 0.25; 1/80000];


%% Simulate and Plot the Model

[t,x] = ode45(@sirODE,times,x0fcn(params),[],params);
y = yfcn(x,params);

figure(1)
    set(gca,'LineWidth',1,'FontSize',16,'FontName','Arial')
    hold on
    plot(t,y,'b','LineWidth',2);
    plot(times,data,'ko','LineWidth',2);
    %legend('Initial Simulation with Starting Parameters'); 
    ylabel('Infected Population');  
    xlabel('Time (days)');

%% Parameter Estimation (OLS + Nelder Mead)

% Estimate the model parameters 
[paramests, fval] = fminsearch(@(p) sirCost(times,p,data,x0fcn,yfcn),params,optimset('Display','iter','MaxFunEvals',5000,'MaxIter',5000));
paramests = abs(paramests);

% Re-simulate the model with the final parameter estimates
[test,xest] = ode45(@sirODE,times,x0fcn(paramests),[],paramests);
yest = yfcn(xest,paramests);

% Plot results
figure(1)  % Model Fit
    set(gca,'LineWidth',1,'FontSize',16,'FontName','Arial')
    hold on
    plot(test,yest,'k','LineWidth',2);
    plot(times,data,'ko','LineWidth',2);
    %legend('Initial Simulation with Starting Parameters','Model with Parameter Estimates'); 
    ylabel('Infected Population');  
    xlabel('Time (days)');
    
figure(2)  %Residuals of the fit
    set(gca,'LineWidth',1,'FontSize',16,'FontName','Arial')
    hold on
    plot(test,(yest-data),'o');
    ylabel('Residuals');  
    xlabel('Time (days)'); 

    
%% Fisher Information Matrix & Confidence Bounds

% Note: you can also do this analytically using the sensitivity equations,
% as you saw in yesterday's lecture/practical! Feel free to use and modify
% that code instead if you prefer.

FIM = MiniFisher(times,paramests,x0fcn,@sirODE,yfcn);

%Calculate & print rank of FIM
rank(FIM)

%Calculate SEs & CVs, print CVs
SEs = sqrt(diag(((length(data) - length(paramests))/fval)*inv(FIM)));
CVs = SEs./paramests
    

%% Other useful things

%%%%% Generating Simulated Data Sets
% For Problem 5, it's handy to have the poissrnd and normrnd functions, 
% which you can use to generate noisy data with a given distribution, mean,
% variance, etc. For example you could do:

poissdata = poissrnd(yest);

% To generate a simulated data set that's given by a Poisson distribution
% with mean = the model output trajectory yest. Similarly, we could make a
% noisy data set where we assume a Gaussian with mean yest and variance,
% say 0.1*yest:

normdata = normrnd(yest,0.1*yest);


%%%%% Plotting scatterplots of your parameter estimates
% Chris wrote a nice function that generates scatterplots of columns--so if
% you do a bunch of parameter estimates and get a column of estimates for
% each parameter, this will generate pair-wise scatterplots for each. 

% For example, here's a made up set of 10 estimates of our 3 parameters 
% (you would actually do the parameter estimation to make this, but just to illustrate):

fakeparams = normrnd(repmat(paramests',50,1),0.001);

% The function is called pairsplot_2, and is posted on the website. The
% function inputs are: 1) array of parameter estimates, 2) list of labels 
% (see example below, and 3) a vector of "true" parameter values which will
% be plotted in a different color.

pairsplot_2(fakeparams,{'\beta','\gamma','k'},paramests)





