Speed of light data revisited.

gistfile1.txt

newcomb <-
c(28,26,33,24,34,-44,27,16,40,-2,
29,22,24,21,25,30,23,29,31,19,
24,20,36,32,36,28,25,21,28,29,
37,25,28,26,30,32,36,26,30,22,
36,23,27,27,28,27,31,27,26,33,
26,32,32,24,39,28,24,25,32,25,
29,27,28,29,16,23)

# Data as a list
a.dat <- list( y=newcomb, I=length(newcomb) )
# Inits as a list of lists
a.ini <- list( list( alpha=20, sigma=8 ),
               list( alpha=23, sigma=10 ),
               list( alpha=26, sigma=12 ) )

# Names of the parameters to monitor
a.par <- c("alpha","sigma", "smallest","mn")

cat( "model
      {
      for( i in 1:I )
         {
         y[i] ~ dt(mu[i],tau,2)
         mu[i] <- alpha
     ##generates predicted values:      
         #y.pred[i] ~ dnorm(mu[i],tau)
          y.pred[i] ~ dt(mu[i],tau,2)
         }
      alpha ~ dnorm(0, 1.0E-6)
      sigma ~ dunif(0,100)
      tau  <- 1/pow(sigma,2)
      smallest <- min(y.pred)
      mn <- mean(y.pred)
      }",
      file="tdistrlight2.jag" )

library(rjags)
a.mod.t <- jags.model( file = "tdistrlight2.jag",
                   data = a.dat,
               n.chains = 3,
                  inits = a.ini,
                n.adapt = 20000,
                quiet=T)

a.res.t <- coda.samples( a.mod.t,
                         var = a.par,
                      n.iter = 40000,
                        thin = 10 )

output<-as.data.frame(as.matrix(a.res.t))

hist(newcomb,freq=FALSE,breaks=50,ylim=c(0,0.25))
lines(density(output$smallest))
lines(density(output$mn),col="red",lwd=4)