在 0/1 矩阵内采样

Question

我有一个由零和一组成的 $N$x$N$ 矩阵 $G$。我只想从伯努利分布中重新采样 $G$ 中的那些，其中概率 $p$ 随原始变化。这个想法是用独立伯努利变量的实现替换 $G$ 中的所有值。我在 R 中试过这个：

for(i in 1:N){ 
for(j in 1:N){ 
G[i,j] <- ifelse(G[i,j] == 1,rbinom(1,1,p[i]),0)}}

但是，为了计算时间，最好避免使用for循环。谁知道如何以矩阵形式写这个？非常感谢！

Answer 1

利用R 的向量循环及其面向列的存储。我会在展示一种解决方案后进行解释。

#
# Generate a 0/1 matrix `G` for testing.
#
N <- 1e3
G <- matrix(rbinom(N^2, 1, 1/2), N)
#
# Generate the vector of probabilities for testing.
#
p <- runif(N)
#
# Here is a solution.
#
U <- runif(N*N) <= p # Random 0/1 values
G.0 <- G * U         # Conditionally replace elements of G

runif 函数为 G 中的每个条目生成一个统一变量。将它与向量 p 进行比较会创建 N 个 0/1 值的向量（1 表示真，0 表示假），所有向量都串在一起，在位置 i、i+N、i+2N、 ..., i+(N-1)N 由 p[i] 的值给出，i=1, 2, ..., N。当它在最后一行乘以 G 时，当旧值为零或新值为零时，新值为零，这是期望的结果，因为 G 的值在内部也按列排列。

这个解决方案比原来的双 for 循环快大约 20 倍。它避免了条件执行 (ifelse) 的缓慢，并利用了 R 基本操作（runif、<= 和 * 的任何内部并行性和优化).您确实付出了（相对较小的）代价：需要 RAM 来存储这些随机结果（此处明确放入数组 U）。

Answer 2

If using Python is not off the table, use Numpy:

import numpy as np
G = np.array([[0,1,0,0], [0,0,0,1], [0,1,1,0], [1,1,1,1]])  
# This constructs G row-by-row 
probs = np.array([0.1, 0.2, 0.3, 0.4])  # Replace with your row-wise probabilities
samples = np.random.binomial(n=1, p=probs, size=(4,4))  # Binomial with n=1 is Bernoulli
desired_matrix = np.multiply(G, samples)

np.multiply 函数按元素乘以数组，因此零条目保持为零，而 1 条目乘以 0 或 1，具体取决于 np.random.binomial 返回的实现。一个关键点是我相信 Numpy 足够聪明，可以在指定 size 时将 probs 中的每个概率广播到右行，但我不是肯定是因为文档没有提到它。

在 0/1 矩阵内采样

Sampling within a 0/1 matrix

r

matrix

sampling