wolffd@0: function [samples, energies, diagn] = metrop(f, x, options, gradf, varargin)
wolffd@0: %METROP	Markov Chain Monte Carlo sampling with Metropolis algorithm.
wolffd@0: %
wolffd@0: %	Description
wolffd@0: %	 SAMPLES = METROP(F, X, OPTIONS) uses the Metropolis algorithm to
wolffd@0: %	sample from the distribution P ~ EXP(-F), where F is the first
wolffd@0: %	argument to METROP.   The Markov chain starts at the point X and each
wolffd@0: %	candidate state is picked from a Gaussian proposal distribution and
wolffd@0: %	accepted or rejected according to the Metropolis criterion.
wolffd@0: %
wolffd@0: %	SAMPLES = METROP(F, X, OPTIONS, [], P1, P2, ...) allows additional
wolffd@0: %	arguments to be passed to F().  The fourth argument is ignored, but
wolffd@0: %	is included for compatibility with HMC and the optimisers.
wolffd@0: %
wolffd@0: %	[SAMPLES, ENERGIES, DIAGN] = METROP(F, X, OPTIONS) also returns a log
wolffd@0: %	of the energy values (i.e. negative log probabilities) for the
wolffd@0: %	samples in ENERGIES and DIAGN, a structure containing diagnostic
wolffd@0: %	information (position and acceptance threshold) for each step of the
wolffd@0: %	chain in DIAGN.POS and DIAGN.ACC respectively.  All candidate states
wolffd@0: %	(including rejected ones) are stored in DIAGN.POS.
wolffd@0: %
wolffd@0: %	S = METROP('STATE') returns a state structure that contains the state
wolffd@0: %	of the two random number generators RAND and RANDN. These are
wolffd@0: %	contained in fields randstate,  randnstate.
wolffd@0: %
wolffd@0: %	METROP('STATE', S) resets the state to S.  If S is an integer, then
wolffd@0: %	it is passed to RAND and RANDN. If S is a structure returned by
wolffd@0: %	METROP('STATE') then it resets the generator to exactly the same
wolffd@0: %	state.
wolffd@0: %
wolffd@0: %	The optional parameters in the OPTIONS vector have the following
wolffd@0: %	interpretations.
wolffd@0: %
wolffd@0: %	OPTIONS(1) is set to 1 to display the energy values and rejection
wolffd@0: %	threshold at each step of the Markov chain. If the value is 2, then
wolffd@0: %	the position vectors at each step are also displayed.
wolffd@0: %
wolffd@0: %	OPTIONS(14) is the number of samples retained from the Markov chain;
wolffd@0: %	default 100.
wolffd@0: %
wolffd@0: %	OPTIONS(15) is the number of samples omitted from the start of the
wolffd@0: %	chain; default 0.
wolffd@0: %
wolffd@0: %	OPTIONS(18) is the variance of the proposal distribution; default 1.
wolffd@0: %
wolffd@0: %	See also
wolffd@0: %	HMC
wolffd@0: %
wolffd@0: 
wolffd@0: %	Copyright (c) Ian T Nabney (1996-2001)
wolffd@0: 
wolffd@0: if nargin <= 2
wolffd@0:   if ~strcmp(f, 'state')
wolffd@0:     error('Unknown argument to metrop');
wolffd@0:   end
wolffd@0:   switch nargin
wolffd@0:     case 1
wolffd@0:       % Return state of sampler
wolffd@0:       samples = get_state(f);	% Function defined in this module
wolffd@0:       return;
wolffd@0:     case 2
wolffd@0:       % Set the state of the sampler
wolffd@0:       set_state(f, x);		% Function defined in this module
wolffd@0:       return;
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: if 0
wolffd@0: seed = 42;
wolffd@0: randn('state', seed);
wolffd@0: rand('state', seed)
wolffd@0: end
wolffd@0: 
wolffd@0: display = options(1);
wolffd@0: if options(14) > 0
wolffd@0:   nsamples = options(14);
wolffd@0: else
wolffd@0:   nsamples = 100;
wolffd@0: end
wolffd@0: if options(15) >= 0
wolffd@0:   nomit = options(15);
wolffd@0: else
wolffd@0:   nomit = 0;
wolffd@0: end
wolffd@0: if options(18) > 0.0
wolffd@0:   std_dev = sqrt(options(18));
wolffd@0: else
wolffd@0:   std_dev = 1.0;   % default
wolffd@0: end			
wolffd@0: nparams = length(x);
wolffd@0: 
wolffd@0: % Set up string for evaluating potential function.
wolffd@0: f = fcnchk(f, length(varargin));
wolffd@0: 
wolffd@0: samples = zeros(nsamples, nparams);		% Matrix of returned samples.
wolffd@0: if nargout >= 2
wolffd@0:   en_save = 1;
wolffd@0:   energies = zeros(nsamples, 1);
wolffd@0: else
wolffd@0:   en_save = 0;
wolffd@0: end
wolffd@0: if nargout >= 3
wolffd@0:   diagnostics = 1;
wolffd@0:   diagn_pos = zeros(nsamples, nparams);
wolffd@0:   diagn_acc = zeros(nsamples, 1);
wolffd@0: else
wolffd@0:   diagnostics = 0;
wolffd@0: end
wolffd@0: 
wolffd@0: % Main loop.
wolffd@0: n = - nomit + 1;
wolffd@0: Eold = feval(f, x, varargin{:});	% Evaluate starting energy.
wolffd@0: nreject = 0;				% Initialise count of rejected states.
wolffd@0: while n <= nsamples
wolffd@0: 
wolffd@0:   xold = x;
wolffd@0:   % Sample a new point from the proposal distribution
wolffd@0:   x = xold + randn(1, nparams)*std_dev;
wolffd@0:   %fprintf('netlab propose: xold = %5.3f,%5.3f, xnew = %5.3f,%5.3f\n',...
wolffd@0:   %	xold(1), xold(2), x(1), x(2));
wolffd@0: 
wolffd@0:   % Now apply Metropolis algorithm.
wolffd@0:   Enew = feval(f, x, varargin{:});	% Evaluate new energy.
wolffd@0:   a = exp(Eold - Enew);			% Acceptance threshold.
wolffd@0:   if (diagnostics & n > 0)
wolffd@0:     diagn_pos(n,:) = x;
wolffd@0:     diagn_acc(n,:) = a;
wolffd@0:   end
wolffd@0:   if (display > 1)
wolffd@0:     fprintf(1, 'New position is\n');
wolffd@0:     disp(x);
wolffd@0:   end
wolffd@0: 
wolffd@0:   r = rand(1);
wolffd@0:   %fprintf('netlab: n=%d, a=%f/%f=%5.3f (%5.3f), r=%5.3f\n',...
wolffd@0:   %	  n, exp(-Enew), exp(-Eold), a, exp(-Enew)/exp(-Eold), r);
wolffd@0:   if a > r	% Accept the new state.
wolffd@0:     Eold = Enew;
wolffd@0:     if (display > 0)
wolffd@0:       fprintf(1, 'Finished step %4d  Threshold: %g\n', n, a);
wolffd@0:     end
wolffd@0:   else			% Reject the new state
wolffd@0:     if n > 0
wolffd@0:       nreject = nreject + 1;
wolffd@0:     end
wolffd@0:     x = xold;	% Reset position 
wolffd@0:     if (display > 0)
wolffd@0:       fprintf(1, '  Sample rejected %4d.  Threshold: %g\n', n, a);
wolffd@0:     end
wolffd@0:   end
wolffd@0:   if n > 0
wolffd@0:     samples(n,:) = x;			% Store sample.
wolffd@0:     if en_save 
wolffd@0:       energies(n) = Eold;		% Store energy.
wolffd@0:     end
wolffd@0:   end
wolffd@0:   n = n + 1;
wolffd@0: end
wolffd@0: 
wolffd@0: if (display > 0)
wolffd@0:   fprintf(1, '\nFraction of samples rejected:  %g\n', ...
wolffd@0:           nreject/(nsamples));
wolffd@0: end
wolffd@0: 
wolffd@0: if diagnostics
wolffd@0:   diagn.pos = diagn_pos;
wolffd@0:   diagn.acc = diagn_acc;
wolffd@0: end
wolffd@0: 
wolffd@0: % Return complete state of the sampler.
wolffd@0: function state = get_state(f)
wolffd@0: 
wolffd@0: state.randstate = rand('state');
wolffd@0: state.randnstate = randn('state');
wolffd@0: return
wolffd@0: 
wolffd@0: % Set state of sampler, either from full state, or with an integer
wolffd@0: function set_state(f, x)
wolffd@0: 
wolffd@0: if isnumeric(x)
wolffd@0:   rand('state', x);
wolffd@0:   randn('state', x);
wolffd@0: else
wolffd@0:   if ~isstruct(x)
wolffd@0:     error('Second argument to metrop must be number or state structure');
wolffd@0:   end
wolffd@0:   if (~isfield(x, 'randstate') | ~isfield(x, 'randnstate'))
wolffd@0:     error('Second argument to metrop must contain correct fields')
wolffd@0:   end
wolffd@0:   rand('state', x.randstate);
wolffd@0:   randn('state', x.randnstate);
wolffd@0: end
wolffd@0: return