m = 10^5;  n = 4;  alp = 2
 x = rgamma(n*m, alp, 1)
 DTA = matrix(x, nrow=m, byrow=T)
 s = rowSums(DTA);  u = DTA[,1]/s
 mean(u);  var(u)
 ## 0.2506025    ## approximates E(U)
 ## 0.02084331   ## approximates V(U)

 hist(u, prob=T, col="wheat", main="")
 curve(dbeta(x, alp, (n-1)*alp), n = 1001, lwd=2, col="blue", add=T)