gamma.sim <- function( m, alpha, beta, n.sim, seed ) { set.seed( seed ) theta.out <- matrix( 0, n.sim, 2 ) for ( i in 1:n.sim ) { theta.sample <- rgamma( m, shape = alpha, scale = 1 / beta ) theta.out[ i, 1 ] <- mean( theta.sample ) theta.out[ i, 2 ] <- sqrt( var( theta.sample ) / m ) } return( theta.out ) } m <- 1000 alpha <- 29.001 beta <- 14.001 n.sim <- 500 seed <- 9626954 theta.out <- gamma.sim( m, alpha, beta, n.sim, seed ) theta.out[ 1:10, ] # [1,] 2.050852 0.01217430 # [2,] 2.081629 0.01209275 # [3,] 2.067799 0.01243688 # [4,] 2.084047 0.01201150 # [5,] 2.063885 0.01219349 # [6,] 2.093904 0.01214949 # [7,] 2.081116 0.01203641 # [8,] 2.056941 0.01220216 # [9,] 2.086504 0.01242793 # [10,] 2.073801 0.01210534 # postscript( "gamma-sim1.ps" ) theta.bar <- theta.out[ , 1 ] qqnorm( ( theta.bar - mean( theta.bar ) ) / sqrt( var( theta.bar ) ) ) abline( 0, 1 ) # dev.off( ) # null device # 1 truth <- alpha / beta theta.bar.SE <- theta.out[ , 2 ] qnorm( 0.025 ) # [1] -1.959964 sum( ( theta.bar - 1.96 * theta.bar.SE < truth ) * ( truth < theta.bar + 1.96 * theta.bar.SE ) ) / n.sim # [1] 0.954 print( theta.bar <- gamma.sim( m, alpha, beta, 1, seed ) ) # [,1] [,2] # [1,] 2.050852 0.0121743 print( theta.bar <- gamma.sim( 569358, alpha, beta, 1, seed ) ) # [,1] [,2] # [1,] 2.071157 0.0005097478