meta.analysis.gibbs <- function( mu.0, sigma2.mu.0, nu.0, sigma2.sigma.0, 
  mu.initial, sigma2.initial, theta.initial, y, V, M, B ) {

  k <- length( y )

  mu <- rep( 0, M + B + 1 )

  sigma2 <- rep( 0, M + B + 1 )

  theta <- matrix( 0, M + B + 1, k )

  mu[ 1 ] <- mu.initial

  sigma2[ 1 ] <- sigma2.initial

  theta[ 1, ] <- theta.initial

  for ( m in 2:( M + B + 1 ) ) {

    mu[ m ] <- mu.full.conditional( mu.0, sigma2.mu.0, sigma2[ m - 1 ], 
      theta[ m - 1, ], y )

    sigma2[ m ] <- sigma2.full.conditional( nu.0, sigma2.sigma.0, 
      mu[ m ], theta[ m - 1, ], y )

    theta[ m, ] <- theta.full.conditional( mu[ m ], sigma2[ m ], y, V )

    if ( m %% 1000 == 0 ) print( m )

  }

  return( cbind( mu, sigma2, theta ) )

}

mu.full.conditional <- function( mu.0, sigma2.mu.0, sigma2.current, 
  theta.current, y ) {

  k <- length( y )

  k.0 <- sigma2.current / sigma2.mu.0

  theta.bar <- mean( theta.current )

  mu.k <- ( k.0 * mu.0 + k * theta.bar ) / ( k.0 + k )

  sigma2.k <- sigma2.current / ( k.0 + k )

  mu.star <- rnorm( n = 1, mean = mu.k, sd = sqrt( sigma2.k ) )

  return( mu.star )

}

sigma2.full.conditional <- function( nu.0, sigma2.sigma.0, 
  mu.current, theta.current, y ) {

  k <- length( y )

  nu.k <- nu.0 + k

  v <- mean( ( theta.current - mu.current )^2 )

  sigma2.k <- ( nu.0 * sigma2.sigma.0 + k * v ) / ( nu.0 + k )
 
  sigma2.star <- rsichi2( 1, nu.k, sigma2.k )

  return( sigma2.star )

}

rsichi2 <- function( n, nu, sigma2 ) {

  sigma2.star <- 1 / rgamma( n, shape = nu / 2, 
    rate = nu * sigma2 / 2 )

  return( sigma2.star )

}

rsichi2.2 <- function( n, nu, sigma2 ) {

  sigma2.star <- nu * sigma2 / rchisq( n, nu )

  return( sigma2.star )

}

theta.full.conditional <- function( mu.current, sigma2.current, y, V ) {

  k <- length( y )

  theta.star <- ( V * mu.current + sigma2.current * y ) / 
    ( V + sigma2.current )

  sigma2.star <- V * sigma2.current / ( V + sigma2.current )

  theta.sim <- rnorm( n = k, mean = theta.star, 
    sd = sqrt( sigma2.star ) )

  return( theta.sim )

}

mu.0 <- 0.0

sigma2.mu.0 <- 100^2

nu.0 <- 0.001

sigma2.sigma.0 <- 1.53

mu.initial <- 1.45

sigma2.initial <- 1.53

theta.initial <- c( 1.92, 1.94, 1.53, 1.84, 1.69, -0.252 )

y <- c( 2.77, 2.50, 1.84, 2.56, 2.32, -1.15 )

V <- c( 1.65, 1.31, 2.34, 1.67, 1.98, 0.90 )^2

M <- 100000

B <- 1000

mcmc.data.set <- meta.analysis.gibbs( mu.0, sigma2.mu.0, nu.0, 
  sigma2.sigma.0, mu.initial, sigma2.initial, theta.initial, 
  y, V, M, B )

mcmc.data.set <- cbind( mcmc.data.set[ , 1:2 ], 
  sqrt( mcmc.data.set[ , 2 ] ), mcmc.data.set[ , 3:8 ] )

apply( mcmc.data.set[ 1001:101001, ], 2, mean )

        mu     sigma2                                                        
1.33013835 2.24106295 1.12196766 1.68639681 1.67526967 1.38514567 1.62389213 
1.51615795 0.09356775 

apply( mcmc.data.set[ 1001:101001, ], 2, sd )

       mu    sigma2                                                  
0.9042468 4.4707971 0.9910910 1.1576621 1.0311309 1.2381000 1.1391841 
1.1917662 0.9944885 

mu.star <- mcmc.data.set[ 1001:101001, 1 ]

sum( mu.star > 0 ) / length( mu.star )

sigma.star <- mcmc.data.set[ 1001:101001, 3 ]

par( mfrow = c( 2, 1 ) )

hist( mu.star, nclass = 100, probability = T, xlab = 'mu',
  main = '' )

hist( sigma.star[ sigma.star < 5 ], nclass = 100, probability = T, 
  xlab = 'sigma', main = '' )