考虑向量s如下：
s=seq(0.01, 0.99, 0.01)
> s
[1] 0.01 0.02 0.03 0.04 0.05 0.06 0.07
0.08 0.09 .......... 0.89 0.90 0.91 0.92 0.93 0.94 0.95 0.96 0.97 0.98 0.99
现在给定s一个固定的长度m，我想为所有可能的长度排列提供一个矩阵，m使得矩阵的每一行总和为1（不包括蛮力方法）。
例如，如果m=4（即列数），所需的矩阵将是这样的：
0.01 0.01 0.01 0.97
0.02 0.01 0.01 0.96
0.03 0.01 0.01 0.95
0.04 0.01 0.01 0.94
0.05 0.01 0.01 0.93
0.06 0.01 0.01 0.92
.
.
.
0.53 0.12 0.30 0.05
.
.
.
0.96 0.02 0.01 0.01
0.97 0.01 0.01 0.01
.
.
.
0.01 0.97 0.01 0.01
.
.
.
以下是使用递归的方法：
permsum <- function(s,m) if (m==1L) matrix(s) else do.call(rbind,lapply(seq_len(s-m+1L),function(x) unname(cbind(x,permsum(s-x,m-1L)))));
res <- permsum(100L,4L);
head(res);
## [,1] [,2] [,3] [,4]
## [1,] 1 1 1 97
## [2,] 1 1 2 96
## [3,] 1 1 3 95
## [4,] 1 1 4 94
## [5,] 1 1 5 93
## [6,] 1 1 6 92
tail(res);
## [,1] [,2] [,3] [,4]
## [156844,] 95 2 2 1
## [156845,] 95 3 1 1
## [156846,] 96 1 1 2
## [156847,] 96 1 2 1
## [156848,] 96 2 1 1
## [156849,] 97 1 1 1
您可以除以 100 得到分数，而不是整数：
head(res)/100;
## [,1] [,2] [,3] [,4]
## [1,] 0.01 0.01 0.01 0.97
## [2,] 0.01 0.01 0.02 0.96
## [3,] 0.01 0.01 0.03 0.95
## [4,] 0.01 0.01 0.04 0.94
## [5,] 0.01 0.01 0.05 0.93
## [6,] 0.01 0.01 0.06 0.92
首先让我们定义输入：
s这是输出矩阵中每一行应求和的目标值。m这是要在输出矩阵中生成的列数。与浮点算术相比，使用整数算术计算结果更有效、更可靠，因此我将解决方案设计为仅使用整数。因此s是一个标量整数，表示目标整数和。
seq_len()现在让我们检查由非基本情况生成的序列：
seq_len(s-m+1L)
这会生成一个从 1 到最高可能值的序列，该序列可能是剩余列的总和s的一部分。m例如，考虑以下情况s=100,m=4：我们可以使用的最大数字是 97，参与的总和为 97+1+1+1。剩余的每一列将可能的最大值减 1，这就是为什么我们在计算序列长度时必须减去m的原因。s
生成序列的每个元素都应被视为求和中加数的一种可能“选择”。
do.call(rbind,lapply(seq_len(s-m+1L),function(x) ...))
对于每个选择，我们必须递归。我们可以用它lapply()来做到这一点。
为了向前跳转，lambda 将进行一次递归调用permsum()，然后cbind()返回当前选择的值。这将产生一个矩阵，对于这个级别的递归，它总是具有相同的宽度。因此，该lapply()调用将返回一个矩阵列表，所有矩阵的宽度都相同。然后我们必须将它们行绑定在一起，这就是为什么我们必须在do.call(rbind,...)这里使用这个技巧。
unname(cbind(x,permsum(s-x,m-1L)))
lambda 的主体相当简单；我们cbind()将当前选择x与递归调用的返回值结合起来，完成此子矩阵的求和。不幸的是，我们必须调用unname()，否则最终从参数设置的每一列都x将具有列名x。
这里最重要的细节是递归调用的参数选择。首先，由于 lambda 参数x刚刚在当前递归评估期间被选中，我们必须将其减去s以获得新的求和目标，即将到来的递归调用将负责实现该目标。因此第一个参数变为s-x。其次，因为 的选择x占用一列，我们必须从 中减去 1 m，这样递归调用将在其输出矩阵中产生少一列。
if (m==1L) matrix(s) else ...
最后，让我们检查一下基本情况。在递归函数的每次求值中，我们必须检查是否m已达到 1，在这种情况下，我们可以简单地返回所需的总和s本身。
我调查了我的结果和 psidom 的结果之间的差异。例如：
library(data.table);
bgoldst <- function(s,m) permsum(s,m)/s;
psidom <- function(ss,m) { raw <- do.call(data.table::CJ,rep(list(ss),m)); raw[rowSums(raw)==1,]; };
## helper function to sort a matrix by columns
smp <- function(m) m[do.call(order,as.data.frame(m)),];
s <- 100L; m <- 3L; ss <- seq_len(s-1L)/s;
x <- smp(bgoldst(s,m));
y <- smp(unname(as.matrix(psidom(ss,m))));
nrow(x);
## [1] 4851
nrow(y);
## [1] 4809
所以我们的两个结果之间存在 42 行的差异。我决定尝试用以下代码行准确找出哪些排列被省略了。基本上，它比较两个矩阵的每个元素并将比较结果打印为逻辑矩阵。我们可以向下扫描回滚以找到第一个不同的行。以下是摘录的输出：
x==do.call(rbind,c(list(y),rep(list(NA),nrow(x)-nrow(y))));
## [,1] [,2] [,3]
## [1,] TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE
## [3,] TRUE TRUE TRUE
## [4,] TRUE TRUE TRUE
## [5,] TRUE TRUE TRUE
##
## ... snip ...
##
## [24,] TRUE TRUE TRUE
## [25,] TRUE TRUE TRUE
## [26,] TRUE TRUE TRUE
## [27,] TRUE TRUE TRUE
## [28,] TRUE TRUE TRUE
## [29,] TRUE FALSE FALSE
## [30,] TRUE FALSE FALSE
## [31,] TRUE FALSE FALSE
## [32,] TRUE FALSE FALSE
## [33,] TRUE FALSE FALSE
##
## ... snip ...
所以它在第 29 行，我们有第一个差异。这是每个排列矩阵中该行周围的一个窗口：
win <- 27:31;
x[win,]; y[win,];
## [,1] [,2] [,3]
## [1,] 0.01 0.27 0.72
## [2,] 0.01 0.28 0.71
## [3,] 0.01 0.29 0.70 (missing from y)
## [4,] 0.01 0.30 0.69 (missing from y)
## [5,] 0.01 0.31 0.68
## [,1] [,2] [,3]
## [1,] 0.01 0.27 0.72
## [2,] 0.01 0.28 0.71
## [3,] 0.01 0.31 0.68
## [4,] 0.01 0.32 0.67
## [5,] 0.01 0.33 0.66
有趣的是，当您手动计算总和时，缺失的排列通常确实总和为 1。起初我认为是 data.table 的CJ()函数对浮点数做了一些奇怪的事情，但进一步的测试似乎表明它rowSums()正在做一些事情：
0.01+0.29+0.70==1;
## [1] TRUE
ss[1L]+ss[29L]+ss[70L]==1;
## [1] TRUE
rowSums(CJ(0.01,0.29,0.70))==1; ## looks like CJ()'s fault, but wait...
## [1] FALSE
cj <- CJ(0.01,0.29,0.70);
cj$V1+cj$V2+cj$V3==1; ## not CJ()'s fault
## [1] TRUE
rowSums(matrix(c(0.01,0.29,0.70),1L,byrow=T))==1; ## rowSums()'s fault
## [1] FALSE
rowSums()我们可以通过在浮点比较中应用手动（并且有些随意）容差来解决这个怪癖。为此，我们需要获取绝对差，然后与容差进行小于比较：
abs(rowSums(CJ(0.01,0.29,0.70))-1)<1e-10;
## [1] TRUE
因此：
psidom2 <- function(ss,m) { raw <- do.call(data.table::CJ,rep(list(ss),m)); raw[abs(rowSums(raw)-1)<1e-10,]; };
y <- smp(unname(as.matrix(psidom2(ss,m))));
nrow(y);
## [1] 4851
identical(x,y);
## [1] TRUE
感谢 Joseph Wood 指出这确实是排列。我最初命名了我的函数combsum()，但我将其重命名permsum()以反映这一启示。而且，正如约瑟夫所建议的，可以修改算法以产生组合，可以按如下方式完成，现在名副其实combsum()：
combsum <- function(s,m,l=s) if (m==1L) matrix(s) else do.call(rbind,lapply(seq((s+m-1L)%/%m,min(l,s-m+1L)),function(x) unname(cbind(x,combsum(s-x,m-1L,x)))));
res <- combsum(100L,4L);
head(res);
## [,1] [,2] [,3] [,4]
## [1,] 25 25 25 25
## [2,] 26 25 25 24
## [3,] 26 26 24 24
## [4,] 26 26 25 23
## [5,] 26 26 26 22
## [6,] 27 25 24 24
tail(res);
## [,1] [,2] [,3] [,4]
## [7148,] 94 3 2 1
## [7149,] 94 4 1 1
## [7150,] 95 2 2 1
## [7151,] 95 3 1 1
## [7152,] 96 2 1 1
## [7153,] 97 1 1 1
这需要 3 次更改。
首先，我添加了一个新参数l，它代表“限制”。基本上，为了保证每个递归生成一个唯一的组合，我强制每个选择必须小于或等于当前组合中的任何先前选择。这需要将当前上限作为参数l。在顶层调用l可以默认为s，对于 的情况，这实际上太高了m>1，但这不是问题，因为它只是将在序列生成期间应用的两个上限之一。
第二个更改当然是在lambda中进行递归调用时将最新选择x作为参数传递给。llapply()
最后的变化是最棘手的。现在必须按如下方式计算选择序列：
seq((s+m-1L)%/%m,min(l,s-m+1L))
下限必须从使用的 1 提高permsum()到仍然允许降幅组合的最低可能选择。当然，可能的最低选择取决于尚未生产多少列；列越多，我们留给未来选择的“空间”就越多。公式是对son进行整数除法m，但我们还必须有效地“四舍五入”，这就是我m-1L在除法之前添加的原因。我还考虑过进行浮点除法然后调用as.integer(ceiling(...))，但我认为全整数方法要好得多。
例如，考虑 的情况s=10,m=3。要在剩余 3 列的情况下产生 10 的总和，我们不能选择小于 4 的选择，因为如果不沿着组合上升，我们将没有足够的数量来产生 10。在这种情况下，公式将 12 除以 3 得到 4。
上限可以从 中使用的相同公式计算permsum()，除了我们还必须l使用调用来应用当前限制min()。
我已经使用以下代码验证了对于许多随机测试用例，我的新功能combsum()与 Joseph 的功能相同：IntegerPartitionsOfLength()
## helper function to sort a matrix within each row and then by columns
smc <- function(m) smp(t(apply(m,1L,sort)));
## test loop
for (i in seq_len(1000L)) {
repeat {
s <- sample(1:100,1L);
m <- sample(2:5,1L);
if (s>=m) break;
};
x <- combsum(s,m);
y <- IntegerPartitionsOfLength(s,m);
cat(paste0(s,',',m,'\n'));
if (!identical(smc(x),smc(y))) stop('bad.');
};
常见的自包含测试代码：
library(microbenchmark);
library(data.table);
library(partitions);
library(gtools);
permsum <- function(s,m) if (m==1L) matrix(s) else do.call(rbind,lapply(seq_len(s-m+1L),function(x) unname(cbind(x,permsum(s-x,m-1L)))));
combsum <- function(s,m,l=s) if (m==1L) matrix(s) else do.call(rbind,lapply(seq((s+m-1L)%/%m,min(l,s-m+1L)),function(x) unname(cbind(x,combsum(s-x,m-1L,x)))));
IntegerPartitionsOfLength <- function(n, Lim, combsOnly = TRUE) { a <- 0L:n; k <- 2L; a[2L] <- n; MyParts <- vector("list", length=P(n)); count <- 0L; while (!(k==1L) && k <= Lim + 1L) { x <- a[k-1L]+1L; y <- a[k]-1L; k <- k-1L; while (x<=y && k <= Lim) {a[k] <- x; y <- y-x; k <- k+1L}; a[k] <- x+y; if (k==Lim) { count <- count+1L; MyParts[[count]] <- a[1L:k]; }; }; MyParts <- MyParts[1:count]; if (combsOnly) {do.call(rbind, MyParts)} else {MyParts}; };
GetDecimalReps <- function(s,m) { myPerms <- permutations(m,m); lim <- nrow(myPerms); intParts <- IntegerPartitionsOfLength(s,m,FALSE); do.call(rbind, lapply(intParts, function(x) { unique(t(sapply(1L:lim, function(y) x[myPerms[y, ]]))); })); };
smp <- function(m) m[do.call(order,as.data.frame(m)),];
smc <- function(m) smp(t(apply(m,1L,sort)));
bgoldst.perm <- function(s,m) permsum(s,m)/s;
psidom2 <- function(ss,m) { raw <- do.call(data.table::CJ,rep(list(ss),m)); raw[abs(rowSums(raw)-1)<1e-10,]; };
joseph.perm <- function(s,m) GetDecimalReps(s,m)/s;
bgoldst.comb <- function(s,m) combsum(s,m)/s;
joseph.comb <- function(s,m) IntegerPartitionsOfLength(s,m)/s;
排列
## small scale
s <- 10L; m <- 3L; ss <- seq_len(s-1L)/s;
ex <- smp(bgoldst.perm(s,m));
identical(ex,smp(unname(as.matrix(psidom2(ss,m)))));
## [1] TRUE
identical(ex,smp(joseph.perm(s,m)));
## [1] TRUE
microbenchmark(bgoldst.perm(s,m),psidom2(ss,m),joseph.perm(s,m));
## Unit: microseconds
## expr min lq mean median uq max neval
## bgoldst.perm(s, m) 347.254 389.5920 469.1011 420.383 478.7575 1869.697 100
## psidom2(ss, m) 702.206 830.5015 1007.5111 907.265 1038.3405 2618.089 100
## joseph.perm(s, m) 1225.225 1392.8640 1722.0070 1506.833 1860.0745 4411.234 100
## large scale
s <- 100L; m <- 4L; ss <- seq_len(s-1L)/s;
ex <- smp(bgoldst.perm(s,m));
identical(ex,smp(unname(as.matrix(psidom2(ss,m)))));
## [1] TRUE
identical(ex,smp(joseph.perm(s,m)));
## [1] TRUE
microbenchmark(bgoldst.perm(s,m),psidom2(ss,m),joseph.perm(s,m),times=5L);
## Unit: seconds
## expr min lq mean median uq max neval
## bgoldst.perm(s, m) 1.286856 1.304177 1.426376 1.374411 1.399850 1.766585 5
## psidom2(ss, m) 6.673545 7.046951 7.416161 7.115375 7.629177 8.615757 5
## joseph.perm(s, m) 5.299452 10.499891 13.769363 12.680607 15.107748 25.259117 5
## very large scale
s <- 100L; m <- 5L; ss <- seq_len(s-1L)/s;
ex <- smp(bgoldst.perm(s,m));
identical(ex,smp(unname(as.matrix(psidom2(ss,m)))));
## Error: cannot allocate vector of size 70.9 Gb
identical(ex,smp(joseph.perm(s,m)));
## [1] TRUE
microbenchmark(bgoldst.perm(s,m),joseph.perm(s,m),times=1L);
## Unit: seconds
## expr min lq mean median uq max neval
## bgoldst.perm(s, m) 28.58359 28.58359 28.58359 28.58359 28.58359 28.58359 1
## joseph.perm(s, m) 50.51965 50.51965 50.51965 50.51965 50.51965 50.51965 1
组合
## small-scale
s <- 10L; m <- 3L;
ex <- smc(bgoldst.comb(s,m));
identical(ex,smc(joseph.comb(s,m)));
## [1] TRUE
microbenchmark(bgoldst.comb(s,m),joseph.comb(s,m));
## Unit: microseconds
## expr min lq mean median uq max neval
## bgoldst.comb(s, m) 161.225 179.6145 205.0898 187.3120 199.5005 1310.328 100
## joseph.comb(s, m) 172.344 191.8025 204.5681 197.7895 205.2735 437.489 100
## large-scale
s <- 100L; m <- 4L;
ex <- smc(bgoldst.comb(s,m));
identical(ex,smc(joseph.comb(s,m)));
## [1] TRUE
microbenchmark(bgoldst.comb(s,m),joseph.comb(s,m),times=5L);
## Unit: milliseconds
## expr min lq mean median uq max neval
## bgoldst.comb(s, m) 409.0708 485.9739 556.4792 591.4774 627.419 668.4548 5
## joseph.comb(s, m) 2164.2134 3315.0138 3317.9725 3540.6240 3713.732 3856.2793 5
## very large scale
s <- 100L; m <- 6L;
ex <- smc(bgoldst.comb(s,m));
identical(ex,smc(joseph.comb(s,m)));
## [1] TRUE
microbenchmark(bgoldst.comb(s,m),joseph.comb(s,m),times=1L);
## Unit: seconds
## expr min lq mean median uq max neval
## bgoldst.comb(s, m) 2.498588 2.498588 2.498588 2.498588 2.498588 2.498588 1
## joseph.comb(s, m) 12.344261 12.344261 12.344261 12.344261 12.344261 12.344261 1
m=4例如，内存密集型方法是：
raw <- data.table::CJ(s,s,s,s)
result <- raw[rowSums(raw) == 1, ]
head(result)
V1 V2 V3 V4
1: 0.01 0.01 0.01 0.97
2: 0.01 0.01 0.02 0.96
3: 0.01 0.01 0.03 0.95
4: 0.01 0.01 0.04 0.94
5: 0.01 0.01 0.05 0.93
6: 0.01 0.01 0.06 0.92
这是一个将返回纯的算法combinations（顺序无关紧要）。它基于 Jerome Kelleher (链接) 构建的整数分区算法。
library(partitions)
IntegerPartitionsOfLength <- function(n, Lim, combsOnly = TRUE) {
a <- 0L:n
k <- 2L
a[2L] <- n
MyParts <- vector("list", length=P(n))
count <- 0L
while (!(k==1L) && k <= Lim + 1L) {
x <- a[k-1L]+1L
y <- a[k]-1L
k <- k-1L
while (x<=y && k <= Lim) {a[k] <- x; y <- y-x; k <- k+1L}
a[k] <- x+y
if (k==Lim) {
count <- count+1L
MyParts[[count]] <- a[1L:k]
}
}
MyParts <- MyParts[1:count]
if (combsOnly) {do.call(rbind, MyParts)} else {MyParts}
}
system.time(res <- combsum(100L,5L))
user system elapsed
0.75 0.00 0.77
system.time(a <- IntegerPartitionsOfLength(100, 5))
user system elapsed
1.36 0.37 1.76
identical(smc(a),smc(res))
[1] TRUE
head(a)
[,1] [,2] [,3] [,4] [,5]
[1,] 1 1 1 1 96
[2,] 1 1 1 2 95
[3,] 1 1 1 3 94
[4,] 1 1 1 4 93
[5,] 1 1 1 5 92
[6,] 1 1 1 6 91
一个非常大的例子（注意使用smc@bgoldst 创建的函数）：
system.time(a <- IntegerPartitionsOfLength(100L,6L))
user system elapsed
4.57 0.36 4.93
system.time(res <- combsum(100L,6L))
user system elapsed
3.69 0.00 3.71
identical(smc(a),smc(res))
[1] TRUE
## this would take a very long time with GetDecimalReps below
注意：IntegerPartitionsOfLength只返回combinations一组特定数字的 ，而不是permutations一组数字的 （顺序很重要）。例如对于 set s = (1, 1, 3)， 的组合s是 精确s的，而 的排列s是：(1, 1, 3), (1, 3, 1), (3, 1, 1)。
如果您想要 OP 要求的答案，您将不得不做这样的事情（这绝不是最好的方法，也没有上面@bgoldst 的效率那么permsum高）：
library(gtools)
GetDecimalReps <- function(n) {
myPerms <- permutations(n,n); lim <- nrow(myPerms)
intParts <- IntegerPartitionsOfLength(100,n,FALSE)
do.call(rbind, lapply(intParts, function(x) {
unique(t(sapply(1L:lim, function(y) x[myPerms[y, ]])))
}))
}
system.time(a <- GetDecimalReps(4L))
user system elapsed
2.85 0.42 3.28
system.time(res <- combsum(100L,4L))
user system elapsed
1.35 0.00 1.34
head(a/100)
[,1] [,2] [,3] [,4]
[1,] 0.01 0.01 0.01 0.97
[2,] 0.01 0.01 0.97 0.01
[3,] 0.01 0.97 0.01 0.01
[4,] 0.97 0.01 0.01 0.01
[5,] 0.01 0.01 0.02 0.96
[6,] 0.01 0.01 0.96 0.02
tail(a/100)
[,1] [,2] [,3] [,4]
[156844,] 0.25 0.26 0.24 0.25
[156845,] 0.25 0.26 0.25 0.24
[156846,] 0.26 0.24 0.25 0.25
[156847,] 0.26 0.25 0.24 0.25
[156848,] 0.26 0.25 0.25 0.24
[156849,] 0.25 0.25 0.25 0.25
identical(smp(a),smp(res)) ## using the smp function created by @bgoldst
[1] TRUE
上面的@bgoldst 算法对于两种返回类型（即组合/排列）都具有优势。另请参阅上面@bgoldst 的出色基准。作为结束语，您可以通过简单地更改to和设置以返回列表来轻松修改IntegerPartionsOfLength以获得1:100该总和的100所有组合。干杯!k <= mk==Limk <= LimcombsOnly = FALSE