    A LINEAR REGRESSION MODEL

    I generated 100 cases of synthetic data for this example, each case
    consisting of two input (predictor) values and one target (response)
    value.  The first input (i0) was randomly generated from the normal
    distribution with mean zero and standard deviation one.  The second
    input (i1) was set to the first input plus an offset that was normally
    distributed with mean zero and standard deviation 0.1.  The two inputs
    were therefore highly correlated.  The target value (t) was set to
    0.5 + 2.5*i0 - 0.5*i1 + e, where e was normally distributed noise with
    mean zero and standard deviation 0.1.  This data was stored in the
    file "rdata".

    Specifying the model.

    We can specify a Bayesian linear model for this data with a command
    such as the following:

        > dist-spec rlog.met \
        >   "u ~ Normal(0,10^2) + w0 ~ Normal(0,10^2) + w1 ~ Normal(0,10^2) \
        >      + v ~ ExpGamma(1,0.2)" \
        >   "t ~ Normal (u + w0i0 + w1i1, Exp(-v))"

    Here, the "\" characters at the ends of the lines tell the Unix shell
    that the command is continued on the next line.  The "rlog.met"
    argument of 'dist-spec' is the name of the log file to use.  The next
    argument (which takes up two lines) is the prior specification.  It
    must be put in quotes to keep the shell from interpreting the special
    characters used.  The last argument is the likelihood specification,
    which will also usually have to be quoted.

    Let's first examine the likelihood above:

        t ~ Normal (u + w0i0 + w1i1, Exp(-v))

    This says that the target value for a case, t, is modeled by a normal
    distribution with mean given by u + w0i0 + w1i1 and variance given by
    Exp(-v).  Here, u, w0, w1, and v are parameters of the model.  Model
    parameters must have names starting with "u", "v", "w", "x", "y", or
    "z", possibly followed by a single digit.  The two inputs for the case
    are called i0 and i1.  If there were more than one target, they would
    be called t0, t1, etc., but here the name "t" is used, which is a
    synonym for "t0" (similarly "i" is a synonym for "i0").

    The notation in the likelihood specification using "~" is actually an
    abbreviation for the following:

        Log(2*Pi)/2 + [t - (u+w0i0+w1i1)]^2 / [2Exp(-v)]

    which is minus the log of the normal probability density for t with
    the given mean and variance.  In general, the likelihood is specified
    by a formula giving minus the log of the probability or probability
    density for the target or targets in one case, given the values of the
    inputs.  Since the cases are assumed to be independent, minus the
    total log likelihood is found by just adding these terms for all
    cases.

    The prior specification above is constructed with similar notation:

        u ~ Normal(0,10^2) + w0 ~ Normal(0,10^2) + w1 ~ Normal(0,10^2)
          + v ~ ExpGamma(1,0.2)

    This gives each of u, w0, and w1 an independent normal prior with mean
    of 0 and standard deviation of 10 (ie, variance of 10^2).  The prior
    specification is a formula for minus the log of the prior density, so
    terms pertaining to different parameters are simply added together.
    The terms involving "~" are abbreviations for minus the log prior
    density for a known distribution.  Here, all the parameters are
    independent under the prior, but it would be possible to make one
    parameter depend on another, for instance by changing the term above
    involving w1 to w1 ~ Normal(w0,10^2).

    The "v" parameter is the log of the "precision" (inverse variance) for
    the error term in the model.  The Markov chain sampling procedures
    assume that the range of all variables being sampled is the entire
    real line.  Since the precision is constrained to be positive, we must
    represent it with some unconstrained variable that is transformed to
    produce the actual precision.  An exponential transformation is often
    convenient, and in this case is probably desirable in any case, in
    order to make the Markov chain sampling more efficient.

    Here, the prior for the precision is specified to be Gamma with shape
    parameter of 1 and scale of 0.2 (which gives a mean value of 5).
    Since we will generally want to represent Gamma variables by their
    logs, the built-in ExpGamma distribution used here is for a variable
    whose exponential has a Gamma distribution with the indicated shape
    and scale parameters.

    Specifying the data source.

    After specifying the model with 'dist-spec', we need to say where the
    data is found using 'data-spec', as follows:

        > data-spec rlog.met 2 1 / rdata .

    The arguments after the name of the log file say that there are two
    input values and one target value.  After the slash, the files where
    the inputs and targets are found are given.  The "." says that the
    targets are found in the same file as the inputs, following them on
    the same line.  See data-spec.doc for more details and other options.

    Sampling with Metropolis updates.

    We can now sample from the posterior distribution of the model
    parameters, after specifying how we want this to be done.  Here we'll
    try sampling by Metropolis updates to all variables simultaneously:

        > mc-spec rlog.met repeat 45 metropolis 0.02
        > dist-mc rlog.met 500

    Each iteration will consist of 45 Metropolis updates, done using a
    proposal distribution that adds a normally-distributed offset to each
    variable with standard deviation 0.02.  The 'dist-mc' command does 500
    such iterations (for a total of 22500 Metropolis updates).  This takes
    4.4 seconds on the system used (see Ex-system.doc).

    We can see what the average rejection rate was for these updates as
    follows:

        > dist-tbl r rlog.met | series m

        Number of realizations: 1  Total points: 500

        Mean: 0.708978

    A rejection rate of 0.7 is about right, so the stepsize of 0.02 looks
    to be about optimal.  If you had started by trying a much larger or
    smaller stepsize, you would need to adjust it by trial and error to
    get a rejection rate in the range of about 0.3 to 0.8.

    We can look at how well the chain is sampling by plotting the four
    state variables over time:

        > dist-plt t uw0w1v rlog.met | plot

    From this plot, it appears that the equilibrium distribution is
    reached by about iteration 100, and that the chain is sampling
    somewhat adequately from then on, though with a substantial degree of
    autocorrelation.  We can quantify the inefficiency in estimation of
    the log precision, v, that results from these autocorrelations as
    follows:

        > dist-tbl v rlog.met 100: | series mac 30

        Number of realizations: 1  Total points: 401

        Mean: 4.461414  S.E. from correlations: 0.026037

          Lag  Autocorr.  Cum. Corr.

            1   0.893653    2.787307
            2   0.807894    4.403095
            3   0.709395    5.821885
            4   0.621722    7.065329
            5   0.539335    8.143999
            6   0.462553    9.069106
            7   0.379245    9.827595
            8   0.311431   10.450456
            9   0.244173   10.938802
           10   0.198226   11.335254
           11   0.164614   11.664482
           12   0.126405   11.917293
           13   0.092779   12.102850
           14   0.074293   12.251436
           15   0.049729   12.350895
           16   0.037565   12.426024
           17   0.028010   12.482045
           18   0.026148   12.534341
           19   0.043167   12.620674
           20   0.062433   12.745541
           21   0.083616   12.912773
           22   0.108346   13.129465
           23   0.112400   13.354265
           24   0.111866   13.577997
           25   0.092298   13.762593
           26   0.076360   13.915313
           27   0.051526   14.018364
           28   0.038505   14.095374
           29   0.018084   14.131542
           30   0.006680   14.144903

    The cumulative correlation up to lag 30, where the autocorrelation is
    approximately zero, indicates that about 14 times as many points will
    be needed to get an estimate of a given accuracy as would be needed if
    the points were independent.  This inefficiency results from the high
    dependency between w0 and w1, which can be seen in the following plot:

        > dist-plt w0 w1 rlog.met 100: | plot-points

    This dependency is why the stepsize must be so small - proposed
    changes must be small, as otherwise the new point is likely to be
    incompatible with this dependency between w0 and w1.

    Estimates for various functions of state can be found using
    'dist-est'.  For example, we can predict the target value for input
    values of i0=1.9 and i1=2.1 as follows:

        > dist-est "u + w0*1.9 + w1*2.1" rlog.met 100:

        Number of sample points: 401

        Estimates for u + w0*1.9 + w1*2.1:

          Mean:    4.18959  (standard error 0.00150075)
          Std.dev: 0.0300526

        NOTE: The standard errors assume points are independent

    Note that since the points aren't close to being independent, the
    standard error shown is not correct.  Note also that the standard
    deviation shown is the posterior standard deviation of the linear fit
    at this point.  The prediction for a target value in a case with these
    inputs would have additional uncertainty due to the noise.

    Sampling with hybrid Monte Carlo.

    The posterior dependency between w0 and w1 makes Markov chain sampling
    somewhat difficult for this problem.  One could try to eliminate this
    dependency by a reparameterization, which would certainly be possible
    for this fairly simple problem, but in general this may be tedious or
    infeasible.  However, the inefficiencies caused by dependencies may be
    reduced by suppressing the random walk aspect of the sampling, which
    can be done by using hybrid Monte Carlo rather than the simple
    Metropolis updates used above.

    After 'dist-spec' and 'data-spec' commands that are the same as those
    described above (apart from the name of the log file), we can sample
    using hybrid Monte Carlo with commands such as the following:

        > mc-spec rlog.hmc heatbath hybrid 25 0.01
        > dist-mc rlog.hmc 500

    One hybrid Monte Carlo update is done each iteration, consisting of a
    "heatbath" operation that picks new values for the momentum variables,
    followed by a dynamical trajectory of 25 leapfrog steps, with a
    stepsize of 0.01.  This stepsize was picked by trial and error, and
    results in a rejection rate of about 0.13.  In general, the stepsize
    for this standard form of hybrid Monte Carlo should chosen so that the
    rejection rate is somewhere between about 0.05 and 0.25.  The 500
    iterations take 4.7 seconds on the system used (see Ex-system.doc),
    about the same as the Metropolis run above.

    From examining plots of the state variables over time, it seems that
    the hybrid Monte Carlo chain reaches equilibrium much more quickly
    than the Metropolis chain, and samples more efficiently once
    equilibrium is reached.  For example, here are the autocorrelations
    for the precision parameter, which we also looked at for the
    Metropolis chain:

        > dist-tbl v rlog.hmc 20: | series mac 5

        Number of realizations: 1  Total points: 481

        Mean: 4.451747  S.E. from correlations: 0.005397

          Lag  Autocorr.  Cum. Corr.

            1  -0.077573    0.844854
            2   0.048169    0.941192
            3  -0.059657    0.821878
            4   0.021492    0.864861
            5  -0.071308    0.722244

    The estimated autocorrelations above do not differ significantly from
    zero.  It therefore appears that hybrid Monte Carlo is about ten times
    more efficient than the simple Metropolis algorithm for this problem.

    Hybrid Monte Carlo with windows.

    Many other Markov chain methods could be applied to this problem, but
    I will mention just one more - the variation of hybrid Monte Carlo in
    which acceptance is based on "windows" of states at the start and at
    the end of a trajectory.  This variation is described in my paper on
    "An improved acceptance procedure for the hybrid Monte Carlo
    algorithm" (Journal of Computational Physics, vol. 111, pp. 194-203,
    1994), and more briefly in my book on Bayesian Learning for Neural
    Networks.

    This method can be applied to this problem as follows (after suitable
    'dist-spec' and 'data-spec' commands):

        > mc-spec rlog.whmc heatbath hybrid 25:2 0.01
        > dist-mc rlog.whmc 500

    The only difference compared to the specification for standard hybrid
    Monte Carlo is the ":2" after the number of leapfrog steps (25).  This
    specifies that windows of two states at the start and end of the
    trajectory should be used.  Acceptance is based on the difference in
    the total probability of all states in the window at the start and in
    the window at the end.  The next state is picked from the end window,
    if the trajectory was accepted, or from the start window, if the
    trajectory was rejected, with state selected from the appropriate
    window according to their relative probabilities.

    Note that when windows are used the state may change even when there
    was a "rejection", as long as one of the states in the start window is
    of comparable probability to the starting state.  This provides some
    insurance against getting stuck at a point where the rejection rate is
    very high.  The use of windows also tends to reduce the rejection
    rate.  In this problem, the rejection rate with windows is only 0.05,
    less than the rate of 0.13 seen above for the standard hybrid Monte
    Carlo algorithm, while the time is only slightly greater (due to a
    slight overhead from using windows).  The benefit is small here, since
    the standard hybrid Monte Carlo method was doing very well in any
    case, but use of windows can make a bigger difference for some other
    problems, and is almost cost-free if the window size is a small
    fraction of the trajectory length.