##### sampling from a Dirichlet distribution

rdirichlet <- function( n.sim, alpha ) {    

  J <- length( alpha )

  theta <- matrix( 0, n.sim, J )

  for ( j in 1:J ) {

    theta[ , j ] <- rgamma( n.sim, alpha[ j ], 1 )

  }

  theta <- theta / apply( theta, 1, sum )

  return( theta )

}

##### checking the rdirichlet function

alpha <- c( 5.0, 1.0, 2.0 )

alpha.0 <- sum( alpha )

test <- rdirichlet( 100000, alpha ) # 7 seconds (1.6 Unix GHz)

print( apply( test, 2, mean ) )

# [1] 0.6258544 0.1247550 0.2493905 (here and below, you should 
# get results like the numbers given, apart from Monte Carlo noise)

print( alpha / alpha.0 )

# [1] 0.625 0.125 0.250

print( apply( test, 2, var ) )

# [1] 0.02603293 0.01216358 0.02071587

print( alpha * ( alpha.0 - alpha ) / ( alpha.0^2 * ( alpha.0 + 1 ) ) )

# [1] 0.02604167 0.01215278 0.02083333

print( cov( test ) )

#              [,1]         [,2]         [,3]
# [1,]  0.026032929 -0.008740319 -0.017292610
# [2,] -0.008740319  0.012163577 -0.003423259
# [3,] -0.017292610 -0.003423259  0.020715869

print( - outer( alpha, alpha, "*" ) / ( alpha.0^2 * ( alpha.0 + 1 ) ) )

#              [,1]         [,2]         [,3]
# [1,] -0.043402778 -0.008680556 -0.017361111
# [2,] -0.008680556 -0.001736111 -0.003472222 # ignore diagonals
# [3,] -0.017361111 -0.003472222 -0.006944444

##### BQQI analysis of the IHGA data

alpha.C <- c( 138.001, 77.001, 46.001, 12.001, 8.001, 4.001, 0.001, 2.001 )

alpha.E <- c( 147.001, 83.001, 37.001, 13.001, 3.001, 1.001, 1.001, 0.001 )

theta.C <- rdirichlet( 100000, alpha.C ) # 8 sec (1.6 Unix GHz)

theta.E <- rdirichlet( 100000, alpha.E ) # also 8 sec

print( post.mean.theta.C <- apply( theta.C, 2, mean ) )

# [1] 4.808015e-01 2.683458e-01 1.603179e-01 4.176976e-02 
# [5] 2.784911e-02 1.395287e-02 3.180905e-06 6.959859e-03

print( post.SD.theta.C <- apply( theta.C, 2, sd ) )

# [1] 0.0294142963 0.0261001259 0.0216552661 0.0117925465 
# [5] 0.0096747630 0.0069121507 0.0001017203 0.0048757485

print( post.mean.theta.E <- apply( theta.E, 2, mean ) )

# [1] 5.156872e-01 2.913022e-01 1.298337e-01 4.560130e-02 
# [5] 1.054681e-02 3.518699e-03 3.506762e-03 3.356346e-06

print( post.SD.theta.E <- apply( theta.E, 2, sd ) )

# [1] 0.029593047 0.026915644 0.019859213 0.012302252 
# [5] 0.006027157 0.003501568 0.003487824 0.000111565

mean.effect.C <- theta.C %*% ( 0:7 )

mean.effect.E <- theta.E %*% ( 0:7 )

mult.effect <- mean.effect.E / mean.effect.C

print( post.mean.mult.effect <- mean( mult.effect ) )

# [1] 0.8189195

print( post.SD.mult.effect <- sd( mult.effect ) )

# [1] 0.08998323

print( quantile( mult.effect, probs = c( 0.0, 0.025, 0.5, 0.975, 1.0 ) ) )

#        0%      2.5%       50%     97.5%      100% 
# 0.5037150 0.6571343 0.8138080 1.0093222 1.3868332 

print( mean( mult.effect < 1 ) )

# [1] 0.9706

pdf( "bqqi-mult-effect.pdf" )

plot( density( mult.effect, n = 2048 ), type = 'l', 
  cex.lab = 1.25, cex.axis = 1.25, main = '',
  xlab = 'Multiplicative Treatment Effect' )

dev.off( )