R 从不允许相邻元素的向量中采样

Question

假设我被允许将 100% 的权重分布在一个 5 长度的向量上。但是，我不能将权重放入两个相邻的值中，并且任何值都不能超过 50%。

例如，

[0, .5, 0, 0, .5] is good
[.5, .5, 0, 0,0] is not good
[.2, 0, .2, 0, .6] is good
[.2, 0, .2, .2, .2] is not good

我想生成 10,000 个这样的向量来运行蒙特卡罗模拟。

我想我可以用 expand.grid 做到这一点，但我不太确定如何。

我可以随机生成一个然后：

nonzero_weights = which(starting_weights>0)
grid_positions = expand.grid(startingPos = nonzero_weights, endingPos = nonzero_weights)

然后进行一些过滤和删除，但这看起来很乱。如果我不需要它们，为什么要生成。有没有更清洁的方法来做到这一点？

Answer 1

如果我们没有邻接限制，使用R 中当前可用的工具，这个问题不会那么困难（请参阅this answer 更多信息）。有了邻接限制，我们必须做更多的工作才能得到我们想要的结果。

我们首先注意到，由于在具有 n 列的向量的一行中不能有 2 个连续的数字（OP 在 cmets 中澄清了他们需要 n = 11所以我们将使用它作为我们的测试用例），具有值的最大列数等于11 - floor(11 / 2) = 6。当值存在于列1 3 5 7 9 11 中时会发生这种情况。我们还应该注意，由于最大值的上限为 0.5，并且我们需要行总和为 1，因此具有值的最小列数等于 2，因为 ceiling(1 / 0.5) = 2。有了这些信息，我们就可以开始攻击了。

我们首先生成 11 个选择 2 到 6 的每个组合。然后我们筛选出违反邻接限制的组合。后一部分可以通过获取每一行的diff 并检查是否有任何结果差异等于 1 来轻松实现。观察（注意，我们使用RcppAlgos（我是作者）进行所有计算）：

library(RcppAlgos)

vecLen <- 11L
lowComb <- as.integer(ceiling(1 / 0.5))
highComb <- 6L
numCombs <- length(lowComb:highComb)

allCombs <- lapply(lowComb:highComb, function(x) {
    comboGeneral(vecLen, x)
})

validCombs <- lapply(allCombs, function(x) {
    which(apply(x, 1, function(y) {
        !any(diff(y) == 1L)
    }))
})

combLen <- lengths(validCombs)
combLen
[1] 45 84 70 21  1

## subset each matrix of combinations using the
## vector of validCombs obtained above
myCombs <- lapply(seq_along(allCombs), function(x) {
    allCombs[[x]][validCombs[[x]], ]
})

我们现在需要找到所有seq(0.05, 0.5, 0.05) 的组合，对于上面计算的每个可能的长度，它们的总和为 1。使用comboGeneral的约束功能，这很容易：

combSumOne <- lapply(lowComb:highComb, function(x) {
    comboGeneral(seq(5L,50L,5L), x, TRUE, 
                 constraintFun = "sum", 
                 comparisonFun = "==", 
                 limitConstraints = 100L) / 100
})

groupLen <- sapply(combSumOne, nrow)
groupLen
1 13 41 66 78

现在，我们使用上面的myCombs 创建一个包含所需列数的矩阵并填充所有可能的组合，以确保满足邻接要求。

myCombMat <- matrix(0L, nrow = sum(groupLen * combLen), ncol = vecLen)
s <- g <- 1L
e <- combRow <- nrow(combSumOne[[1L]])

for (a in myCombs[-numCombs]) {
    for (i in 1:nrow(a)) {
        myCombMat[s:e, a[i, ]] <- combSumOne[[g]]
        s <- e + 1L
        e <- e + combRow
    }
    e <- e - combRow
    g <- g + 1L
    combRow <- nrow(combSumOne[[g]])
    e <- e + combRow
}

## the last element in myCombs is simply a
## vector, thus nrow would return NULL
myCombMat[s:e, myCombs[[numCombs]]] <- combSumOne[[g]]

这里是输出的一瞥：

head(myCombMat)
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
[1,]  0.5    0  0.5  0.0  0.0  0.0  0.0  0.0    0     0     0
[2,]  0.5    0  0.0  0.5  0.0  0.0  0.0  0.0    0     0     0
[3,]  0.5    0  0.0  0.0  0.5  0.0  0.0  0.0    0     0     0
[4,]  0.5    0  0.0  0.0  0.0  0.5  0.0  0.0    0     0     0
[5,]  0.5    0  0.0  0.0  0.0  0.0  0.5  0.0    0     0     0
[6,]  0.5    0  0.0  0.0  0.0  0.0  0.0  0.5    0     0     0

tail(myCombMat)
        [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
[5466,] 0.10    0 0.10    0 0.20    0 0.20    0 0.20     0  0.20
[5467,] 0.10    0 0.15    0 0.15    0 0.15    0 0.15     0  0.30
[5468,] 0.10    0 0.15    0 0.15    0 0.15    0 0.20     0  0.25
[5469,] 0.10    0 0.15    0 0.15    0 0.20    0 0.20     0  0.20
[5470,] 0.15    0 0.15    0 0.15    0 0.15    0 0.15     0  0.25
[5471,] 0.15    0 0.15    0 0.15    0 0.15    0 0.20     0  0.20

set.seed(42)
mySamp <- sample(nrow(myCombMat), 10)
sampMat <- myCombMat[mySamp, ]
rownames(sampMat) <- mySamp

sampMat
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
5005 0.00 0.05 0.00 0.05 0.00 0.15 0.00 0.35 0.00   0.4  0.00
5126 0.00 0.15 0.00 0.15 0.00 0.20 0.00 0.20 0.00   0.0  0.30
1565 0.10 0.00 0.15 0.00 0.00 0.00 0.25 0.00 0.00   0.5  0.00
4541 0.05 0.00 0.05 0.00 0.00 0.15 0.00 0.00 0.25   0.0  0.50
3509 0.00 0.00 0.15 0.00 0.25 0.00 0.25 0.00 0.00   0.0  0.35
2838 0.00 0.10 0.00 0.15 0.00 0.00 0.35 0.00 0.00   0.0  0.40
4026 0.05 0.00 0.10 0.00 0.15 0.00 0.20 0.00 0.50   0.0  0.00
736  0.00 0.00 0.10 0.00 0.40 0.00 0.00 0.00 0.00   0.0  0.50
3590 0.00 0.00 0.15 0.00 0.20 0.00 0.00 0.30 0.00   0.0  0.35
3852 0.00 0.00 0.00 0.05 0.00 0.20 0.00 0.30 0.00   0.0  0.45

all(rowSums(myCombMat) == 1)
[1] TRUE

如您所见，每一行的总和为 1，并且没有相邻的值。

如果您真的想要排列，我们可以生成所有可能的长度为 1 的 seq(0.05, 0.5, 0.05) 排列（就像我们为组合所做的那样）：

permSumOne <- lapply(lowComb:highComb, function(x) {
    permuteGeneral(seq(5L,50L,5L), x, TRUE, 
                   constraintFun = "sum", 
                   comparisonFun = "==", 
                   limitConstraints = 100L) / 100
})

groupLenPerm <- sapply(permSumOne, nrow)
groupLenPerm
[1]     1    63   633  3246 10872

并使用这些来创建我们所有可能的排列矩阵，总和为 1 并满足我们的邻接要求：

myPermMat <- matrix(0L, nrow = sum(groupLenPerm * combLen), ncol = vecLen)
s <- g <- 1L
e <- permRow <- nrow(permSumOne[[1L]])

for (a in myCombs[-numCombs]) {
    for (i in 1:nrow(a)) {
        myPermMat[s:e, a[i, ]] <- permSumOne[[g]]
        s <- e + 1L
        e <- e + permRow
    }
    e <- e - permRow
    g <- g + 1L
    permRow <- nrow(permSumOne[[g]])
    e <- e + permRow
}

## the last element in myCombs is simply a
## vector, thus nrow would return NULL
myPermMat[s:e, myCombs[[numCombs]]] <- permSumOne[[g]]

再一次，这里是输出的一瞥：

head(myPermMat)
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
[1,]  0.5    0  0.5  0.0  0.0  0.0  0.0  0.0    0     0     0
[2,]  0.5    0  0.0  0.5  0.0  0.0  0.0  0.0    0     0     0
[3,]  0.5    0  0.0  0.0  0.5  0.0  0.0  0.0    0     0     0
[4,]  0.5    0  0.0  0.0  0.0  0.5  0.0  0.0    0     0     0
[5,]  0.5    0  0.0  0.0  0.0  0.0  0.5  0.0    0     0     0
[6,]  0.5    0  0.0  0.0  0.0  0.0  0.0  0.5    0     0     0

tail(myPermMat)
          [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
[128680,] 0.15    0 0.20    0 0.20    0 0.15    0 0.15     0  0.15
[128681,] 0.20    0 0.15    0 0.15    0 0.15    0 0.15     0  0.20
[128682,] 0.20    0 0.15    0 0.15    0 0.15    0 0.20     0  0.15
[128683,] 0.20    0 0.15    0 0.15    0 0.20    0 0.15     0  0.15
[128684,] 0.20    0 0.15    0 0.20    0 0.15    0 0.15     0  0.15
[128685,] 0.20    0 0.20    0 0.15    0 0.15    0 0.15     0  0.15

all(rowSums(myPermMat) == 1)
[1] TRUE

而且，正如 OP 所述，如果我们想随机选择其中的 10000 个，我们可以使用 sample 来实现：

set.seed(101)
mySamp10000 <- sample(nrow(myPermMat), 10000)
myMat10000 <- myPermMat[mySamp10000, ]
rownames(myMat10000) <- mySamp10000

head(myMat10000)
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
47897 0.00  0.0 0.00 0.50  0.0 0.25  0.0 0.00 0.05   0.0  0.20
5640  0.25  0.0 0.15 0.00  0.1 0.00  0.5 0.00 0.00   0.0  0.00
91325 0.10  0.0 0.00 0.15  0.0 0.40  0.0 0.00 0.20   0.0  0.15
84633 0.15  0.0 0.00 0.35  0.0 0.30  0.0 0.10 0.00   0.1  0.00
32152 0.00  0.4 0.00 0.05  0.0 0.00  0.0 0.25 0.00   0.3  0.00
38612 0.00  0.4 0.00 0.00  0.0 0.35  0.0 0.10 0.00   0.0  0.15

由于RcppAlgos 的效率很高，上述所有步骤都会立即返回。在我的 2008 Windows 机器 i5 2.5 GHz 上，整个生成（包括排列）需要不到 0.04 秒。

Answer 2

首先，您可以通过从上一个样本中删除采样索引来生成二进制样本。然后生成要分配给这些二进制样本的权重：

idx <- 1:11

system.time(
    binsampl <- t(replicate(10000L, {
        x <- rep(0L, length(idx))
        while(length(idx) > 0L) {
            chosen <- if (length(idx) > 1L) sample(idx, 1L) else idx
            idx <- setdiff(idx, chosen + -1L:1L)
            x[chosen] <- 1L
        }
        x
    }))
)

system.time(
    weights <- t(apply(binsampl, 1, function(s) {
        y <- runif(sum(s))
        s[s==1L] <- y/sum(y) 
        s
    }))
)
head(weights)

输出：

            [,1]       [,2]      [,3]      [,4]       [,5]      [,6]       [,7]      [,8]       [,9]
[1,] 0.114636912 0.00000000 0.1136963 0.0000000 0.00000000 0.1938791 0.00000000 0.3495739 0.00000000
[2,] 0.267907091 0.00000000 0.1487623 0.0000000 0.21628596 0.0000000 0.08326985 0.0000000 0.03803797
[3,] 0.000000000 0.06195168 0.0000000 0.0000000 0.07972502 0.0000000 0.00000000 0.3749550 0.00000000
[4,] 0.083384611 0.00000000 0.0000000 0.3867607 0.00000000 0.0000000 0.16300188 0.0000000 0.00000000
[5,] 0.005233208 0.00000000 0.4106275 0.0000000 0.15796746 0.0000000 0.10168549 0.0000000 0.00000000
[6,] 0.188153707 0.00000000 0.1867017 0.0000000 0.29426748 0.0000000 0.00000000 0.2962538 0.00000000
         [,10]     [,11]
[1,] 0.2282138 0.0000000
[2,] 0.0000000 0.2457368
[3,] 0.0000000 0.4833683
[4,] 0.3668528 0.0000000
[5,] 0.3244863 0.0000000
[6,] 0.0000000 0.0346233

在我使用 R-3.5.1 Windows x64 8GB RAM 2.8GHz 处理器的机器上生成 10k 个样本所需的时间不到 1 秒。