    SAMPLING FROM A UNIVARIATE NORMAL DISTRIBUTION

    I will start with a trivial example showing how the "Metropolis"
    Markov chain method can be used to sample from a univariate normal
    distribution.  This will illustrate the basic facilities for
    specifying distributions, for specifying the Markov chain operations
    to use, for looking at log files, and for estimating the expectations
    of functions.

    Specifying the distribution.

    We begin by specifying the distribution we want to sample from.  This
    is done by giving a formula for the "energy function", which is minus
    the log of the probability density, plus any arbitrary constant.  To
    sample from a univariate normal distribution for a variable called
    "x", with mean 10 and variance 1, we can type the following command to
    the Unix command interpreter (the shell):

        > dist-spec nlog "(x-10)^2/2"

    Here, "nlog" is the name of a "log file", in which this specification
    is saved, and in which the results of sampling will later be stored.
    The energy formula "(x-10)^2/2" is minus the log of the probability
    density for the desired distribution over "x", with the constant term
    Log(2*Pi)/2 omitted, as it is not necessary for most purposes.  Such
    formulas must usually be put in quotes, since some characters such as
    parentheses will otherwise have special meaning to the shell.  See
    formula.doc for the syntax of formulas, which is fairly conventional,
    except perhaps that function names must be capitalized (eg, "Sin").

    Since only a single variable, "x", is mentioned in the specification,
    this is a univariate distribution.  We could have mentioned other
    "state variables" in the specification, as in the example in the next
    section (see Ex-dist-g.doc).  Valid state variable names start with
    one of "u", "v", "w", "x", "y", or "z", which may optionally be
    followed by a single digit.  See dist-spec.doc for further details.

    The density functions for some distributions, including the normal,
    are pre-defined.  We could have used the following command instead of
    the one above:

        > dist-spec nlog "Normal(x,10,1)"

    The same result can also be obtained as follows:

        > dist-spec nlog "x ~ Normal(10,1)"

    This last form is particularly useful for specifying Bayesian models
    (see Ex-bayes.doc).

    Sampling using Metropolis updates.

    After specifying the distribution, we specify what Markov chain
    operations should be used to sample from this distribution, using a
    command such as:

        > mc-spec nlog metropolis 1

    This specifies that each Markov chain iteration should consist of a
    single Metropolis operation, with a "stepsize" of 1.  In a Metropolis
    operation, a new state is proposed by randomly drawing from the normal
    distribution centred at the current state, with standard deviation
    given by the stepsize.  This proposed state is then accepted or
    rejected based on the change in the energy (ie, on the change in
    probability density).  If the proposed state is rejected, the new
    state is the same as the old state.

    This Markov chain specification is saved in the log file, after the
    distribution specification.  We could also specify an initial state
    for the Markov chain (see dist-initial.doc), but here we will let the
    initial state default to x=0.  To actually do some Markov chain
    iterations, we use a command such as:

        > dist-mc nlog 1000

    This performs 1000 Markov chain updates, as specified by the last
    mc-spec command, and saves the state after each update in the log
    file.  This takes only a fraction of a second on our machine, but
    Markov chain sampling runs for more difficult problems can take much
    longer, so one would often wish to run the 'dist-mc' command in the
    background, by putting an "&" at the end of the command line.

    If we later decided that we wanted to continue Markov chain sampling
    for more iterations, we could just use another 'dist-mc' command with
    a larger iteration limit.  For example, the command

        > dist-mc nlog 10000

    would produce another 9000 iterations, for a total of 10000.  We could
    issue another mc-spec command before this, if we wished to use
    different Markov chain operations for these further iterations.

    Checking how well the sampling worked.

    After or during the Markov chain sampling, we can look at how the
    state has changing during the run, and at certain properties of the
    Markov chain methods.

    The 'dist-display' command lets us see the state at any given
    iteration.  For example, the state at iteration 10 of the Markov chain
    run above might look like this:

        > dist-display nlog 10

        STATE VARIABLES IN FILE "nlog" WITH INDEX 10

            x = 1.45784

    However, what you see might not be exactly the same as this, due to
    differences in random number generators or in floating-point roundoff
    errors.  If the iteration number is omitted, 'dist-display' shows the
    last iteration.

    The 'dist-plt' command is usually more useful in getting a picture of
    how well the chain is sampling.  The following command will show how
    the state changes over the course of the run:

        > dist-plt t x nlog | plot

    If 'plot' is an appropriate plotting program, this will display a
    graph of the state variable "x" vs. the iteration number (the "t"
    quantity).  From this plot, you will likely see that up to about
    iteration 50, the values of "x" are not typical of those seen later in
    the run.  These early iterations should be discarded when estimating
    functions of state (see below).

    After iteration 50, the chain seems to move around the region that has
    high probability fairly rapidly.  It should therefore be possible to
    estimate expectations of functions of the state with reasonable
    accuracy.  If instead, the chain moved very slowly, it would be
    necessary to run it for many more iterations, or to use better Markov
    chain operations.

    Other interesting quantities can also be plotted using 'dist-plt'.
    For instance, the energy can be monitored with the following command:

        > dist-plt t E nlog | plot

    The change in energy on the basis of which the Metropolis proposals
    were accepted or rejected can be examined as follows:

        > dist-plt t D nlog | plot-points

    Here, it is best if the plot program used is set up to plot individual
    points, rather than lines.  For this chain, the energy difference is
    often less than one, so many proposals will be accepted.  If instead
    the energy change was almost always large, it would be necessary to
    reduce the "stepsize" argument following the "metropolis" operation in
    'mc-spec'.  On the other hand, if the energy change is usually very
    close to zero, a larger stepsize would produce better sampling.

    Other quantities that can be plotted are documented in quantities.doc,
    mc-quantities.doc, and dist-quantities.doc.

    Estimating expectations of functions.

    Finally, we can use the states from the Markov chain to estimate the
    expectations of functions of this state.  Since we decided above that
    the first 50 iterations were not typical of the chain's equilibrium
    distribution, we will use only states after these when estimating
    expectations.

    The 'dist-est' command is one way of estimating expectations.  The
    following command estimates the expectation of the "x" itself:

        > dist-est x nlog 51:

        Number of sample points: 950

        Estimates for x:

          Mean:    10.1496  (standard error 0.0336591)
          Std.dev: 1.03744

        NOTE: The standard errors assume points are independent

    The arguments to 'dist-est' are a formula for what we want to
    estimate, the log file for the run, and the range of iterations to use
    from that log file (here, from 51 on).  The output (which might be a
    bit different when you run the programs) gives the estimated
    expectation (mean) for the requested function of state, along with its
    estimated standard deviation.

    A standard error is also given for the estimated mean, but it is valid
    only if the points are independent, which is generally not true if
    they were obtained using Markov chain sampling.  In the output above,
    the estimate differs from the true mean of 10 by over four times the
    standard error, which illustrates that the standard error cannot be
    trusted when the points are not independent.

    Ideally, the 'dist-est' program would automatically compensate for
    this lack of independence, and give correct standard errors, but it
    doesn't yet.  If you are interested in the expectation of a state
    variable, however, you can get correct standard errors using the
    'series' program, which adjusts the standard errors based on estimated
    autocorrelations, if told how far out to look.  The following command
    illustrates the procedure:

        > dist-tbl x nlog 51: | series msac 15

        Number of realizations: 1  Total points: 950

        Mean: 10.149550  S.E. from correlations: 0.096872

        Standard deviation: 1.037442

          Lag  Autocorr.  Cum. Corr.

            1   0.793346    2.586691
            2   0.620884    3.828460
            3   0.492481    4.813421
            4   0.406282    5.625986
            5   0.332223    6.290432
            6   0.261235    6.812903
            7   0.195664    7.204231
            8   0.162258    7.528748
            9   0.129544    7.787835
           10   0.104300    7.996435
           11   0.079715    8.155866
           12   0.055909    8.267684
           13   0.032055    8.331795
           14   0.011870    8.355534
           15  -0.020351    8.314831

    The 'dist-tbl' command produces a list of values for "x" at iterations
    from 51 on.  This list is piped into the 'series' program.  With the
    options "msac 15", the 'series' program calculates estimates for the
    mean and standard deviation, which match those calculated above by
    'dist-est'.  It also estimates the autocorrelations for "x" at lags up
    to 15.  From these autocorrelation estimates, 'series' finds the
    "cumulative correlations" at the various lags, which are defined as
    one plus twice the sum of the autocorrelations up to that lag.  If the
    autocorrelations from that point on are approximately zero, as seems
    to be the case above, the cumulative correlations to that point will
    approximate the autocorrelation time for the quantity, which is the
    factor by which the effective sample size is less than the number of
    points used.  In calculating the standard error for the mean, 'series'
    assumes that this is the case, and it therefore divides the sample
    size by the cumulative correlation for the last lag requested when
    computing the standard error.

    In the example above, the standard error calculated by 'series'
    appears to be realistic.  This will be true only if the maximum lag
    given as an argument to series is appropriate, however.  This maximum
    lag should be large enough that the autocorrelations at larger lags
    are close to zero, but not too much larger than this, since including
    autocorrelations at many lags introduces extra noise into the estimate
    of the autocorrelation time.

    An alternative approach is to use 'series' to find out the lag past
    which the autocorrelations are almost zero, and then use 'dist-est' to
    find the expectation, telling it to look at iterations separated by
    that lag.  This is somewhat wasteful of data, but is at present the
    only easy way to get correct standard errors when estimating some
    function of state rather than a state variable itself.  For example,
    from the output of 'series' above, it seems that autocorrelations for
    "x" are almost zero at about lag 10.  This is consistent with states
    at that lag being almost independent (although lack of correlation
    does not guarantee independence).  On that basis, we might estimate
    the expectation of Sin(x) with the following command:

        > dist-est "Sin(x)" nlog 51:%10

        Number of sample points: 95

        Estimates for Sin(x):

          Mean:    -0.358403  (standard error 0.0654498)
          Std.dev: 0.637926

        NOTE: The standard errors assume points are independent

    These results are consistent with the true expectation of Sin(x) with
    respect to the Normal(10,1) distribution, which is -0.330.