% simple Metropolis algorithm
% run in Matlab as "p = metropolis2"

function r = metropolis2

% parameter specific entries
 
par = [ 3 3 3 ];  % initial values for parameters 
LW = -100; % lower bounds for parameters 
UP = 100; % upper bounds for parameters

% algorithm specific entries

iter = 40000; % number of iterations  
step = .5;    % average step size (var)

Npar = length(par);
j = 1;  % counter for accepted parameters
ind_best_err = 1; % index of best error

current_err = compute_error( par );

for i = 1:iter

  par2 = par + step * randn(1,Npar); % choose new parameters

  if ( max(par2) < UP & min(par2) > LW )

    new_err = compute_error( par2 );
    err_diff = ( current_err - new_err ) / current_err;
    exp_err_diff = exp( err_diff );

    if ( exp_err_diff > 1 | rand < exp_err_diff )
      par = par2;
      current_err = new_err;

      err(j) = current_err; % record stuff
      param1(j) = par(1);
      param2(j) = par(2);
      param3(j) = par(3);

      if ( current_err < err(ind_best_err) ) % record index of best error
	ind_best_err = j;
	best_error = current_err
      end
      
      j = j + 1;

    end
  end
end

subplot(2,2,1); plot(param1); title('par 1 (6)');
subplot(2,2,2); plot(param2); title('par 2 (sqrt 2)');
subplot(2,2,3); plot(param3); title('par 3 (pi)');
subplot(2,2,4); plot(err); title('L2 error');

ind_best_err
r = [ param1(ind_best_err) param2(ind_best_err) param3(ind_best_err)]  
err(ind_best_err)