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 个这样的向量，从中运行进行 monte carlo 模拟。

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

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

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

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

Answer 1

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

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 秒。

Answer 2

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

我们首先注意到，因为我们不能在具有 n 列的向量的一行中有 2 个连续的数字（OP 在评论中澄清他们需要 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，可以轻松实现后一部分。观察（N.B。我们使用 RcppAlgos（我是author) 用于所有计算):

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，并且没有相邻值。

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

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效率很高，return以上的所有步骤都立即完成。在我的 2008 Windows 机器 i5 2.5 GHz 上，整个生成过程（包括排列）不到 0.04 秒。

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

R sampling from a vector where adjacent elements are not allowed

r

combinations

permutation

montecarlo