close all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% P. M. Berthouex and L. C. Brown:
% Statistics for Environmental Engineers_, CRC Press, 2002.
%
% Fit the Monod model
%
%  y = theta_1 * x /(theta_2 + x)  + eps,  epsilon ~ N(0,I\sigma^2) 
%
% to observations
%
%   x (mg / L COD):  28    55    83    110   138   225   375   
%   y (1 / h):       0.053 0.060 0.112 0.105 0.099 0.122 0.125
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xdata = [ 28    55    83    110   138   225   375  ];
ydata = [ 0.053 0.060 0.112 0.105 0.099 0.122 0.125];

%first plot the data,  
plot(xdata,ydata,'o-');

% the value of the Monod model for large 'x' is close to 'theta_1', so form
% the plot get the initial guess theta_1=0.1;
% for theta_2 just take something of the size of 'xdata'. So the initial guess:

par_ini = [.1 100];
% call of the optimizer:
[par0,rss]   = fminsearch('monodss',par_ini,[],xdata,ydata);

% get the estimate for sigma2:
nobs = length(xdata);
npar = 2;                            %  n of parameters,th1 ja th2
sigma2 = rss/(nobs-npar);            %  variance of meas.error

%As for choice of the proposal, try the following:
cov_choice = input(' Enter 1 for identity covariance and 0 for chain updated covariance. ');
if cov_choice == 1
  qcov   = 1e-4*eye(npar);              %  covariance for Gaussian proposal
else
  qcov   = cov(chain);                 %  covariance for Gaussian proposal
end

%%%%%%%%%%%%%%%MCMC sampling starts here %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
R      = chol(qcov);
oldpar = par0;
ss     = monodss(oldpar,xdata,ydata); % first SS value
rej        = 0;                     % initialized count for rejections
nsimu  = 10000;                     % length of chain
chain = zeros(nsimu,length(par0));
chain(1,:) = oldpar;
for i=2:nsimu % simulation loop
   newpar = oldpar+randn(1,npar)*R;     % new candidate
   ss2    = ss;                         % old SS
   ss1    = monodss(newpar,xdata,ydata);  % new SS
   if (ss1<ss2 | rand < exp(-0.5*(ss1-ss2)/sigma2))
     chain(i,:) = newpar;              % accept
     oldpar     = newpar;
     ss         = ss1;   
   else   
     chain(i,:) = oldpar;              % reject
     rej        = rej+1;  
   end  
end
accept = 1-rej./nsimu;                  % acceptance rate

%%%%%%%%%%%MCMC sampling done. Now plot.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(2), subplot(211), plot(chain(:,1)), subplot(212), plot(chain(:,2))
sample_plot(chain',par0,3)