这是社区生态学中的一个经典（也是困难的）问题。picante软件包有一个基于Miklos & Podani（2004）的蒙特卡罗算法，该算法在保留边缘的情况下随机化二进制矩阵。然而，该算法假设您从保留约束的二进制矩阵开始;然后，它会提供尽可能多的具有这些约束的随机矩阵。
我想不出一个简单的算法来生成一个二进制矩阵（即使是非随机的），并满足作为起始值的约束条件，所以我使用了蛮力-我使用r2dtable()生成了很多的矩阵，并认为其中一些将是二进制矩阵。
我们走吧

r2 dtable生成起始矩阵

rt <- c(4, 2, 3, 5, 3)
ct <- c(5, 1, 5, 2, 4)
set.seed(101)
system.time(tt <- r2dtable(100000, rt, ct))  ## 0.06 seconds
w <- which(sapply(tt, max) == 1)
length(w) ## 36/100000 (0.036%) are binary
m0 <- tt[[w[1]]]  ## pick the first one

如果你不关心效率，你可以停在那里。如果你需要成千上万的矩阵满足条件，但是，第二阶段是更好的...

试换 Shuffle

Miklos和Podani的算法进行“试验交换”（交换保留行/列总数的行/列对）;默认情况下，它进行1000次交换以随机化矩阵。

library(picante)
results <- list(m0)
nm <- 100000
pb <- txtProgressBar(max = nm, style = 3)
system.time(
  for (i in 1:nm) {
   setTxtProgressBar(pb, i)
   results[[i+1]] <- randomizeMatrix(results[[i]], null.model = "trialswap")
  }
)

在我的机器上，这在7秒内生成了10^5个矩阵（根据我的计算，通过r2dtable()，这将花费大约20倍的时间）。

然而……

以这种方式生成的所有10^5矩阵，* 和 * 满足约束的所有r2dtable矩阵（只有36个），是相同的！（我们通过将整个矩阵折叠成二进制字符串来进行比较......）可能只有一个矩阵满足这些约束，或者在一个很大的空间里有一个很小的数字，所以很难从一个到另一个。

table(sapply(results, paste, collapse = ""))
## 1111100010111111001010111 
##                    100001 
table(sapply(tt[w], paste, collapse = ""))
## 1111100010111111001010111 
##                       36

最后的想法

• 如果你已经有了一个满足约束的矩阵（可能是非随机的），或者你有一个合理的算法来生成一个示例，你就可以开始了
• 试用交换算法应该可以很好地扩展（扩展到更大的大小）;典型的问题可能是30 x30左右。然而，r2dtable()黑客，如果你的真实的问题比5x 5大得多，那就不切实际了。（模拟退火或遗传算法可能会找到一个初始矩阵-从满足行约束的随机矩阵开始，然后用一个等于列约束的平方偏差的目标函数进行 Shuffle -但我没有花时间去弄清楚。）

Miklos I. & Podani J.2004.Randomization of presence-absence matrix：注解和新算法。生态学85：86-92.