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.