comparison toolboxes/FullBNT-1.0.7/netlab3.3/hmc.m @ 0:e9a9cd732c1e tip

first hg version after svn

+function [samples, energies, diagn] = hmc(f, x, options, gradf, varargin)
+%HMC	Hybrid Monte Carlo sampling.
+%
+%	Description
+%	SAMPLES = HMC(F, X, OPTIONS, GRADF) uses a  hybrid Monte Carlo
+%	algorithm to sample from the distribution P ~ EXP(-F), where F is the
+%	first argument to HMC. The Markov chain starts at the point X, and
+%	the function GRADF is the gradient of the `energy' function F.
+%
+%	HMC(F, X, OPTIONS, GRADF, P1, P2, ...) allows additional arguments to
+%	be passed to F() and GRADF().
+%
+%	[SAMPLES, ENERGIES, DIAGN] = HMC(F, X, OPTIONS, GRADF) also returns a
+%	log of the energy values (i.e. negative log probabilities) for the
+%	samples in ENERGIES and DIAGN, a structure containing diagnostic
+%	information (position, momentum and acceptance threshold) for each
+%	step of the chain in DIAGN.POS, DIAGN.MOM and DIAGN.ACC respectively.
+%	All candidate states (including rejected ones) are stored in
+%	DIAGN.POS.
+%
+%	[SAMPLES, ENERGIES, DIAGN] = HMC(F, X, OPTIONS, GRADF) also returns
+%	the ENERGIES (i.e. negative log probabilities) corresponding to the
+%	samples.  The DIAGN structure contains three fields:
+%
+%	POS the position vectors of the dynamic process.
+%
+%	MOM the momentum vectors of the dynamic process.
+%
+%	ACC the acceptance thresholds.
+%
+%	S = HMC('STATE') returns a state structure that contains the state of
+%	the two random number generators RAND and RANDN and the momentum of
+%	the dynamic process.  These are contained in fields  randstate,
+%	randnstate and mom respectively.  The momentum state is only used for
+%	a persistent momentum update.
+%
+%	HMC('STATE', S) resets the state to S.  If S is an integer, then it
+%	is passed to RAND and RANDN and the momentum variable is randomised.
+%	If S is a structure returned by HMC('STATE') then it resets the
+%	generator to exactly the same state.
+%
+%	The optional parameters in the OPTIONS vector have the following
+%	interpretations.
+%
+%	OPTIONS(1) is set to 1 to display the energy values and rejection
+%	threshold at each step of the Markov chain. If the value is 2, then
+%	the position vectors at each step are also displayed.
+%
+%	OPTIONS(5) is set to 1 if momentum persistence is used; default 0,
+%	for complete replacement of momentum variables.
+%
+%	OPTIONS(7) defines the trajectory length (i.e. the number of leap-
+%	frog steps at each iteration).  Minimum value 1.
+%
+%	OPTIONS(9) is set to 1 to check the user defined gradient function.
+%
+%	OPTIONS(14) is the number of samples retained from the Markov chain;
+%	default 100.
+%
+%	OPTIONS(15) is the number of samples omitted from the start of the
+%	chain; default 0.
+%
+%	OPTIONS(17) defines the momentum used when a persistent update of
+%	(leap-frog) momentum is used.  This is bounded to the interval [0,
+%	1).
+%
+%	OPTIONS(18) is the step size used in leap-frogs; default 1/trajectory
+%	length.
+%
+%	See also
+%	METROP
+%
+
+%	Copyright (c) Ian T Nabney (1996-2001)
+
+% Global variable to store state of momentum variables: set by set_state
+% Used to initialise variable if set
+global HMC_MOM
+if nargin <= 2
+  if ~strcmp(f, 'state')
+    error('Unknown argument to hmc');
+  end
+  switch nargin
+    case 1
+      samples = get_state(f);
+      return;
+    case 2
+      set_state(f, x);
+      return;
+  end
+end
+
+display = options(1);
+if (round(options(5) == 1))
+  persistence = 1;
+  % Set alpha to lie in [0, 1)
+  alpha = max(0, options(17));
+  alpha = min(1, alpha);
+  salpha = sqrt(1-alpha*alpha);
+else
+  persistence = 0;
+end
+L = max(1, options(7)); % At least one step in leap-frogging
+if options(14) > 0
+  nsamples = options(14);
+else
+  nsamples = 100;	% Default
+end
+if options(15) >= 0
+  nomit = options(15);
+else
+  nomit = 0;
+end
+if options(18) > 0
+  step_size = options(18);	% Step size.
+else
+  step_size = 1/L;		% Default  
+end
+x = x(:)';		% Force x to be a row vector
+nparams = length(x);
+
+% Set up strings for evaluating potential function and its gradient.
+f = fcnchk(f, length(varargin));
+gradf = fcnchk(gradf, length(varargin));
+
+% Check the gradient evaluation.
+if (options(9))
+  % Check gradients
+  feval('gradchek', x, f, gradf, varargin{:});
+end
+
+samples = zeros(nsamples, nparams);	% Matrix of returned samples.
+if nargout >= 2
+  en_save = 1;
+  energies = zeros(nsamples, 1);
+else
+  en_save = 0;
+end
+if nargout >= 3
+  diagnostics = 1;
+  diagn_pos = zeros(nsamples, nparams);
+  diagn_mom = zeros(nsamples, nparams);
+  diagn_acc = zeros(nsamples, 1);
+else
+  diagnostics = 0;
+end
+
+n = - nomit + 1;
+Eold = feval(f, x, varargin{:});	% Evaluate starting energy.
+nreject = 0;
+if (~persistence | isempty(HMC_MOM))
+  p = randn(1, nparams);		% Initialise momenta at random
+else
+  p = HMC_MOM;				% Initialise momenta from stored state
+end
+lambda = 1;
+
+% Main loop.
+while n <= nsamples
+
+  xold = x;		    % Store starting position.
+  pold = p;		    % Store starting momenta
+  Hold = Eold + 0.5*(p*p'); % Recalculate Hamiltonian as momenta have changed
+
+  if ~persistence
+    % Choose a direction at random
+    if (rand < 0.5)
+      lambda = -1;
+    else
+      lambda = 1;
+    end
+  end
+  % Perturb step length.
+  epsilon = lambda*step_size*(1.0 + 0.1*randn(1));
+
+  % First half-step of leapfrog.
+  p = p - 0.5*epsilon*feval(gradf, x, varargin{:});
+  x = x + epsilon*p;
+  
+  % Full leapfrog steps.
+  for m = 1 : L - 1
+    p = p - epsilon*feval(gradf, x, varargin{:});
+    x = x + epsilon*p;
+  end
+
+  % Final half-step of leapfrog.
+  p = p - 0.5*epsilon*feval(gradf, x, varargin{:});
+
+  % Now apply Metropolis algorithm.
+  Enew = feval(f, x, varargin{:});	% Evaluate new energy.
+  p = -p;				% Negate momentum
+  Hnew = Enew + 0.5*p*p';		% Evaluate new Hamiltonian.
+  a = exp(Hold - Hnew);			% Acceptance threshold.
+  if (diagnostics & n > 0)
+    diagn_pos(n,:) = x;
+    diagn_mom(n,:) = p;
+    diagn_acc(n,:) = a;
+  end
+  if (display > 1)
+    fprintf(1, 'New position is\n');
+    disp(x);
+  end
+
+  if a > rand(1)			% Accept the new state.
+    Eold = Enew;			% Update energy
+    if (display > 0)
+      fprintf(1, 'Finished step %4d  Threshold: %g\n', n, a);
+    end
+  else					% Reject the new state.
+    if n > 0 
+      nreject = nreject + 1;
+    end
+    x = xold;				% Reset position 
+    p = pold;   			% Reset momenta
+    if (display > 0)
+      fprintf(1, '  Sample rejected %4d.  Threshold: %g\n', n, a);
+    end
+  end
+  if n > 0
+    samples(n,:) = x;			% Store sample.
+    if en_save 
+      energies(n) = Eold;		% Store energy.
+    end
+  end
+
+  % Set momenta for next iteration
+  if persistence
+    p = -p;
+    % Adjust momenta by a small random amount.
+    p = alpha.*p + salpha.*randn(1, nparams);
+  else
+    p = randn(1, nparams);	% Replace all momenta.
+  end
+
+  n = n + 1;
+end
+
+if (display > 0)
+  fprintf(1, '\nFraction of samples rejected:  %g\n', ...
+    nreject/(nsamples));
+end
+if diagnostics
+  diagn.pos = diagn_pos;
+  diagn.mom = diagn_mom;
+  diagn.acc = diagn_acc;
+end
+% Store final momentum value in global so that it can be retrieved later
+HMC_MOM = p;
+return
+
+% Return complete state of sampler (including momentum)
+function state = get_state(f)
+
+global HMC_MOM
+state.randstate = rand('state');
+state.randnstate = randn('state');
+state.mom = HMC_MOM;
+return
+
+% Set complete state of sampler (including momentum) or just set randn
+% and rand with integer argument.
+function set_state(f, x)
+
+global HMC_MOM
+if isnumeric(x)
+  rand('state', x);
+  randn('state', x);
+  HMC_MOM = [];
+else
+  if ~isstruct(x)
+    error('Second argument to hmc must be number or state structure');
+  end
+  if (~isfield(x, 'randstate') | ~isfield(x, 'randnstate') ...
+      | ~isfield(x, 'mom'))
+    error('Second argument to hmc must contain correct fields')
+  end
+  rand('state', x.randstate);
+  randn('state', x.randnstate);
+  HMC_MOM = x.mom;
+end
+return