[R] Re : Generate a random bistochastic matrix

[Florent Bresson]

Thanks, I tried someting like this, but computation takes times for large matrices

btransf    <-    function(y,X=length(y)^4)    { 
        N    <-    length(y)
        bm    <-    matrix(rep(1/N,N^2),N,N)
        for(j in 1:X){
                    coord    <-    sample(1:N,4,replace=T)
                    d    <-    runif(1,0,min(bm[coord[1],coord[2]],bm[coord[3],coord[4]]))
                    bm[coord[1],coord[2]]    <-    bm[coord[1],coord[2]]-d 
                    bm[coord[3],coord[4]]    <-    bm[coord[3],coord[4]]-d
                    bm[coord[1],coord[4]]    <-    bm[coord[1],coord[4]]+d 
                    bm[coord[2],coord[3]]    <-    bm[coord[2],coord[3]]+d    }    
                    y.btransf    <-    bm%*%y
                    y.btransf    <-    y.btransf+(mean(y)-mean(y.btransf))
                    as.vector(y.btransf)    }

the fonction is designed to perform a mean-preserving transformation of a vector.

----- Message d'origine ----
De : Richard M. Heiberger <rmh at temple.edu>
À : Florent Bresson <f_bresson at yahoo.fr>; r-help at stat.math.ethz.ch
Envoyé le : Lundi, 16 Octobre 2006, 14h58mn 13s
Objet : Re: [R] Generate a random bistochastic matrix

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