r - 在给定的长度上，所有可能的十进制数字(百分数)的总和为1

Question

考虑向量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
.
.
.

Answer 1

这是使用递归操作的方法：

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的选择占用一列，因此必须从 m减去1，以便递归调用将在其输出矩阵中减少产生一列。

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()。而且，正如约瑟夫（Joseph）所建议的，可以修改算法以产生组合，可以按以下步骤进行操作，直至达到名称 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的情况来说实际上还是太高了，但这不是问题，因为它只是在调用过程中将应用的两个上限之一序列生成。

第二个更改当然是在 x lambda中进行递归调用时将最新选择 l作为参数传递给 lapply()。

最后的改变是最棘手的。现在必须按如下方式计算选择顺序：

seq((s+m-1L)%/%m,min(l,s-m+1L))

必须将下限从 permsum()中使用的1提高到可能仍允许下降幅度组合的最低可能选择。当然，最低可能的选择取决于尚未生产多少根色谱柱。列越多，我们必须留出更多的“空间”以供将来选择。公式是对 s进行 m的整数除法，但是我们也必须有效地“舍入”，这就是为什么我在进行除法之前添加 m-1L的原因。我还考虑过进行浮点除法，然后调用 as.integer(ceiling(...))，但是我认为全整数方法要好得多。

例如，考虑 s=10,m=3的情况。要产生10剩余的总和，剩下3列，我们不能选择小于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