Simulating data: learning about the
                        degrees of freedom of a
                        t-distribution

This example uses simulated data: we generate n = 1000 observations from a (standard) t distribution with d = 4 degrees of freedom:

       model {
          d <- 4
          for (i in 1 : 1000) {
             y[i] ~ dt(0, 1, d)
          }
       }

We let the sampler run for a while (just to check that we're simulating from the correct distribution) and then we save a set of simulated data using the 'save state' facility (select Save State from the Model menu).

Check: should have mean 0, variance 2, etc.

      mean   sd   MC_error   val2.5pc   median   val97.5pc   start   sample
   y[1]   0.01129   1.383   0.01415   -2.645   -0.01056   2.843   1001   10000
   

   
To analyse the simulated data: first we try a model for learning about the degrees of freedom as a continuous quantity:

Model:

       model {
          for (i in 1:1000) {
             y[i] ~ dt(0, 1, d)
          }
          d ~ dunif(2, 100)         # degrees of freedom must be at least two
       }

      mean   sd   MC_error   val2.5pc   median   val97.5pc   start   sample
   d   4.282   0.4398   0.01448   3.493   4.252   5.217   1001   10000



Now we attempt to model the degrees of freedom parameter with a discrete uniform prior on {2, 3, 4, ..., 50}. The sampler soon converges to d = 4 but mixes poorly.


Model:

       model {
          for (j in 1:49) {
             p[j] <- 1 / 49
             d.value[j] <- j + 1
          }
          for (i in 1:1000) {
             y[i] ~ dt(0, 1, d)
          }
          K ~ dcat(p[])
          d <- d.value[K]
       }
      
      mean   sd   MC_error   val2.5pc   median   val97.5pc   start   sample
   d   4.204   0.4102   0.03144   4.0   4.0   5.0   1001   10000




We should get better mixing if we specify the prior for d on a finer grid, e.g. {2.0, 2.1, 2.2, ..., 6.0}:

Model:

    model {
       for (j in 1:41) {
          p[j] <- 1 / 41
          d.value[j] <- 2 + (j - 1) * 0.1
       }
       for (i in 1:1000) {
          y[i] ~ dt(0, 1, d)
       }
       K ~ dcat(p[])
       d <- d.value[K]
    }

      mean   sd   MC_error   val2.5pc   median   val97.5pc   start   sample
   d   4.277   0.4162   0.01194   3.5   4.3   5.1   1001   10000



Simulated data:

   click to opensimulated t4 data