set.seed(123) nrep<- 105000 # Total number of MC replications nb<- 5000 # Number of observations for the "Burn-in" n<- 10 # Sample size tau<- array(,nrep) # Set up vectors for storing results mu<- array(,nrep) y<- rnorm(n,mean=1,sd=1) # Create a sample of data: N[1,1] # True values of Mu and Tau are each 1 ybar<- mean(y) yy<- sum(y^2) lambda<- 1/(0.5*n*var(y)) ttau<- rgamma(1, shape = n/2, rate = lambda) #initialize Tau #START OF THE MCMC LOOP: for (i in 1:nrep) { mmu<-rnorm(1,mean = ybar,sd = 1/sqrt(n*ttau)) scale<- (0.5*(yy+n*mmu^2-2*n*mmu*ybar)) ttau<- rgamma(1, shape=n/2, scale) tau[i]<- ttau mu[i]<- mmu } # END OF THE MCMC LOOP # Drop the first "nb" repetitions for the "Burn-in" # We have 100,000 values for the marginal posteriors # Let's see if the results seem to be accurate: nb1<-nb+1 taub<-tau[nb1:nrep] mub<- mu[nb1:nrep] # The marginal posteriors for Mu and Tau should be Student-t (n-1), and Gamma, respectively library(moments) summary(mub) ; var(mub) ybar # The mean of the poterior for Mu should be ybar ( = 1.0746) skewness(mub) # the skewness of Student-t is zero kurtosis(mub) # The EXCESS kurtosis for Student-t (n-1) is 6/(n-5)=1.2; so kurtosis = 4.2 summary(taub) ; var(taub) skewness(taub) # the skewness of Gamma is (2/sqrt(shape))= (2/sqrt(n/2)=0.8944 kurtosis(taub) # excess kurtosis for Gamma is (6/shape) = 6/(n/2)=1.2 hist (mub, prob=T,col=2, main="Marginal Posterior for Mu: Student-t", xlab="Mu", ylab="Marginal PDF for Mu") hist (taub,prob=T,col=4, main="Marginal Posterior for Tau: Gamma", xlab="Tau", ylab="Marginal PDF for Tau")