代码之家 › 专栏 › 技术社区 › user1357015

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

montecarlo permutation combinations r

1

user1357015 · 技术社区 · 6 年前

假设允许我沿着一个5长度的向量分布100%的权重。但是,我不能将权重放入两个相邻的值中,并且任何值都不能超过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

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

我想我可以用它 expand.grid 但我不太确定怎么做。

我可以随机生成一个,然后:

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

然后做一些过滤和删除,但这看起来很混乱。如果我不需要它们,为什么要生成它们。有更干净的方法吗?

2 回复 | 直到 6 年前

1

Joseph Wood 6 年前

如果我们没有邻接限制,这个问题对于当前可用的工具来说就不会那么困难了 R (见 this answer 更多信息)。有了邻接限制,我们必须做更多的工作才能得到我们想要的结果。

我们首先注意到,因为在一个向量的行中不能有两个连续的数字 n个列(操作在它们需要的注释中澄清 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,并且没有相邻的值。

如果你真的想要置换,我们可以生成 顺序(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公司 高效,以上所有步骤立即返回。在我的2008 Windows机器i5 2.5 GHz上,整个一代(包括排列)所需的时间不到0.04秒。

2

1

chinsoon12 6 年前

首先,可以通过从prev samples中移除采样索引来生成二进制样本。然后生成要分配给这些二进制样本的权重:

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