# main simulation function

rdiric_function(n,a) {
        p <- length(a)
        m <- matrix(nrow=n,ncol=p)
        for (i in 1:p) {
                m[,i] <- rgamma(n,a[i])
        }
        sumvec <- m %*% rep(1,p)
	m / as.vector(sumvec)
}

# examples
op_par(mfrow=c(2,2))

x_rdiric(5000,c(1,1,1))
plot(x[,1],x[,2],main="(1,1,1)",xlim=c(0,1),ylim=c(0,1))

x_rdiric(5000,c(1,1,1,1,1))
plot(x[,1],x[,2],main="(1,1,1,1,1)",xlim=c(0,1),ylim=c(0,1))

# mean of (0.3,0.4,0.3)
x_rdiric(5000,100*c(0.3,0.4,0.3))
plot(x[,1],x[,2],main="50*(0.3,0.4,0.3)",xlim=c(0,1),ylim=c(0,1))

x_rdiric(5000,100*c(0.5,0.4,0.1))
plot(x[,1],x[,2],main="20*(0.5,0.4,0.1)",xlim=c(0,1),ylim=c(0,1))

par(op)