function[chain,acceptance_rate]=Metropolis(M,param,C,sigsq,cost_fn,...
                                 cost_params,bounds,adapt_int)
% This function implements the Metropolis algorithm for sampling from the
% posterior distribution of the parameters in the param vector.
% M                 = the length of the Markov chain;
% param             = the initial parameter (row) vector in the chain;
% C                 = covariance of the Gaussian proposal;
% sigsq             = (estimated) sample variance.
% cost_fun          = the likelihood of the parameters given the data;
% cost_params       = parameters needed to evaluate cost_fn.
% bounds            = constraints for the parameter values. Row 1 contains 
%                     the lower bounds and row 2 the upper bounds.                      
% adapt_int         = interval for adaptive update of the covariance. 0 =
%                             no adaptation.
% First randomize the seed of the random number generator.
rand('state',sum(100*clock));
% Now implement the Metropolis algorithm.
nparams = length(param);
chain = zeros(M,nparams);
old_param = param;
chain(1,:) = old_param;
SS_old = feval(cost_fn,old_param,cost_params);
R = chol(C);
naccept = 0;
bound_flag = 0;
for i=1:M
  hh = waitbar(i/M);
  new_param = old_param + randn(1,nparams)*R;
  % Don't accept the step if it is outside the constraints or if the
  % objective function isn't sufficiently decreased.
  bound_flag = sum((new_param < bounds(1,:)) + (new_param > bounds(2,:)));
  if bound_flag == 0
    SS_new = feval(cost_fn,new_param,cost_params);
    if (min(1,exp(-(SS_new-SS_old)/(2*sigsq))) > rand)
      % Accept
      naccept = naccept + 1;
      old_param = new_param;
      SS_old = SS_new;
    end
  end
  chain(i,:) = old_param;
  if adapt_int > 0
    if i/adapt_int == fix(i/adapt_int)
      Cadapt = cov(chain(1:i,:));
      R = chol(Cadapt);
    end
  end
end
close(hh)
acceptance_rate = naccept/M;