以不同的概率重复采样答案

【问题标题】：Sampling repeatedly with different probability以不同的概率重复采样
【发布时间】：2013-12-28 19:22:44
【问题描述】：

在下面的代码中，“权重”是一个权重集的大矩阵。该矩阵由假设 1000 行和 4 列组成。每行是一组权重（每行元素之和等于1）。

另外，有四个对象，我想根据每个权重集选择其中一个。换句话说，应该对所有权重集重复这种随机选择。

现在我已经用for 解决了这个问题。但是有没有更有效的方法在 R 中对其进行编码？

y <- c("a", "b", "c", "d")
for(i in 1:nrow(Weight)){
  selection[i] <- sample(y, 1, prob=Weight[i,]) #selection is a vector with the same number of rows as Weight
}

【问题讨论】：

如何定义“高效”？

标签： r

【解决方案1】：

一种更有效的方法是首先计算权重的逐行累积总和，然后在 0 和 1 之间绘制一个数字，然后查看它在该累积总和中的位置。这样，您只需一次调用 runif 即可获取随机数据，而不是使用其他方法调用 1000。

Weight <- matrix(sample(1:100, 1000 * 4, TRUE), 1000, 4)

x <- runif(nrow(Weight))
cumul.w <- Weight %*% upper.tri(diag(ncol(Weight)), diag = TRUE) / rowSums(Weight)
i <- rowSums(x > cumul.w) + 1L
selection <- y[i]

还要注意我是如何通过乘以三角矩阵而不是使用较慢的apply(Weight, 1, cumsum) 来计算累积和的。一切都是矢量化的，所以它应该比使用 apply 或 for 循环更快。

与apply 和for 的基准比较：

f_runif <- function(Weight, y) {
  x <- runif(nrow(Weight))
  cumul.w <- Weight %*% upper.tri(diag(ncol(Weight)), diag = TRUE) /
    rowSums(Weight)
  i <- rowSums(x > cumul.w) + 1L
  y[i]
}

f_for <- function(Weight, y) {
  selection <- rep(NA, nrow(Weight))
  for(i in 1:nrow(Weight)){
    selection[i] <- sample(y, 1, prob=Weight[i,])
  }
}

f_apply <- function(Weight, y) {
  apply(Weight, 1, function(w)sample(y, 1, prob=w))
}

y <- c("a", "b", "c", "d")
Weight <- matrix(sample(1:100, 1000 * 4, TRUE), 1000, 4)

library(microbenchmark)
microbenchmark(f_runif(Weight, y),
               f_for  (Weight, y),
               f_apply(Weight, y))

# Unit: microseconds
#                expr       min        lq    median         uq       max neval
#  f_runif(Weight, y)   223.635   231.111   274.531   281.2165  1443.208   100
#    f_for(Weight, y) 10220.674 11238.660 11574.039 11917.1610 14583.028   100
#  f_apply(Weight, y)  9006.974 10016.747 10509.150 10879.9245 27060.189   100

【讨论】：

【解决方案2】：

将您的 sample 包装到一个函数中，该函数只允许您传递一个参数，即来自 Weight 的一行：

myfun <- function(w) {
    sample(y, 1, prob=w)
}

那么你可以使用其中一个apply族：

apply(Weight, 1, myfun)

但是，只要您预先分配了selection，您的方法就不会非常低效。

【讨论】：

或者，你可以直接使用apply而不创建额外的函数：apply(Weight, 1, sample, x=y, size=1, replace=FALSE)