Essentially the same idea of Bill Huberwhuber, nowthis time using Dirichlets: the lines of the matrix are independent Dirichlets, with each Dirichlet putting more mass on the corresponding diagonal element.
rdirichlet <- function(a) {
x <- sapply(a, function(a) rgamma(1, a, 1))
return(x / sum(x))
}
n <- 5
A0 <- 10
a <- matrix(data = 1, nrow = n, ncol = n)
diag(a) <- rep(A0, n)
x <- matrix(nrow = n, ncol = n)
for (i in 1:n) x[i,] <- rdirichlet(a[i,])
Reject if x
does not satisfy the constraints. If the rejection rate is too high, increase the value of A0
. Here is a sample x
.
> x
[,1] [,2] [,3] [,4] [,5]
[1,] 0.484437196 0.03511667 0.04919437 0.12718905 0.304062711
[2,] 0.008626386 0.76520244 0.03231747 0.02454215 0.169311552
[3,] 0.176631389 0.01971251 0.55780424 0.07952712 0.166324747
[4,] 0.003109732 0.12056624 0.09335086 0.77330892 0.009664246
[5,] 0.037015097 0.02485376 0.04536731 0.05083834 0.841925490
P.S. VectorizePlease, vectorize that pesky for
.