5

ベクトル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たので、マトリックスの各行の合計が(ブルート フォース アプローチを除く)になるように、長さのすべての可能な順列のマトリックスが必要です。m1

たとえば、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
.
.
.
4

3 に答える 3

7

再帰を使用してこれを行う方法は次のとおりです。

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

生成されたシーケンスの各要素は、加算における加数の 1 つの可能な「選択」と見なす必要があります。

do.call(rbind,lapply(seq_len(s-m+1L),function(x) ...))

選択ごとに、再帰する必要があります。これを行うために使用できますlapply()。

先に進むために、ラムダは単一の再帰呼び出しをpermsum()行いcbind()、現在の選択で戻り値を返します。これにより、このレベルの再帰に対して常に同じ幅の行列が生成されます。したがって、lapply()呼び出しはすべて同じ幅の行列のリストを返します。次に、それらを一緒に行バインドする必要があるため、do.call(rbind,...)ここでトリックを使用する必要があります。

unname(cbind(x,permsum(s-x,m-1L)))

ラムダの本体はかなり単純です。再帰呼び出しの戻り値でcbind()現在の選択xを行い、この部分行列の合計を完成させます。残念ながら、 を呼び出す必要があります。そうしないと、引数unname()から設定される各列に列名が付けられます。xx

ここで最も重要な詳細は、再帰呼び出しに対する引数の選択です。まず、x現在の再帰評価中にラムダ引数が選択されたばかりなのでs、差し迫った再帰呼び出しが達成する責任を持つ新しい合計ターゲットを取得するために、ラムダ引数を減算する必要があります。したがって、最初の引数は になりs-xます。第 2 に、 の選択はx1 列を占めるため、 から 1 を引く必要があります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

したがって、2 つの結果の間には 42 行の不一致があります。次のコード行を使用して、省略された順列を正確に見つけようとしました。基本的に、2 つの行列の各要素を比較し、比較結果を論理行列として出力します。スクロールバックを下にスキャンして、最初の異なる行を見つけることができます。以下は、抜粋された出力です。

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シーケンス生成中に適用される 2 つの上限のうちの 1 つにすぎないため、問題ではありません。

2 番目の変更はもちろん、ラムダで再帰呼び出しを行うときに、最新の選択xを引数として渡すことでした。llapply()

最後の変更は最も厄介です。選択シーケンスは、次のように計算する必要があります。

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

下限は、使用されている 1 からpermsum()、降順の組み合わせを可能にする可能な限り低い選択まで引き上げる必要がありました。もちろん、可能な限り低い選択は、まだ作成されていない列の数によって異なります。列が多ければ多いほど、将来の選択のために残さなければならない「余地」が増えます。s式はonの整数除算を取ることですが、m効果的に「切り上げる」必要もありm-1Lます。浮動小数点除算を行ってから を呼び出すことも検討しましたas.integer(ceiling(...))が、すべて整数のアプローチの方がはるかに優れていると思います。

たとえば、 の場合を考えてみましょうs=10,m=3。残りの 3 列で 10 の合計を生成するには、4 未満の選択を行うことはできません。これは、組み合わせに沿って昇順でなければ 10 を生成するのに十分な数量がないためです。この場合、式は 12 を 3 で割り、4 を求めます。

の呼び出しを使用しpermsum()て現在の制限も適用する必要があることを除いて、上限は で使用したのと同じ式から計算できます。lmin()

次のコードを使用して、多くのランダムなテスト ケースで、新しい関数combsum()がジョセフの関数と同じように動作することを確認しました。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
于 2016-06-07T19:05:52.533 に答える
3

以下は、pure を返すアルゴリズムです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のみを返し、数値セットの は返しません(順序は重要です)。たとえば、集合の場合、 の組み合わせは正確にであり、 の順列は です。combinationspermutationss = (1, 1, 3)sss(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:100100k <= mk==Limk <= LimcombsOnly = FALSE

于 2016-06-07T21:33:29.040 に答える