# Problem 2

rdirichlet <- function( n, alpha ) {
  p <- length( alpha )
  g <- matrix( NA, n, p )
  for ( i in 1:n ) {
    for ( j in 1:p ) {
      g[ i, j ] <- rgamma( 1, rate = 1, shape = alpha[ j ] )
    }
    g[ i, ] <- g[ i, ] / sum( g[ i, ] )
  }
  return( g )
}

n1 = 727
n2 = 583
n3 = 137
n = n1 + n2 + n3

nh1 = n1 / n
nh2 = n2 / n
nh3 = n3 / n

# Standard deviation is sqrt(V/n)
v1 = nh1*(1 - nh1)/n
v2 = nh2*(1 - nh2)/n
v3 = nh3*(1 - nh3)/n
se1 = sqrt(v1/n)
se2 = sqrt(v2/n)
se3 = sqrt(v3/n)

a = c(1, -1, 0)
sigma = diag(c(v1, v2, v3))
sigma[1,2] = -nh1 * nh2 / n
sigma[2,1] = -nh1 * nh2 / n
sigma[1,3] = -nh1 * nh3 / n
sigma[3,1] = -nh1 * nh3 / n
sigma[2,3] = -nh2 * nh3 / n
sigma[3,2] = -nh2 * nh3 / n

gamma = nh1 - nh2
vgamma = t(a) %*% sigma %*% a
hat.ci = c(
  gamma - 1.96*sqrt(vgamma/n),
  gamma + 1.96*sqrt(vgamma/n)
)

alpha = c(0.01, 0.01, 0.01)
M = 500
g = rdirichlet(M, alpha)
x = g[,1] - g[,2]
sum(x > 0)