    A PROBABILITY ESTIMATION PROBLEM WITH BINARY DATA

    As a first illustration of the mixture model and Dirichlet diffusion
    tree software, I generated a dataset in which each case was composed
    of ten binary attributes.  The distribution of these binary vectors
    was a mixture of four component distributions, in each of which the
    ten attributes were independent.  The four components each had
    probability 0.25 of being chosen to generate a case.  The
    probabilities for each component for each binary attributes were as
    shown below:

              1    2    3    4    5    6    7    8    9   10

         1   0.1  0.2  0.2  0.2  0.2  0.8  0.8  0.8  0.8  0.7
         2   0.1  0.8  0.8  0.8  0.8  0.2  0.2  0.2  0.2  0.7
         3   0.9  0.2  0.2  0.8  0.8  0.2  0.2  0.8  0.8  0.7
         4   0.9  0.8  0.8  0.2  0.2  0.8  0.8  0.2  0.2  0.7

    Each row gives the probabilities of each of the attributes being '1'
    for one component of the mixture.  The columns are for the ten binary
    attributes in each case.  The vectors generated in this way can be
    seen as coming from one of four "patterns": 0000011111, 0111100001,
    1001100111, and 1110011001, but with each bit of the chosen pattern
    having a small probability of being switched (ranging from 0.1 to 0.3)
    in any particular case.

    I generated 1000 cases from this distribution, which are stored in the
    file 'bdata'.  The first 500 are used for training, the rest for
    testing.

    A finite mixture model for the binary data.

    We will first see what happens when we model this data as a mixture of
    four components, which is the true number of components for this
    problem.  Each component of the mixture will give each target
    attribute a certain probability of being '1' rather than '0'.  These
    probabilities are determined from the "offset" parameters for each
    target, by means of the "logistic" function:

       Pr(target=1)  =  1 / (1 + exp(-offset))

    The offset parameters for each attribute for each component are given
    Gaussian prior distributions, with means and standard deviations that
    are hyperparameters.  These hyperparameters could be fixed, but here
    we will make them variable, separately for each target attribute, but
    with each of them linked to a top-level hyperparameter, common to all
    targets.

    We set up this mixture model with the 'mix-spec' command, as follows:

        > mix-spec blog.4 0 10 4 / x1 0.05:0.5:0.2 10

    Here, "blog.4" is the name of the log file used for this run, which is
    created by the 'mix-spec' command.  The arguments of 'mix-spec' that
    follow the name of the log file are the number of input attributes
    (always 0 for mixture models at present), the number of target
    attributes (10 for this problem), and the number of mixture components
    (set to the true value of 4 in the command above).  The specifications
    for priors are given following these arguments, after a "/" separator.

    The first prior specification gives the concentration parameter of the
    Dirichlet distribution for the mixing proportions for the components.
    If this specification starts with "x", this value is automatically
    divided by the number of components, which is the scaling required for
    the model to reach a sensible limit as the number of components goes
    to infinity.  The specification of "x1" above mildly favours unequal
    mixing proportions (recall that the true proportions are equal).

    The second argument after the "/" specifies the prior for the standard
    deviations of the "offsets" for components (which determine the
    probabilities for the binary target attributes).  In the hierarchical
    scheme specified here (described in general in prior.doc), there is a
    top-level standard deviation hyperparameter, given a prior in which
    the corresponding "precision" (standard deviation to the power -2) has
    a gamma distribution with mean 0.05 and shape parameter 0.5.  There is
    a lower-level standard deviation hyperparameter for each target; for
    each of these, the corresponding precision has a gamma distribution
    with mean given by the higher-level precision, and shape parameter of
    0.2.  Both these priors are rather vague (but the shape parameter of
    0.2 is vaguer than 0.5), so these hyperparameters can adapt to the
    structure of the data.

    The last argument of 'mix-spec' is the standard deviation for the
    means of the offsets for each of the targets, here set to 10.
    Currently, this standard deviation for the mean offset is the same for
    all targets, and is fixed in the 'mix-spec' command (it cannot be a
    variable hyperparameter).

    We can look at these specifications by calling 'mix-spec' with just
    the log file as an argument:

        > mix-spec blog.4

        Number of inputs:      0
        Number of targets:     10

        Number of components:  4

        Dirichlet concentration parameter: x1.000
        Prior for SD hyperparameters:       0.050:0.50:0.20
        Prior for mean component offsets:   10.000

    After the 'mix-spec' command, we need to specify that the targets are
    binary using the 'model-spec' command, as follows:

        > model-spec blog.4 binary

    We can then specify where the data comes from using 'data-spec':

        > data-spec blog.4 0 10 2 / bdata@1:500 . bdata@501:1000 .

    The number of inputs is here specified to be 0 (as it must be at
    present for mixture models), and the number of targets is specified to
    be 10.  These must match the values given to 'mix-spec'.  We also
    specify that the targets must be the integers 0 or 1 putting in a
    following argument of 2 (meaning that there are 2 possible values,
    which start at 0).  After a slash, we say where the inputs and targets
    come from.  Note that we have to say where the inputs come from even
    though there are 0 of them, but fortunately, we can then say that the
    targets come from the same place by just using the "." argument.  In
    the specification above, the data comes from the first 500 lines of
    the file 'bdata', with one case per line.  The remaining 500 lines are
    specified to be used for testing.  See data-spec.doc for more details.

    Finally, we need to specify how Markov chain sampling is to be done
    for this model.  At present, none of the standard Markov chain methods
    are allowed for mixture models, only the specialized procedures
    documented in mix-mc.doc.  Each of these procedures updates one of the
    three parts of the state - the values of the hyperparameters, the
    values of the parameters for the various mixture components, and the
    values of the indicator variables saying which mixture component is
    currently associated with each training case.  The values in each of
    these parts of the state can be updated by Gibbs sampling, using the
    corresponding procedure.  The following call of 'mc-spec' sets this
    up:

        > mc-spec blog.4 repeat 20 gibbs-indicators gibbs-params gibbs-hypers

    Here, the three Gibbs sampling operations are repeated 20 times each
    iteration, just to cut down on the number of states saved in the log
    file.

    We can now run the Markov chain simulation for 100 iterations with the
    following command:

        > mix-mc blog.4 100

    The simulation starts with a state in which all the training cases are
    associated with the same component, whose parameters are set to their
    means, as specified by hyperparameter values that are also set to
    their prior means.  The 'mix-gen' program could be used to initialize
    things differently (see mix-gen.doc).

    The above run takes 29 seconds on the system used (see Ex-system.doc).
    It can be run in the background, in which case progress can be
    monitored with the 'mix-display' and 'mix-plt' programs.  For example,
    after a while, you might see the following:

        > mix-display blog.4

        MIXTURE MODEL IN FILE "blog.4" WITH INDEX 50

        HYPERPARAMETERS

        Standard deviations for component offsets:

            0.326:     1.407    1.159    1.213    3.331    1.417
                       1.557    3.433    2.203    4.990    0.192

        Means for component offsets:

                      -0.370   -0.126   -0.538   +0.521   -0.620
                      -1.642   +1.510   -1.206   -1.502   +0.941

        PARAMETERS AND FREQUENCIES FOR COMPONENTS OF THE MIXTURE

           1: 0.286   -1.799   +0.931   +1.753   +1.481   +0.853
                      -1.503   -1.433   -0.942   -1.245   +1.165

           2: 0.252   +1.905   +1.054   +1.749   -2.588   -1.594
                      +1.062   +1.625   -1.532   -1.227   +1.252

           3: 0.250   -2.109   -1.505   -2.220   -0.961   -1.402
                      +2.265   +1.635   +1.463   +2.015   +1.082

           4: 0.212   +1.743   -1.050   -1.470   +1.502   +1.553
                      -1.593   -1.376   +1.569   +1.025   +0.864

    This displays the state at the last iteration computed so far (here
    iteration 50).  The hyperparameter values are shown first, with
    top-level hyperparameters on the left, before the colon.  The
    lower-level hyperparameters follow, ordered by the target attribute
    that they pertain to.  Values for the component parameters follow,
    sorted by the fraction of training cases that they are currently
    associated with.  These fractions are the first numbers after 1:, 2:,
    etc.; note that the actual mixing probabilities are not explicitly
    represented, but are instead always integrated over.  Any components
    that are not associated with any training case are omitted from the
    output, and are not explicitly represented in the state.  The other
    numbers shown for each component are the "offset" parameters for each
    target attribute.

    In this example, the simulation appears to have learned the essence of
    the distribution by the iteration shown above.  Recall that a positive
    offset corresponds to the probability of a 1 being greater than 1/2, a
    negative offset to the probability of a 1 being less than 1/2.  The
    four components shown above can thus be seen to correspond to the four
    patterns used to generate the cases, as described earlier.  Note that
    the last of the offset standard deviation hyperparameters (for the
    last target attribute) is quite small.  This indicates that the model
    has "learned" that the probabilities for the last target are the same
    for all components, and hence can be modeled most efficiently at the
    hyperparameter level, by the mean for that offset, rather than
    separately for each component.

    The indicators of which components are associated with each training
    case can also be examined with 'mix-display', and of course iterations
    other than the last can be viewed.  See mix-display.doc for details.

    One can also look at the progress of the simulation using 'mix-plt'.
    For example, the following will produce a time-plot of the cumulative
    proportions of the training cases associated with the four components
    (as sorted by decreasing frequency):

        > mix-plt t C1C2C3C4 blog.4 | plot

    Here, 'plot' is some suitable plotting program.  One can also just let
    the output of 'mix-plt' go to standard output, and then examine the
    numbers manually.  Other quantities can also be plotted, as described
    in mix-quantities.doc.

    As well as examining quantities with 'mix-display' and 'mix-plt', one
    can also produce a sample of cases from the predictive distribution
    using 'mix-cases'.  Predictive probabilities for test cases can be
    found using 'mix-pred'.  The following command gives the average log
    probability for all the test cases, using the last 80 iterations:

        > mix-pred pa blog.4 21:

        Number of iterations used: 80

        Number of test cases: 500

        Average log probability of targets:    -6.097+-0.069

    For comparison, the average log probability for these test cases
    using the true mixture distribution is -6.034.

    A countably infinite mixture model for the binary data.

    Even though the true distribution for this example is a mixture of
    four components, good results can nevertheless be obtained using a
    mixture with a countably infinite number of components.  The prior for
    mixture proportions used with such a model is designed so that a few
    of the infinite number of components have substantial probability, so
    the model does not result in an "overfitted" solution, in which every
    training case is "explained" as coming from a different component.

    We specify a countably infinite mixture model by simply omitting the
    argument to 'mix-spec' that gives the number of components.  For this
    example, we change the 'mix-spec' command used above to the following:

        > mix-spec blog.inf 0 10 / x1 0.05:0.5:0.2 10

    The Dirichlet prior specification for mixing proportions of "x1" is
    the same as for the mixture model with four components.  The "x"
    specifies scaling downward with the number of components, which
    produces a sensible limit as the number of components goes to
    infinity, as here.  The "x" is therefore required for infinite
    mixtures.

    The other arguments of 'mix-spec' are as described for the finite
    mixture model.  The 'model-spec' and 'data-spec' commands used are
    also the same:

        > model-spec blog.inf binary
        > data-spec blog.inf 0 10 2 / bdata@1:500 . bdata@501:1000 .

    We must change the 'mc-spec' command, however, since it is impossible
    to do Gibbs sampling for the component indicators associated with
    training cases - since their number is infinite, we can't compute all
    the required conditional probabilities.  However, we can instead use a
    specialized Metropolis-Hastings update, in which a new component to go
    with a training case is proposed with probability determined by the
    frequencies with which components are associated with other training
    cases.  The proposal is accepted or rejected based on the resulting
    change in likelihood.  The process can then be repeated any desired
    number of times, to produce an approximation to Gibbs sampling (the
    "approximation" is only with respect to convergence speed, the answer
    is exact, asymptotically).  With this change, we can use the following
    'mc-spec' command for this model:

        > mc-spec blog.inf repeat 20 met-indicators 10 gibbs-params gibbs-hypers

    This is the same as for the finite mixture model, except that 10
    repetitions of the Metropolis-Hastings update for the indicators are
    specified using 'met-indicators 10', in place of 'gibbs-indicators'.

    As before, we can now run the simulation with a command such as:

        > mix-mc blog.inf 100

    We can examine the states with 'mix-display'.  For example, once the
    'mix-mc' command completes, which takes 63 seconds on the system used
    (see Ex-system.doc), the state should be something like the following:

        > mix-display blog.inf

        MIXTURE MODEL IN FILE "blog.inf" WITH INDEX 100

        HYPERPARAMETERS

        Standard deviations for component offsets:

            0.359:     2.747    4.294    1.280    3.530    1.353
                       2.461    5.661    2.928    1.506    0.080

        Means for component offsets:

                      +0.023   +0.452   +0.792   +0.507   -0.084
                      -0.634   -4.794   +4.560   -0.293   +1.124

        PARAMETERS AND FREQUENCIES FOR COMPONENTS OF THE MIXTURE

           1: 0.310   -1.565   +0.617   +1.467   +1.191   +0.468
                      -0.569   -1.192   -0.397   -0.749   +1.003

           2: 0.254   -1.374   -1.134   -2.009   -1.008   -1.825
                      +2.384   +2.377   +1.579   +2.553   +1.080

           3: 0.238   +2.949   +0.979   +2.419   -2.694   -1.331
                      +1.597   +1.590   -1.689   -1.190   +1.124

           4: 0.168   +1.565   -1.687   -1.958   +1.762   +1.633
                      -1.365   -6.917   +1.001   +1.180   +0.919

           5: 0.012   +1.780   -3.038   -0.664   +5.257   -3.499
                      -8.315   -5.443   +0.816   -0.774   +1.143

           6: 0.012   +1.892   -0.991   +3.387   +0.513   +1.287
                      -4.129   +1.283   +8.346   +1.452   +0.882

           7: 0.004   +5.850   +1.412   +0.295   -5.867   +0.701
                      -0.403   -9.392   +7.748   -2.211   +1.163

           8: 0.002   +1.028   -3.102   +2.383   -1.070   +1.691
                      +0.497   +0.062   +0.505   -0.762   +1.226

    The output is similar to that seen for the mixture model with four
    components, except that eight components are associated with at least
    one training case at iteration 100, as shown above.  (The infinite
    number of remaining components are not explicitly represented, and are
    not displayed by 'mix-display'.)  Four of these components are
    associated with very few training cases, however (as seen from the
    fractions of the training set after the component number).  This is to
    be expected when the true distribution can be expressed using only
    four components.  The number of such low-frequency components will
    vary from iteration to iteration, as will other aspects of the state,
    in accordance with the variability in the posterior distribution.

    Since there are exactly four components in the real distribution, one
    would expect that the model with four components would perform better
    than a model with an infinite number of components.  Any penalty is
    quite small in this example, however, as can be seen by looking at the
    predictive probabilities:

        > mix-pred pa blog.inf 21:

        Number of iterations used: 80

        Number of test cases: 500

        Average log probability of targets:    -6.105+-0.069

    This is almost the same as the average log probability for the four
    component model.

    A number of other ways of sampling for the indicators for countably
    infinite mixture models are also available.  Here are two:

        > mc-spec blog.inf2 repeat 20 \
                    gibbs-ext-indicators 2 gibbs-params gibbs-hypers

        > mc-spec blog.inf3 repeat 20 \
                    met1-indicators gibbs1-indicators gibbs-params gibbs-hypers

    These methods may be more efficient than using the "met-indicators"
    operation for some models.  See mix-mc.doc for descriptions of these
    and other Markov chain operations for mixture models.

    A Dirichlet diffusion tree model for the binary data.

    Binary data can also be modeled using Dirichlet diffusion trees.  This
    is not the simplest example of how Dirichlet diffusion trees can be
    used, since the trees naturally produce real valued data.  To model
    binary data, these real values are treated as latent values, which are
    put through the logistic function to produce the probabilities for the
    binary data to be "1".  The binary data example here does not have the
    sort of hierarchical structure that would make use of a Dirichlet
    diffusion tree model desirable.  For both these reasons, the following
    example of how to model binary data with a Dirichlet diffusion tree
    may not be the best introduction to such models - you might want to
    read the real-valued example in Ex-mixdft-r.doc first.

    Nevertheless, this example does demonstrate how binary data can be
    modeled, and shows that even when there is no hierarchical structure
    to the data, a Dirichlet diffusion tree model can do quite well.

    We start by specifying a Dirichlet diffusion tree model, as follows:

        > dft-spec blog.dft 0 10 / 2 - 1

    As for a mixture mdoel specification, this says that there are 0
    "input" variables (as there must be at present), and 10 "target"
    variables.  The argument of "2" is the standard deviation for the
    diffusion process that produces the latent variables used to model the
    binary variables (this could have been a prior specification rather
    than a fixed value).  The last two arguments are the parameters of the
    divergence function used to determine (stochastically) when the
    diffusion paths for different cases diverge.  The general form of the
    divergence function is

       a(t) = c0 + c1/(1-t) + c2/(1-t)^2

    The parameters c0, c1, c2 are given at the end of the 'dft-spec'
    command, with "-" being equivalent to zero.  Any zero parameters at
    the end can just be omitted.  The specification above therefore
    corresponds to the divergence function a(t) = 1/(1-t).  Any of c0, c1,
    and c2 can be a prior specification, instead of a constant.  See
    dft-spec.doc for details, and Ex-mixdft-r.doc for an example.

    We next specify the data model, which is "binary", indicating binary
    data modeled using the logistic function:

        > model-spec blog.dft binary

    We go on to specify the source of the data, as for the mixture example
    above:

        > data-spec blog.dft 0 10 2 / bdata@1:500 . bdata@501:1000 .

    Finally, we need to specify how to do Markov chain sampling for the
    posterior distribution of the trees and latent values that underly the
    data.

        > mc-spec blog.dft repeat 10 gibbs-latent slice-positions met-terminals

    Since the hyperparameters are all fixed for this model, there is no
    need for Markov chain operations to update them.  The operations above
    will lead to the state of the Markov chain being only the latent
    values and the tree structure, with divergence times, but without the
    locations of non-terminal nodes (which will be integrated over).  The
    'gibbs-latent' operation does a Gibbs sampling scan over the latent
    values associated with the targets in training cases.  (Initially, the
    latent values are taken to be +1 for targets of "1", and -1 for
    targets of "0".)  The 'slice-positions' and 'met-terminals' operations
    both update the tree structure and divergence times.  See dft-mc.doc
    for details.

    We can run the Markov chain for 100 iterations with the following
    command:

        > dft-mc blog.dft 100

    We can see the progress of the Markov chain by examining various
    quantities, such as the latent values associated with training cases.
    The following command will plot the latent values for the first target
    for the first six training cases:

        > dft-plt t o1@1:6 blog.dft | plot

    The first four of these cases have "0" for their first target value,
    the last two have "1".  The plot should show that the latent values
    for the first four cases tend to be negative, whereas the latent
    values for the last two cases tend to be positive, at least once the
    Markov chain has approached convergence.

    The 'dft-mc' command above takes 11 minutes of time on the system used
    (see Ex-system.doc).  Once it is done, we can see how well the model
    predicts new data using the following command, which looks at the
    Dirichlet diffusion tree models produced by the last 60 iterations of
    the Markov chain:

        > dft-pred pa blog.dft 41:

        Number of iterations used: 60

        Number of test cases: 500

        Average log probability of targets:    -6.102+-0.067

    This takes 3.5 minutes on the system used (see Ex-system.doc).
    Performance is not much worse than for the four-component mixture
    model above, which is based on the true structure of the data.