MC-SPEC: Specify how to do the Markov chain simulation.
MC-spec writes records to a log file that specify the operations
making up a single Markov chain iteration, and that also specify how
dynamic trajectories are to be computed. When invoked with just a log
file as argument, it displays the specifications stored at the given
index, or the last specifications stored if no index is given.
MC-spec is often invoked several times during a single simulation run.
The last specification in the log file is used for further iterations.
Usage:
mc-spec log-file { operation-spec } [ / trajectory-spec ]
or:
mc-spec log-file [ index | "all" ]
The log file must already exist. For the first form, specifications
are appended to the log file with index equal that of the last record
in the file, or 0 if the log file has no non-negative records.
For the second form, with no index specified, the last specifications
stored in the log file are displayed. If an index is specified in the
second form, the specifications in the record with the given index are
displayed. If the string "all" (no quotes) is given after the log
file, specifications at all indexes are displayed.
The allowed operations can be grouped in the categories general,
metropolis, dynamical, tempering, and slicing. Note that the slicing
operations should not be used when coupling chains, since they can
cause the random number sequence for different chains to get out of
synchronization.
Common argument types:
The following optional arguments are used for many of the operations
described below:
[ stepsize-adjust[:stepsize-alpha] ]
Specifies a stepsize to use, or a distribution of stepsizes
to randomly select from. A minus sign before a stepsize-adjust
value results in the adjustment being applied to uniform
stepsizes of one; otherwise, the adjustment is applied to the
application-specific stepsizes (which may be non-uniform). A
minus sign in front of stepsize alpha means that instead of the
usual Gamma distribution for stepsizes, a distribution over 'alpha'
orders of magnitude, uniform in the log domain, centred at the
value found using stepsize-adjust, is used. The default is
stepsize-adjust of 1 and stepsize-alpha of infinity (so that the
stepsize is a constant).
[ first[:last] ] ]
Specifies a range of coordinates to operate on. The order of
coordinates is specified by the application. Numbering starts
with 0. If omitted, the range defaults to all components. If
first is specified but last is omitted, last defaults to first.
This range can instead be specified by a single upper-case letter,
which is translated by the application into a range, which may
depend on the values of state variables not in the range. This
is done at present only by the 'src' module.
General operations:
repeat times
Repeat the following operations the specified number of times.
The repeated group of operations is terminated by 'end', or
the end of the list.
end
Terminates a group of operations.
plot
Prints the current values of those quantities specified in the
"xxx-mc" command on standard output, preceded by an index that
is reset to zero every iteration.
multiply-stepsizes factor [ first[:last] ] ]
Multiplies the stepsizes for coordinates in the given range by the
given factor. This modification persists for the remainder of the
operations in the current iteraton, unless modified again, or unless
the automatic stepsizes are overridden by a stepsize specification
starting with "-", or unless an application-specific operation is
done. (These last provisions are just to make implementation easier.)
If a multiply-stepsizes operation is inside a "repeat", it will
be done repeatedly. In particular, if factor is less than one,
the stepsizes will be successively smaller for each repetition.
slevel [ move [ random ] ]
Sets the manner in which the [0,1] value used for a Metropolis
accept/reject decision or for the setting the slice level is
obtained. With no arguments, the value is found by generating
a [0,1] value uniformly at random, ignoring the previous
value. If move is specified, the value is obtained from the
previous value by adding move, and also adding a uniform
random value in (-random,+random) if random is specified, and
then reflecting off the boundaries at 0 and 1 (implemented by
confining the value to [-1,+1] by wrapping circularly, and
using its absolute value as the [0,1] value). The value is
also updated appropriately when a move is accepted, to achieve
reversibility.
The default is to generate a [0,1] value randomly, ignoring
the previous value. Also, if no slevel operation is done, the
value will not be recorded in the log file, and hence will not
be available as the "s" quantity.
(any operation not otherwise defined here) [ number [ number ] ]
Invoke the application-specific update procedure, passing the
name of the operation, and the numerical parameters following
it, which default to zero.
These operations may or may not be allowed when coupling.
Metropolis/Gibbs operations:
The "-b" option for the Metropolis operations below causes the
"Barker" or "Boltzmann" acceptance probability to be used, rather than
the usual "Metropolis" form. If R is the ratio of probability
densities for the new and old states, the "-b" option acceptance
probability is R/(1+R), whereas the Metropolis form is min(1,R).
metropolis [ -b ] [ stepsize-adjust[:stepsize-alpha] ]
Do a simple Metropolis update, to all components at once, using
a Gaussian proposal distribution.
met [ -b ] [ stepsize-adjust[:stepsize-alpha] ]
Same as "metropolis".
met-1 [ -b ] [ -r ] [ stepsize-adjust[:stepsize-alpha] [ first[:last] ] ]
Do single-component Metropolis updates, for the range of components
given (default all). If no -r argument is given, the components are
all updated in order; with the -r option, just one randomly-chosen
component is updated.
rgrid-met [ -b ] [ stepsize-adjust[:stepsize-alpha] ]
Do a random-grid Metropolis update, to all components at once, using
a uniform proposal distribution centred on the current state and
extending to plus and minus the stepsize specified. This uniform
proposal is found by randomly positioning a grid of hypercubes and
then proposing a move to the centre of the hypercube containing the
current state. This leads to exact coalescence of coupled chains
whenever both chains' states are in the same hypercube.
rgrid-met-1 [ -b ] [ -r ] [ stepsize-adjust[:stepsize-alpha]
[ first[:last] ] ]
Do single-component random-grid Metropolis updates, for the range
of components given (default all). If no -r argument is given, the
components are all updated in order; with the -r option, just one
randomly-chosen component is updated.
binary-gibbs [ -r ] [ first[:last] ] ]
Do Gibbs sampling updates for the range of coordinates given
(default all), or do a Gibbs sampling update for one of these
coordinates chosen at random, if -r is specified. The Gibbs sampling
is done by ASSUMING that the range of the conditional distribution
is the values 0 and 1.
gaussian-gibbs [ -r ] [ first[:last] ] ]
Do Gibbs sampling updates for the range of coordinates given
(default all), or do a Gibbs sampling update for one of these
coordinates chosen at random, if -r is specified. The Gibbs sampling
is done by just ASSUMING that the conditional distribution is Gaussian,
with its mean and variance being found by evaluating the energy at
values of -1, 0, and +1 for the coordinate being updated. This
WILL NOT BE VALID if the conditional distribution is not Gaussian.
Dynamical operations:
heatbath [ decay ]
Do a heatbath update of the momentum variables. If decay is
zero (the default), the current momentum is forgotten, and
new values are picked randomly from their distribution. If
decay is non-zero, the momentum variables are multiplied by
decay, and Gaussian noise with variance 1-decay^2 is then
added.
radial-heatbath
Do a radial heatbath update of the momentum variables, in
which the squared magnitude of the momentum is sampled
from a chi-squared distribution, while the direction is
left unchanged.
mix-momentum factor
Does a random rotation of all momentum variables, by adding
mean zero Gaussian noise to each momentum variable, and then
multiplying all by the factor needed to produce the same
kinetic energy as before. The variance of the noise added is
the average squared magnitude of the current momentum variables
times the specified factor squared.
negate
Negate the momentum variables.
dynamic steps [ stepsize-adjust[:stepsize-alpha] ]
Follow a dynamical trajectory for the given number of steps,
always accepting the result (hence does not leave distribution
exactly invariant).
permuted-dynamic steps [ stepsize-adjust[:stepsize-alpha] ]
Like 'dynamic', but the order of approximations is randomly
permuted.
hybrid steps[:window[:jump]]
[ stepsize-adjust[:stepsize-alpha] [ first[:last] ] ]
OR
hybrid max-steps/max-ok[:jump]
[ stepsize-adjust[:stepsize-alpha] [ first[:last] ] ]
Use the results of following a dynamical trajectory as a
candidate state for a Metropolis update. There are two forms,
differing in the way length of a trajectory and its acceptance
are determined. In both forms, states along a trajectory are
looked at only at every 'jump' steps (and either 'steps' or
'max-steps' must be multiples of 'jump'). The default is a jump
of one. The stepsize is determined as for the 'dynamic' operation.
In both forms, the momentum is negated if the proposal is accepted
in order to make the step reversible.
In the first form, acceptance is based on a 'window' of states
at the beginning and end of the trajectory. In the second form,
acceptance is based on a single state, but the number of steps
in the trajectory is not fixed - instead, the trajectory ends
after max-ok states that would be accepted have been found (looking
only every 'jump' states), or when max-steps states of any sort have
been produced; if the trajectory ends for the first reason, acceptance
is guaranteed. The default if neither 'max-ok' nor 'window' is
specified is standard hybrid Monte Carlo (ie, a window of one).
The first:last option restricts the update to only the specified
range of coordinates (numbered from 0), same as if the others
had a stepsize of zero.
tempered-hybrid temper-factor steps[:window[:jump]]
[ stepsize-adjust[:stepsize-alpha] [ first[:last] ] ]
Like hybrid (first form), but the trajectory is "tempered" by
multiplying the momenta by temper-factor in the first half of
the portion of the trajectory outside any windows, and by
1/temper-factor in the second half.
spiral temper-factor steps [ stepsize-adjust[:stepsize-alpha] ]
Performs spiral dynamics, in which steps of Hamiltonian dynamics
alternate with multiplication/division of the momentum by the
specified tempering factor. The specified number of steps are
randomly divided into steps "before" and "after" the current
state. A state is the selected from all that are generated
according to their probabilities adjusted for expansion or
contraction by multiply/divide.
double-spiral temper-factor steps [ stepsize-adjust[:stepsize-alpha] ]
Like spiral, but with a reversal of multiply vs. divide at a
randomly chosen point, producing a double spiral.
Slicing operations:
slice-1 [ -r ] [ -s [-]factor[/threshold] ]
[ stepsize-adjust[:stepsize-alpha]
[ max-steps [ first[:last] ] ] ]
Do single-component slice sampling, for the range of components
given (default all). If no -r argument is given, the components are
all updated in order; with the -r option, just one randomly-chosen
component is updated.
If the max-steps argument is positive, it gives the maximum number
of intervals to create when using the stepping out procedure to
find an interval around the current point. A value of zero
indicates that stepping out should be done with no limit (this
is the default). If max-steps is negative, the interval is found
by doubling, with minus max-steps again being the limit on the
number of intervals (so setting max-steps to -1 gives the same
results as setting it to 1). (There is no way to double without
limit, but setting max-steps to -1000 is the same for all practical
purposes.)
With no -s option, when a point is rejected, one end of the interval
shrinks to the rejected point. If the -s option is given, with no
"-" before "factor", the interval may then be shrunk further by
dividing it up into "factor" subintervals, and retaining only the
subinterval that contains the current point. This additional
shrinkage is done if the energy is more than "threshold" above the
slice level (the default is a threshold of zero). If "-" is given
before "factor", the initial shrinkage to the rejected point is not
done - only the shrinkage by "factor" is done. As a consequence,
shrinkage of the interval can be disabled entirely using "-s -1".
slice [ -g | -G ] [ stepsize-adjust[:stepsize-alpha] ]
Do multi-dimensional slice sampling with the given stepsize. At
present, the initial hyperrectangle found cannot be expanded (either
by stepping out or by doubling). After a point is randomly chosen
from the hyperrectangle, it is shrunk in certain directions if the
chosen point is not in the slice, and another random point chosen,
until eventually a point in the slice is found. If neither the -g
nor the -G option is specified, all coordinate directions are shrunk.
If -g is specified, only one coordinate direction is shrunk - that
for which the product of the absolute value of the energy gradient
and the current dimension of the hyperrectange in that coordinate
is the maximum. If -G is specified, this coordinate direction
is shrunk (as for -g), along with all other coordinate directions
for which the corresponding product is at least half this maximum.
slice-gaussian [ -e ] [ stepsize-adjust[:stepsize-alpha] ]
Do multi-dimensional slice sampling with Gaussian-distributed
"crumbs". The distribution of the first crumb has diagonal
covariance, with standard deviations given by the stepsizes.
If no option is given, this is also the distribution of later
crumbs. When -e is specified, later crumbs have their standard
deviations rescaled based on the energy of the most recent trial
point, so as to (hopefully) lead to a good chance of the next
trial point lying in the slice.
slice-over [ -r ] [ refinements [ refresh-prob
[ stepsize-adjust[:stepsize-alpha] [ max-steps [ first[:last] ] ] ] ] ]
Do overrelaxed slice sampling, for the range of components given,
or for one such component chosen randomly, if the -r option is given.
The endpoints are computed using the given number of refinements
(default zero). The refresh probability is the probability of doing
an ordinary slice sampling update rather than an overrelaxed one;
it defaults to zero. The meaning of max-steps is as for slice-1.
slice-inside steps [ stepsize-adjust[:stepsize-alpha] ]
Performs multivariate slice sampling by reflection from inside
points, proceeding for the indicated number of steps, with the
indicated stepsize adjustment factor. The momentum is negated
in such a way as to make the operation reversible.
slice-outside steps[/in-steps] [ stepsize-adjust[:stepsize-alpha] ]
Performs multivariate slice sampling by reflection from outside
points, proceeding for the indicated number of steps, with the
indicated stepsize adjustment factor. The momentum is negated
at the end in such a way as to make the operation reversible.
The in-steps argument gives a limit on the number of steps that
can be inside the slice; it defaults to the same as steps (thereby
having no effect). Setting in-steps to less than steps may decrease
the chances of the trajectory ending on an outside point, and being
rejected.
Tempering/annealing operations:
See also the descriptions under "dynamical operations" of the spiral,
double-spiral, and tempered-hybrid operations.
sim-temp
Do a metropolis update of the simulated tempering inverse
temperature, with a proposal of changing the temperature
index in accord with the current tempering direction. The
direction is negated if the proposal is accepted to make the
step reversible.
rand-dir
Randomize the tempering direction.
neg-dir
Negate the tempering direction.
temp-trans
Perform a tempered transition, with components given by the
following operations (terminated by 'end'). The components
are done in forward order for the first half of the tempered
transition, in reverse order for the second. For this to
work correctly, the components must be reversible in themselves.
AIS
Proceed to the next step of an annealed importance sampling
run, with the tempering index increased by one, adjusting the
importance weight appropriately. If the tempering index is
already at the final value (inverse temperature of one), a new
state is randomly generated from the distribution at inverse
temperature zero, and the tempering index is set to the beginning
of the tempering schedule.
Test operations:
The following operations are useful for testing and research only,
or possibly when trying to manually get your Markov chain started
in a reasonable place. They do not leave the desired distribution
invariant, and hence should not be used as part of any final sampling
scheme.
set-temp index
Sets the temperature for simulated tempering to the given
index from the tempering schedule (see mc-temp-sched.doc),
with 0 being the highest temperature. The direction for
temperature change is also set to +1 (to lower temperatures).
multiply-momentum factor
All momentum variables are multiplied by the given factor.
set-momentum value
All momentum variables are set to the given value.
set-value value [ first[:last] ]
Sets the "position" variable from first to last to the specified
value. All variables are set if first:last is omitted, just
first is set if last is omitted. Variable indexes start at 0.
Certain of these operations are normally used in standard combinations,
in particular the following:
heatbath hybrid <steps> Standard hybrid Monte Carlo
heatbath <decay> hybrid <steps> negate Persistent Hybrid Monte Carlo
rand-dir sim-temp Standard simulated tempering
sim-temp neg-dir Persistent Simulated tempering
All operations are reversible (other than 'ais', 'repeat', and 'end'
for which the concept is not applicable), except for 'dynamic',
'permuted-dynamic', 'multiply-momentum', 'set-momentum', and perhaps
the application-specific operations. However, note that in general
sequential combinations of reversible operations are not reversible.
Defaults are all components, decay of zero, stepsize-adjust of one,
stepsize-alpha of infinity, window of one, and jump of one.
The trajectory specification can have one of the following forms:
leapfrog [ halfp | halfq ] [ N-approx | frac ] [ approx-file ]
Use the leapfrog method. The "halfp" or "halfq" option specifies
whether the half-steps at the beginning and end of the trajectory
are for p (momentum) or q (position). The default is "halfp".
The N-approx option results in N-approx leapfrog steps for each
overall step, which use successive approximations to the energy
(whose average must be the true value), in an order randomly selected
at the start of the trajectory. The default for N-approx is one.
The optional approx-file arguement gives the name of a file containing
further information on how the approximations are constructed, whose
contents are specific to the application.
opt2 [ rev | sym ] [ firstq | firstp ] [ frac ]
Apply McLachlan and Atela's "optimal" two stage method, with the
first step being done for the position (firstq) or momentum
(firstp) variables. The default is "firstp". This method is
not symmetrical. The "rev" option gives the reversed version.
The "sym" option gives the ordinary followed by the reversed
version, with half the stepsize for each, thereby producing
a symmetrical method. Only the symmetrical method should be
used in conjuntion with "hybrid" and "tempered hybrid" operations.
gen2 [ rev | sym ] [ firstq | firstp ] a2 [ frac ]
Apply a two stage second-order method with parameter a2, as
described by McLachlan and Atela. Firstq and firstp are as for
opt2. This is probably only interesting for research, since the
most useful values for a2 are covered by "leapfrog" and "opt2".
The "rev" and "sym" options are as for opt2.
opt4 [ rev | sym ] [ frac ]
Apply McLachlan and Atela's "optimal" four stage method for
quadratic kinetic energy. The "rev" and "sym" options are as for opt2.
The frac option for all these specifications causes the gradient of
the energy use to be based on that fraction of some set of terms (eg,
training cases). The default is to use all terms (ie, frac of 1).
This option may be ignored by some applications.
The default if no specification is given is "leapfrog halfp". The
opt2, gen2, and opt4 methods are taken from a paper by Robert I.
McLachlan and Pau Atela, "The accuracy of symplectic integrators",
Nonlinearity, vol. 5, pp. 541-562, 1992. (Note: there is a typo in
their description of opt2 in table 2. b1=1-b2 should be 1-1/sqrt(2).)
The list of operations is stored in a record of type 'o'; the
trajectory specification in a record of type 't'.
Copyright (c) 1995-2019 by Radford M. Neal