% An MCMC example case. 
%
% Here we fit the function (used to model, e.g., the Biological Oxygen Demand, BOD)
%
%  y = b_1 * (1 - exp(-b_2*x)) + eps,     eps ~ N(0,I\sigma^2) 
%
% to observations
%
%   tvec   1      3      5      7      9   
%   yvec   0.076  0.258  0.369  0.492  0.559
%
clear all, close all
%% First fit the model to the data
tvec = [1 3 5 7 9]';  
yvec = [0.076  0.258  0.369  0.492  0.559]'; 
data.tvec     = tvec;  
data.yvec     = yvec;
modelfun      = @(tvec,theta) theta(1)*(1-exp(-theta(2)*tvec));
ssfun         = @(theta,data) sum((data.yvec-modelfun(data.tvec,theta)).^2);
th_init       = [1 0.15]; % ... after some trial & error, by plotting various options
[thmin,ssmin] = fminsearch(ssfun,th_init,[],data)
figure(1);
  plot(tvec,yvec,'s',tvec,modelfun(tvec,thmin));
  title('Data vs the fitted model'); 

%% Next, get things ready for sampling 
n            = length(data.tvec);
p            = 2;
mse          = ssmin/(n-p); % estimate for the error variance
J            = jacob(modelfun,data.tvec,thmin);
init_cov     = inv(J'*J)*mse;
% name, starting value for chain, minimum value, maximum value
params       = {
    {'theta1', thmin(1), 0, Inf}
    {'theta2', thmin(2), 0, Inf}
               };
model.ssfun  = ssfun;
model.sigma2 = mse; 
model.N      = n;
options.nsimu = 10000;
options.qcov = init_cov;

%% Now Sample
%(i) no sigma^2 sampling
[res,chain] = mcmcrun(model,data,params,options);
s2chain = [];
%(ii) with sigma^2 sampling
%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);

%% Finally, plot information from the sampling run.
figure(2); clf
  mcmcplot(chain,[],res,'chainpanel');
figure(3); clf
  mcmcplot(chain,[],res,'pairs');
% with sigma^2 sampling
%figure(4); clf
%  mcmcplot(sqrt(s2chain),[],[],'hist')
%  title('std of error, posterior distribution histogram')


out = mcmcpred(res,chain,s2chain,tvec,modelfun);
obs{1}=[data.tvec data.yvec];
mcmcpredplots(out,obs);
hold on;plot(tvec,yvec,'r*'); 
title('Data and possible model fits')

x = linspace(0,40,100); 
x = x'; %note that x has to be column vector for 'mcmcpred' !
out = mcmcpred(res,chain,s2chain,x,modelfun);
mcmcpredplots(out);
title('Extrapolated model predictions')