如何在R语言中生成满足行列总和约束的随机0-1存在/缺失矩阵？

你遇到的问题是r2dtable()生成的矩阵允许元素大于1，这是因为它的设计目标是生成非负整数构成的列联表，而你需要的是严格的0/1二元矩阵（也叫二元列联表），每个元素只能表示存在/缺失，不能重复计数。下面提供两种可靠的解决方案：

方法一：使用mipfp包生成符合约束的0-1矩阵 #

mipfp包专门用于处理边际约束下的矩阵抽样，支持生成0-1类型的矩阵：

# 安装并加载包
install.packages("mipfp")
library(mipfp)

# 定义你的行列和
row_sums <- c(4, 2, 3, 5, 3)
col_sums <- c(5, 1, 5, 2, 4)
n_rows <- length(row_sums)
n_cols <- length(col_sums)

# 生成初始满足行和的0-1矩阵作为种子
init_matrix <- t(sapply(row_sums, function(r) {
  x <- rep(0, n_cols)
  x[sample(n_cols, r)] <- 1
  x
}))

# 指定边际约束：第1维度（行）的和为row_sums，第2维度（列）的和为col_sums
targets <- list(row_sums, col_sums)
dim_targets <- list(1, 2)

# 抽样生成符合条件的矩阵
set.seed(123) # 设置随机种子保证结果可复现
result <- genpfp(seed = init_matrix, targets = targets, dim.targets = dim_targets, 
                 method = "ipfp", iter = 1000, verbose = FALSE)

# 提取最终的0-1矩阵
binary_matrix <- result$x.hat

# 验证约束是否满足
rowSums(binary_matrix) # 应等于row_sums
colSums(binary_matrix) # 应等于col_sums
print(binary_matrix)

方法二：手动实现MCMC交换法抽样 #

如果不想依赖第三方包，可以用马尔可夫链蒙特卡洛（MCMC）的交换法，通过保持行列和不变的元素交换来生成随机0-1矩阵：

set.seed(123)
row_sums <- c(4, 2, 3, 5, 3)
col_sums <- c(5, 1, 5, 2, 4)
n_rows <- length(row_sums)
n_cols <- length(col_sums)

# 生成初始满足行和的0-1矩阵
init_mat <- matrix(0, n_rows, n_cols)
for (i in 1:n_rows) {
  init_mat[i, sample(n_cols, row_sums[i])] <- 1
}

# 定义MCMC交换函数
mcmc_binary_table <- function(mat, n_iter = 10000) {
  current_mat <- mat
  for (iter in 1:n_iter) {
    # 随机选取不同的两行和两列
    i <- sample(n_rows, 1)
    j <- sample(setdiff(1:n_rows, i), 1)
    k <- sample(n_cols, 1)
    l <- sample(setdiff(1:n_cols, k), 1)
    
    # 检查四个位置的元素是否可以交换（保证交换后行列和不变）
    vals <- current_mat[c(i,j), c(k,l)]
    if (all(vals %in% c(0,1)) && sum(vals) == 2 && vals[1,1] != vals[1,2]) {
      current_mat[c(i,j), c(k,l)] <- 1 - vals
    }
  }
  return(current_mat)
}

# 运行MCMC抽样
result_mat <- mcmc_binary_table(init_mat, n_iter = 10000)

# 验证约束
rowSums(result_mat)
colSums(result_mat)
print(result_mat)

注意事项 #

• 首先要确保行和的总和等于列和的总和，否则不存在符合条件的矩阵（你的例子中两者都是17，是可行的）。
• MCMC方法需要足够多的迭代次数，才能保证抽样结果的随机性。

内容的提问来源于stack exchange，提问作者mycorrhiza