In order to simulate from a mixture of Dirichlet distributions, we can write:
alpha_s=c(1,1,1)
library(DirichletReg)
Simulation <- function(r, alpha_s){
z <- matrix(NA, nrow = r, ncol = length(alpha_s))
for (i in 1:r) {
x1 <- rdirichlet(n=1, alpha=alpha_s)
x1 <- as.numeric(x1)
x2 <- rdirichlet(n=1, alpha=alpha_s)
x2 <- as.numeric(x2)
w <- rbeta(1,3,3)
z[i,] <- w * x1 + (1-w) * x2
}
return(z)
}
# histograms
f <- Simulation(r=1000, alpha_s = c(1,1,1))
colnames(f) <- c('Z1','Z2','Z3')
hist(f[,1], col='gold', breaks=20, main='histogram of z1', xlab='z1')
hist(f[,2], col='gold', breaks=20, main='histogram of z2', xlab='z2')
hist(f[,3], col='gold', breaks=20, main='histogram of z3', xlab='z3')
# scatter plots
plot(f[,1], col='gold', main='scatter plot of z1')
plot(f[,2], col='gold', main='scatter plot of z2')
plot(f[,3], col='gold', main='scatter plot of z3')
# 3D scatter plots
library("scatterplot3d")
scatterplot3d(f, pch = 16, color="gold")