% aidsrun.m
clear all, close all
%  Input aids data logistic model.
load aids
tvec = aids.year(13:24);%(1:12);
yvec = aids.cases(13:24);%(1:12);
nlin_fn = @logistic;
%b=[260 47572 1]'; 
b=[106618 42514 0.1]';
nparams = length(b);

% First run the optimzation routine
[b_opt,r,J,CovDelta_b,mse] = nlinfit(tvec,yvec,@logistic,b);
model_opt = feval(nlin_fn,b_opt,tvec);
figure(1)
  plot(tvec,yvec,'o',tvec,model_opt,'*')
  legend('data','model')

% Compute MCMC chain
params = {
    {'X0',b_opt(1), 0}
    {'M', b_opt(2), 0}
    {'k', b_opt(3), 0}
    };
data.tvec = tvec;
data.yvec = yvec;
model.ssfun  = @ss_logistic;
model.sigma2 = mse; 
model.N = length(yvec);
options.nsimu = 10000;
options.qcov = inv(J'*J)*mse*2.4^2./nparams;
% Sample the noise also. In particular the sigma^2 in N(0,sigma^2)
options.updatesigma = 1;    % to activate sigma2 sampling
model.N0            = 10;   % a prior parameter for 'sigma2'. Try n-p?
               % Small (0 ... 1) N0:large uncertainty in estimated |sigma2|
               % Large (10 --> ) N0:small uncertainty in estimated |sigma2|
[res,chain,s2chain] = mcmcrun(model,data,params,options);

% Plot the samples, compute confidence intervals, and model predictions
figure(2); clf
  mcmcplot(chain,[],res,'chainpanel');
figure(3); clf
  mcmcplot(chain,[],res,'pairs');
figure(4); clf
  mcmcplot(sqrt(s2chain),[],[],'hist')
  title('std of error, posterior distribution histogram')
out = mcmcpred(res,chain,s2chain,tvec,@logistic_switch);
obs{1}=[tvec yvec];
opts.states=[]; % no states plotted, just the response
mcmcpredplots(out,obs,opts);
title('Data vs distributions for model(dark grey) and noise(light grey)')

x = linspace(tvec(1),tvec(length(tvec))+10,100); 
x = x'; %note that x has to be column vector for |mcmcpred| !
out = mcmcpred(res,chain,s2chain,x,@logistic_switch);
mcmcpredplots(out,obs);
hold on;plot(tvec,yvec,'r*'); 
title('Extrapolated model predictions')