#Example. Estimating the speed of light using Simon Newcomb's measurement

 library(MASS)
 newcomb
 [1]  28 -44  29  30  24  28  37  32  36  27  26  28  29  26  27  22  23  20  25
[20]  25  36  23  31  32  24  27  33  16  24  29  36  21  28  26  27  27  32  25
[39]  28  24  40  21  31  32  28  26  30  27  26  24  32  29  34  -2  25  19  36
[58]  29  30  22  28  33  39  25  16  23
 hist(newcomb)
 hist(newcomb,20)
#postscript("fig4.ps", horizontal = FALSE, onefile = FALSE,
#          paper = "special", width = 10, height = 6)
 hist(newcomb,30,main="",xlab="")
#dev.off()

 t.test(newcomb)

################################################################################
        #One Sample t-test

data:  newcomb 
t = 19.8178, df = 65, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 0 
95 percent confidence interval:
 23.57059 28.85365 
sample estimates:
mean of x 
 26.21212 

#################################################################################
	#Exact posterior interval for a new measurement 
     
y <- newcomb
n <- length(y)
c(mean(y)+qt(0.025,df=n-1)*sd(y)*sqrt(1+1/n),mean(y)-qt(0.025,df=n-1)*sd(y)*sqrt(1+1/n))
[1]  4.590262 47.833980
##################################################################################
      #Simulation experiment


s2v<-numeric(1000)
muv<-numeric(1000)
ynew <- numeric(1000)

for(i in 1:1000){
  s2<-65*10.8^2/rchisq(1,df=65)
  mu<-rnorm(1,26.2,sd=sqrt(s2/66))
  s2v[i]<-s2
  muv[i]<-mu
  ynew[i] <- rnorm(1,muv[i],sqrt(s2v[i]))
}


plot(muv,s2v)
plot(density(muv))

#Posterior interval for mu
quantile(muv,c(0.025,0.975))
#Posterior interval for a new observation
quantile(ynew,c(0.025,0.975))