bistochastic.3x3 <- function() { B <- matrix(0, 3, 3) ## 2 df tmp.1 <- runif(3) B[1,] <- tmp.1/sum(tmp.1) ## 1 df tmp.2 <- runif(2) B[2:3, 1] <- (1-B[1,1]) * tmp.2/sum(tmp.2) ## 1 df B[2, 2] <- runif(1, max=min(1-B[1,2], 1-B[2,1])) ## Fill in the rest B[2,3] <- 1-sum(B[2, 1:2]) B[3,2] <- 1-sum(B[1:2, 2]) B[3,3] <- 1-sum(B[1:2, 3]) B }
B <- bistochastic.3x3() apply(B, 1, sum) apply(B, 2, sum) To extend this to larger than 3x3 requires the same kind of conditional generation of alternating rows and columns of the matrix. The hard part is the extension of the two-way conditioning I illustrated in the B[2, 2] line. Rich ______________________________________________ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.