DFT-MC: Do Markov chain sampling for a Dirichlet diffusion tree model.
The dft-mc program is the specialization of xxx-mc to the task of
sampling from the posterior distribution of the overall model
hyperparameters, the diffusion tree parameters, the latent vectors for
training cases, and the structures, divergence times, and node
locations of the diffusion trees. If no training data is specified,
the prior distribution for the overall hyperparameters and diffusion
tree parameters is sampled from, in which case only the gibbs-sigmas
and slice-div operations below do anything (this is useful only for
testing).
The state of the simulation has up to five parts:
1) overall hyperparameter values
2) parameters associated with each diffusion tree
3) diffusion tree structure and divergence times
4) latent vectors for training cases
5) locations for non-terminal nodes in the diffusion trees
The first three parts are always present, once sampling has begun, but
latent vectors and node locations are sometimes optional - if not
present, they are integrated over. Four situations can be
distinguished in this respect, first according to whether the model is
Gaussian (or there is no model - ie, no noise) or non-Gaussian (binary
data or real data with t-distributed noise), and second according to
whether there is only one diffusion tree or there is more than one
diffusion tree. The following table shows for each situation whether
latent vectors and/or node locations are absolutely required (req),
and whether they are needed in practice to do useful sampling (need):
GAUSSIAN MODEL NON-GAUSSIAN MODEL
ONLY ONE TREE req: neither req: latent vectors or
node locations
need: neither need: latent vectors
SEVERAL TREES req: node locations req: node locations
need: node locations need: node locations and
latent vectors
Even when keeping latent vectors and/or node locations around is not
necessary, it is possible that doing so will reduce computation time.
Occasionally, one might wish to discard latent vectors and/or node
locations that are useful for sampling but not absolutely required at
the end of each iteration, in order to conserved on disk space,
resampling them at the beginning of the next iteration.
Note that the locations of terminal nodes are never stored - if there
is only one tree, they are the same as the latent vectors; if there is
more than one tree, these terminal node locations can easily be
integrated over given the latent vectors.
If the overall hyperparameters or tree parameters are not present when
dft-mc is run, they are set to the location parameters of their
priors. If trees are not present to begin, random trees are generated
by sucessively pairing randomly selected terminal or non-terminal
nodes, with divergence times ranging from zero to 0.1. As described
below, latent vectors and node locations are also created as required
by some operations, whereas for other operations, one or the other
must be present beforehand.
The generic features of this program are described in xxx-mc.doc.
Note, however, that the dynamical state for the standard Markov chain
operations is null for this module, so the standard Markov chain
operations (described in mc-spec.doc) are of no use. The state can be
updated only by using the application-specific operations described
below. Time requirements are given in terms of the number of training
cases, N, the number of trees, T, and the number of target variables,
V, and the average depth of a node in the tree, D, which is typically
of order log(N). It is assumed that T is much less than N. The number
of scans, K, that is an argument to some operations defaults to one.
Typical uses of these operations are described after the list.
gibbs-hypers [ K ]
Does K Gibbs sampling scans updating the hyperparameters
controlling the standard deviations of the diffusions for each
variable, for each tree. If node locations and/or latent values
are present, the updates are done conditional on these values;
otherwise values for node locations are generated, as for
sample-locations and then discarded once the hyperparameters
have been updated. Locations for terminal nodes are also
generated before the K Gibbs sampling updates, and discarded
afterwards. The time required for this operation is O(N*T*K).
Latent vectors and node locations must not both be absent when
the data model is non-Gaussian. Node locations must not be
absent when there is more than one tree.
gibbs-noise [ K ]
Does K Gibbs sampling scans updating the noise standard deviations,
if the data model is real; does nothing if the model isn't real.
If latent vectors are present, the updates are done conditional
on these values; otherwise values for latent vectors are
generated, as for sample-latent, and then discarded once the noise
parameters have been updated. Locations for terminal nodes are
also generated before the K Gibbs sampling updates, and discarded
afterwards, as are case-by-case noise variances, if there is a third
level in the prior for the noise. The time required for this
operation is O(N*T*K).
Latent vectors and node locations must not both be absent when the
data model is non-Gaussian. Node locations must not be absent when
there is more than one tree.
gibbs-sigmas [ K ]
Does K Gibbs sampling scans for both diffusion standard deviations
and noise standard deviations, combining the effects of gibbs-hyper
and gibbs-noise. Equivalent to "repeat K gibbs-hypers gibbs-noise".
slice-div [ scale [ K ] ]
Does K slice sampling updates for the parameters of the divergence
functions for each of the diffusion trees. The single variable slice
sampling procedure is used, applied to the logarithms (base e) of the
parameters, using an initial interval of width "scale" (default one),
which is not expanded. This takes O(N*T*K) time.
sample-latent
Samples values for the latent vectors underlying training cases,
independently of any previous such values. If non-terminal node
locations are present, sampling is done conditional on these values,
which is possible for any data model. If node locations are not
present, this operation is possible only when there is only one tree,
and the targets are real-valued, with no data model or with a Gaussian
noise model. The time required is O(N*T).
sample-locations
Samples values for the locations of non-terminal nodes. This is
possible only when there is only one tree. If latent vectors are
present, the sampling is conditional on their values. If latent
vectors are not present, sampling is conditional on the observed
data, which is possible only when the data is real, with no data model
or with a Gaussian noise model. The time required is O(N).
gibbs-latent [ K ]
Does K Gibbs sampling scans for the latent vectors underlying
training cases, updating the vector for each training case in turn.
This operation is possible for any data model. If node locations
are present, the Gibbs updates are conditional on these locations,
and this operation is equivalent to sample-latent. If node locations
are absent, udpates are based on integrating over the node locations,
which is possible only when there is only one tree. The time required
is O(K*N*T).
If latent vectors were not present to begin, they are created,
as if create_latent were done.
gibbs-locations [ K ]
Does K Gibbs sampling scans for the node locations of each tree.
In each update, the entire set of node locations for one tree is
updated at once, conditional on the latent vectors or the data.
If there is only one tree, this operation is equivalent to
sample-locations. This operation is possible for any data model
if latent vectors are present; if latent vectors are absent, it is
possible only when the data is real, with no data model or with a
Gaussian noise model. The time required is O(K*N*T^2).
If node locations were not present to begin, they are created,
as if create_locations were done.
mh-latent [ scale [ K ] ]
This operation does nothing if node locations are present, or if
there is no data model, and is not allowed if there is more than
one tree. Otherwise, it does K Metropolis-Hastings updates for the
latent vectors underlying training cases. Each update proposes a
change to the entire set of latent vectors simultaneously, with
the proposed vectors found by multiplying the current vectors by
sqrt(1-scale^2) and then adding Gaussian noise with mean zero and
covariance given by the covariance of the latent vectors times scale^2.
The default for scale is one, in which case the current vectors do not
affect the proposal, but this is probably not desirable in most cases.
If latent vectors were not present to begin, they are created,
as if create_latent were done.
discard-latent
Discard the values stored for the latent vectors for each training
case. Does nothing if there are no latent vectors.
Note that if latent vectors are discarded during a run, they
should be re-created (if desired) with sample-latent, not with
create-latent, as the later does not preserve the desired distribution.
discard-locations
Discard the values stored for node locations. Does nothing if
node locations aren't present.
Note that if node locations are discarded during a run, they
should be re-created (if desired) with sample-locations, not with
create-locations, as the later does not preserve the desired
distribution.
create-latent
Creates latent vectors if they do not already exist. The new
vectors be equal to the data vectors if the data is real, and
will have values of plus or minus one if the data is binary
(with +1 for 1 and -1 for 0).
Note that this operation should be used only in order to set
up a run properly - discarding and recreating the latent variables
during the course of a run will not preserve the posterior
distribution.
create-locations
Creates node locations if they do not already exist. The new
location vectors will all be zero.
Note that this operation should be used only in order to set
up a run properly - discarding and recreating the node locations
during the course of a run will not preserve the posterior
distribution.
slice-positions [ [-]K ]
Does K passes over the terminal nodes of each tree, performing for
each terminal node an update of the following sort. One of the
direct ancestors of the terminal node is selected uniformly at
random. This ancestral node is moved to a new position along
the path from the root to the terminal node, and a new location
for this ancestral node is drawn, if node locations are being kept.
This new position for the ancestral node is identified by the time of
divergence at this node of the path to the first terminal from the path
to the second terminal. This divergence time is updated using slice
sampling, with the initial interval being from zero to the divergence
time of the ancestor's other child (not on the path to the selected
terminal node). The choice of a new time is made on the basis of the
probability of the new configuration of the tree, with the location of
the ancestral node integrated over.
Latent vectors must be present for this operation if the data model
is non-Gaussian. Node locations must be present if there is more
than one tree.
If the average number of points tested before one is accepted is E,
the time required for this operation is O(K*N*T*D*E).
If "-" precedes K, incremental recomputation of likelihoods (when
node locations are absent) is disabled. This may be useful as a way
of reducing memory usage, at a large cost in increased time, but it
is mainly intended for testing.
slice-positions2 [ [-]K ]
Like slice-positions, except for the method used to choose which
ancestor of the terminal node is to be updated. With slice-positions2,
a different terminal node is chosen at random, and the node updated
is the most recent common ancestor of this terminal node and the
original terminal node.
met-terminals [ [-]K ]
Does a series of Metropolis-Hastings updates for each terminal node,
x, in turn, for all the trees, repeating this K times. Each
update proposes to move x's parent to a segment of the tree and a
time within that segment that are chosen by simulating the generation
of a hypothetical data point from the tree with x omitted. The point
where this generation process diverges defines the new position of x's
parent; the current segment is a possible position, with a new time.
If the move is accepted, the other child of x's present parent becomes
a child of it's grandparent, and the parent of x will have the new
segment's initial node as its parent and the new segment's final node
as its other child.
The change is accepted or rejected based on the change in probability
excluding the probability associated with the generation of the path
in the tree to x, since this is already accounted for by the
probability of proposing a change. If latent vectors and/or node
locations are present, the joint probability including these values
is used, except that the location of the x's parent is integrated over.
A new location for the parent of x is generated afterward (regardless
of whether the move was accepted or rejected). If latent vectors or
node locations are absent, the marginal probability integrating over
these values is used. Note that (as always) the terminal node
locations are integrated over, regardless of whether latent vectors
or non-terminal node locations are present.
Latent vectors must be present for this operation if the data model
is non-Gaussian. Node locations must be present if there is more
than one tree.
The time required for this operation is O(K*N*T*D).
If "-" precedes K, incremental recomputation of likelihoods (when
node locations are absent) is disabled. This may be useful as a way
of reducing memory usage, at a large cost in increased time, but it
is mainly intended for testing.
met-terminals-uniform [ [-]K ]
Like met-terminals, except that the new position of x's parent
is found by first choosing a segment of the tree (without x) uniformly
at random, and then choosing a time within this segment uniformly.
The new position is accepted or rejected based on the change in
probability, accounting for the asymmetry of the proposal probabilities
when the current and proposed segments are of different lengths.
met-nonterminals [ [-]K ]
Does a series of Metropolis-Hastings updates for each nonterminal node,
x, except the root, for all the trees, repeating this K times. Each
update proposes to move x's parent to a segment of the tree and a
time within that segment that are chosen by simulating the generation
of a hypothetical data point from the tree with x omitted. The point
where this generation process diverges defines the new position of x's
parent; the current segment is a possible position, with a new time.
If this new divergence time for x's is later than the divergence time
for x itself, the process is repeated, until a valid divergence time
is generated, or until 100 attempts have failed.
If the move is accepted, the other child of x's present parent becomes
a child of it's grandparent, and the parent of x will have the new
segment's initial node as its parent and the new segment's final node
as its other child.
The change is accepted or rejected based on the change in probability
excluding the probability associated with the generation of the path
in the tree to x, since this is already accounted for by the
probability of proposing a change. If latent vectors and/or node
locations are present, the joint probability including these values
is used, except that the location of the x's parent is integrated over.
A new location for the parent of x is generated afterward (regardless
of whether the move was accepted or rejected). If latent vectors or
node locations are absent, the marginal probability integrating over
these values is used. Note that (as always) the terminal node
locations are integrated over, regardless of whether latent vectors
or non-terminal node locations are present.
Latent vectors must be present for this operation if the data model
is non-Gaussian. Node locations must be present if there is more
than one tree.
The time required for this operation is O(K*N*T*D).
If "-" precedes K, incremental recomputation of likelihoods (when
node locations are absent) is disabled. This may be useful as a way
of reducing memory usage, at a large cost in increased time, but it
is mainly intended for testing.
The 'e' quantity described in mc-quantities.doc is updated by
slice-params. The 'r', 'm', and 'D' quantities are updated by the
Metropolis operations, with 'm' and 'D' being set based on the update
for one training case, selected at random each time.
Here are some typical sequences of operations for each of the four
situations described earlier (others are possible too):
ONLY ONE TREE, GAUSSIAN MODEL:
slice-positions gibbs-sigmas
or: sample-locations slice-positions gibbs-sigmas
ONLY ONE TREE, NON-GAUSSIAN MODEL:
gibbs-latent slice-positions gibbs-sigmas
or: gibbs-latent sample-locations slice-positions gibbs-sigmas
or: mh-latent 0.1 slice-positions gibbs-sigmas
SEVERAL TREES, GAUSSIAN MODEL:
gibbs-locations slice-positions gibbs-sigmas
SEVERAL TREES, NON-GAUSSIAN MODEL:
create-latent gibbs-locations sample-latent slice-positions gibbs-sigmas
The gibbs-sigmas operation is not needed if there are no variable
hyperparameters. A slice-div operation should be added after
gibbs-sigmas if the parameters of the divergence functions are not
fixed. Adding a met-terminals operation after the slice-positions
operation might speed up convergence. If the cumulative divergence
function doesn't tend to infinity at time one a met-terminals
operation is definitely needed, and an mh-latent operation should be
used for non-Gaussian models. It may be advantageous to adjust the
number of scans for various operations, or to repeat certain sections
several times, to achieve the right balance of effort.
Annealed importance sampling (AIS) is implemented (but not very well
tested). The annealed distributions are based on raising the
likelihood of the data given latent values to the power of the inverse
temperature. Hence, annealing can be used only when there is a data
model, and only when latent values are represented (at the point where
importance weights have to be computed - ie, when the "AIS" operation
is done, except when it starts a new annealing sequence).
Tempered transitions are not implemented.
The current Markov chain operations are not suitable for use with
circular coupling.
Copyright (c) 1995-2003 by Radford M. Neal