# unbuffer the output

beta.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 <- rbeta( m, alpha, beta )

    theta.out[ i, 1 ] <- mean( theta.sample )

    theta.out[ i, 2 ] <- sqrt( var( theta.sample ) / m )

  }

  return( theta.out )

}

m <- 100

alpha <- 76.5

beta <- 353.5

n.sim <- 10000

seed <- 1

system.time( {

  print( Sys.time( ) )

  theta.out <- beta.sim( m, alpha, beta, n.sim, seed )

  }

)

# [1] "2018-02-27 10:27:36 PST"
#    user  system elapsed 
#    0.38    0.02    0.39

# start Speccy (if you have it) and Resource Monitor (CPU)

# Intel Core i7-7820HQ CPU @ 2.90GHz (1 thread: not parallelized)

n.sim <- 1000000

system.time( {

  print( Sys.time( ) )

  theta.out <- beta.sim( m, alpha, beta, n.sim, seed )

  }

)

# [1] "2018-02-27 10:28:11 PST"
#    user  system elapsed 
#   37.28    0.02   38.54

theta.out[ 1:5, ]

#           [,1]        [,2]
# [1,] 0.1817399 0.001686891
# [2,] 0.1756339 0.001698506
# [3,] 0.1803205 0.001720211
# [4,] 0.1765938 0.001817393
# [5,] 0.1784614 0.001838266

print( truth <- alpha / ( alpha + beta ) )

# [1] 0.177907

###################################################################

theta.bar <- theta.out[ , 1 ]

qqnorm( ( theta.bar - mean( theta.bar ) ) / sd( theta.bar ),
  xlab = "Quantiles of Standard Normal", main = "", pch = 20 )

abline( 0, 1, lwd = 2, col = 'red' )

###################################################################

theta.bar.SE <- theta.out[ , 2 ]

sum( ( theta.bar - 1.96 * theta.bar.SE < truth ) *
  ( truth < theta.bar + 1.96 * theta.bar.SE ) ) / n.sim

[1] 0.947062

###################################################################

seed <- 5

print( theta.bar <- beta.sim( m, alpha, beta, 1, seed ) )

          [,1]        [,2]
[1,] 0.1794077 0.001893772

###################################################################