MODELING REAL-VALUED DATA WITH A T-DISTRIBUTION
One can also define models for the unconditional distribution of one
or more target values, without reference to any input values. Here,
I will show how one can sample from the posterior distribution for the
parameters of a univariate t-distribution. This example also shows
how to specify different stepsizes for different state variables.
Specifying the model.
Here is a specification for a model in which the degrees of freedom of
the t-distribution are fixed (at two), but the location and scale
parameters are unknown:
> dist-spec tlog1 d=2 "u ~ Normal(0,10^2) + w ~ Normal(0,1)" \
"w + [(d+1)/2] Log { 1 + [(t-u)/Exp(w)]^2/d }"
The "d=2" argument defines a constant (the degrees of freedom), which
may be used in the formulas for the prior and the likelihood. The
location parameter, u, is given a normal prior with mean zero and
standard deviation 10. The log of the scale parameter, w, is given a
normal prior with mean zero and standard deviation 1. The software
does not have a built-in t-distribution, so the likelihood has to be
specified by explicitly writing a formula for minus the log of the
probability of the target value (t) for given values of the model
parameters. This formula must include all terms that depend on the
model parameters, but for sampling from the posterior it is not
necessary to include constant terms. (However, it is necessary to
include all terms in the likelihood if the marginal likelihood for the
model is to be found using Annealed Importance Sampling.)
Specifying the data source.
We can specify that the data comes from the file "tdata" as follows:
> data-spec tlog1 0 1 / tdata .
This says that there are no input values and one target value, and
that the data comes from the file "tdata" (note that one must say
where the input values come from even though there are zero of them,
but that one can then say that the targets come from the same place
using the "." specification).
The contents of "tdata" are as follows:
0.9
1.0
1.1
6.9
7.0
7.1
Sampling with the Metropolis algorithm.
We can now try sampling using the Metropolis algorithm with a stepsize
of 1, repeated 10 times each iteration:
> mc-spec tlog1 repeat 10 metropolis 1
> dist-mc tlog1 5000
This takes 0.16 seconds on the system used (see Ex-test-system.doc).
We can look at plots such as the following to assess convergence:
> dist-plt t uw tlog1 | plot
Equilibrium appears to have been nearly reached within 50 iterations
(or less). We can now see a picture of the posterior distribution
using a command such as
> dist-plt u w tlog1 50:%3 | plot-points
This produces a scatterplot for the values of the two state variables
at those iterations from 50 on that are divisible by three. Looking
at only every third state reduces the number of points that will be
superimposed, due to rejections of the Metropolis proposals. (If some
points are superimposed, the scatterplot might be misleading.)
If you produce such a plot, you should see a tooth-shaped distribution
whose main mass is around u=3.9 and w=1, with two roots descending at
u=1 and u=7.
Varying the stepsizes.
The rejection rate of the metropolis operations with stepsize of 1 can
be estimated as follows:
> dist-tbl r tlog1 | series m
Number of realizations: 1 Total points: 5000
Sample mean: 0.623220
The efficiency of this chain at estimating the mean of u can be
assessed as follows:
> dist-tbl u tlog1 50: | series mac 10
Number of realizations: 1 Total points: 4951
Sample mean: 3.883301 S.E. from correlations: 0.056870
Lag Autocorr. Cum. Corr.
1 0.641608 2.283216
2 0.408967 3.101149
3 0.260927 3.623004
4 0.176241 3.975486
5 0.119894 4.215273
6 0.074613 4.364499
7 0.028414 4.421327
8 0.006842 4.435012
9 -0.018193 4.398626
10 -0.031221 4.336183
An inefficiency factor of 4.3 is not bad, but we still might try to
improve sampling by using a different stepsize. We can use a stepsize
twice as big by changing the 'mc-spec' command as follows (with other
commands the same, except for using tlog2 as the log file):
> mc-spec tlog2 repeat 10 metropolis 2
This produces a higher rejection rate:
> dist-tbl r tlog2 | series m
Number of realizations: 1 Total points: 5000
Sample mean: 0.815860
Despite this, the efficiency of sampling is improved:
> dist-tbl u tlog2 50: | series mac 10
Number of realizations: 1 Total points: 4951
Sample mean: 3.807873 S.E. from correlations: 0.047144
Lag Autocorr. Cum. Corr.
1 0.482591 1.965182
2 0.244778 2.454738
3 0.117978 2.690693
4 0.078094 2.846881
5 0.052391 2.951662
6 0.016818 2.985299
7 0.009607 3.004513
8 0.003207 3.010927
9 -0.002763 3.005401
10 -0.007795 2.989811
From the scatterplot of the posterior distribution, the larger
stepsize seems appropriate for u, but it seems too big for w. We can
specify different stepsizes for u and w as follows:
> dist-stepsizes tlog12 w=1 u=2
> mc-spec tlog12 repeat 10 metropolis 1
The 'dist-stepsizes' command lets you specify individual stepsizes for
the state variables. These values are multiplied by the stepsize
specified for the "metropolis" operation before being used, so you can
change the overall size of the steps by changing this argument, while
keeping the relative stepsizes the same.
Using twice as big a stepsize for u as for w seems to work well, as
seen by looking at the autocorrelations for u:
> dist-tbl u tlog12 50: | series mac 10
Number of realizations: 1 Total points: 4951
Sample mean: 3.872209 S.E. from correlations: 0.038725
Lag Autocorr. Cum. Corr.
1 0.293898 1.587795
2 0.080062 1.747919
3 0.022408 1.792735
4 0.038690 1.870116
5 0.004691 1.879497
6 -0.007760 1.863978
7 0.005857 1.875691
8 0.029724 1.935138
9 0.007939 1.951016
10 -0.008345 1.934326
The autocorrelations have been reduced to the point where the estimate
based on the points from this chain are a factor of only about 1.9
times less efficient than an estimate based on independent points.