# Monte Carlo Integration
# estimate expected value of u1/(u2+u3) 
# where u1 = x1/(x1+x2+x3+x4), u2 = x2/(x1+x2+x3+x4), u3 = x3/(x1+x2+x3+x4) 


result=c()
n <- 10000 # number of simulations

#loop through different values of beta

for (beta in c(1,2,5,10,100))
{
x1 <- rgamma(n, 0.1, beta) # gamma distribution | alpha = 0.1
x2 <- rgamma(n, 0.3, beta) # gamma distribution | alpha = 0.3
x3 <- rgamma(n, 5, beta) # gamma distribution | alpha = 5
x4 <- rgamma(n, 2, beta) # gamma distribution | alpha = 2

s <- x1+x2+x3+x4
u1 <- x1/s
u2 <- x2/s
u3 <- x3/s

result <- rbind( result, c(beta, mean(u1/(u2+u3))))
}
result