[R] Generate a random bistochastic matrix

[Richard M. Heiberger]

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