在 0/1 矩阵内采样答案

【问题标题】：Sampling within a 0/1 matrix在 0/1 矩阵内采样
【发布时间】：2017-10-24 17:23:32
【问题描述】：

我有一个由 0 和 1 组成的 $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 循环。有谁知道如何以矩阵形式写这个？非常感谢你！

【问题讨论】：

我不关注——如果你重新采样“只有那些”，那么你为什么不应该只得到一组呢？你能澄清一下你想要做什么吗？在什么意义上，这段代码还不是“矩阵形式”？您是否想说 (1) 您想用来自独立伯努利变量的实现替换 $G$ 中的所有 1 以及 (2) 您想避免使用 for 循环？
@whuber 抱歉这个模糊的问题。是的（1）和（2）。我已经编辑了帖子。

标签： r matrix sampling

【解决方案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 中）。

【讨论】：

感谢您回答我的问题 :)

【解决方案2】：

如果不使用 Python，请使用 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 中的每个概率广播到右行，但我不确定，因为文档没有提一下。

【讨论】：