ベクトルを作成するために数値をサンプリングしたい場合は、次のようにします。
set.seed(123)
x <- sample(1:100,200, replace = TRUE)
sum(x)
# [1] 10228
合計が100になる20個の乱数をサンプリングし、次に合計が100になる30個の乱数をサンプリングしたい場合はどうなりますか。これは、見た目よりも難しいと思います。?sampleグーグルを検索しても手がかりは得られなかった。そして、サンプリングするループは、希望する合計に十分に近づかない場合（たとえば、5以内）に拒否するため、時間がかかる場合があります。
これを達成するためのより良い方法はありますか？
例は次のとおりです。
foo(10,100) # ten random numbers that sum to 100. (not including zeros)
# 10,10,20,7,8,9,4,10,2,20
Rを使用した試み
# Config
n <- 20L
target <- 100L
vec <- seq(100)
set.seed(123)
# R repeat loop
sumto_repeat <- function(vec,n,target) {
res <- integer()
repeat {
cat("begin:",sum(res),length(res),"\n")
res <- c( res, sample(vec,1) )
if( sum(res)<target & length(res)==(n-1) ) {
res[length(res)+1] <- target - sum(res)
}
# cat("mid:",sum(res),length(res),"\n")
if(sum(res)>target) res <- res[-length(res)]
if( length(res)>n | length(res)<n & sum(res)==target ) {
res <- res[-sample(seq(length(res)),1)]
}
# cat("end:",sum(res),length(res),"\n")
# cat(dput(res),"\n")
if( sum(res)==target & length(res)==n ) break
}
res
}
test <- sumto_repeat(vec=vec,n=n,target=target)
> sum(test)
[1] 100
> length(test)
[1] 20
また、どのディストリビューションから描画したいかについても考えてみます。要素と正確に合計する方法はいくつかあると思います（たとえば、最後の要素を常ににすることができます）。これは、分布に異なる影響を与える場合と持たない場合があります。targetntarget - sum(res)
Rcppの非常によく似たアルゴリズムです。
cpp_src <- '
Rcpp::IntegerVector xa = clone(x); // Vector to be sampled
Rcpp::IntegerVector na(n); // Number of elements in solution
Rcpp::IntegerVector sa(s); // Sum of solution
int nsampled;
int currentSum;
int dropRandomIndex;
int numZeroes;
Rcpp::IntegerVector remainingQuantity(1);
int maxAttempts = 100;
// Create container for our results
Rcpp::IntegerVector res(maxAttempts);
std::fill( res.begin(), res.end(), NA_INTEGER );
// Calculate min/max so that we can draw random integers from within range
Rcpp::IntegerVector::iterator mn = std::min_element(xa.begin(), xa.end()) ;
Rcpp::IntegerVector::iterator mx = std::max_element(xa.begin(), xa.end()) ;
std::cout << "mx = " << *mx << std::endl;
// Now draw repeatedly
nsampled = 0;
for( int i = 0; i < maxAttempts; i++ ) {
std::cout << "\\n" << i;
int r = *mn + (rand() % (int)(*mx - *mn + 1));
res[i] = xa[r+1];
// Calculate n and s for current loop iteration
numZeroes = 0;
for( int j = 0; j < maxAttempts; j++)
if(res[j]==0) numZeroes++;
std::cout << " nz= " << numZeroes ;
nsampled = maxAttempts - sum( is_na(res) ) - numZeroes - 1;
currentSum = std::accumulate(res.begin(),res.begin()+i,0); // Cant just use Rcpp sugar sum() here because it freaks at the NAs
std::cout << " nsamp= " << nsampled << " sum= " << currentSum;
if(nsampled == na[0]-1) {
std::cout << " One element away. ";
remainingQuantity[0] = sa[0] - currentSum;
std::cout << "remainingQuantity = " << remainingQuantity[0];
if( (remainingQuantity[0] > 0) && (remainingQuantity[0]) < *mx ) {
std::cout << "Within range. Prepare the secret (cheating) weapon!\\n";
std::cout << sa[0] << " ";
std::cout << currentSum << " ";
std::cout << remainingQuantity[0] << std::endl;
if( i != maxAttempts ) {
std::cout << "Safe to add one last element on the end. Doing so.\\n";
res[i] = remainingQuantity[0];
}
currentSum = sa[0];
nsampled++;
if(nsampled == na[0] && currentSum == sa[0]) std::cout << "It should end after this...nsamp= " << nsampled << " and currentSum= " << currentSum << std::endl;
break;
} else {
std::cout << "Out of striking distance. Dropping random element\\n";
dropRandomIndex = 0 + (rand() % (int)(i - 0 + 1));
res[dropRandomIndex] = 0;
}
}
if(nsampled == na[0] && currentSum == sa[0]) {
std::cout << "Success!\\n";
for(int l = 0; l <= i+1; l++)
std::cout << res[l] << " " ;
break;
}
if(nsampled == na[0] && currentSum != sa[0]) {
std::cout << "Reached number of elements but sum is ";
if(currentSum > sa[0]) {
std::cout << "Too high. Blitz everything and start over!\\n";
for(int k = 0; k < res.size(); k++) {
res[k] = NA_INTEGER;
}
} else {
std::cout << "Too low. \\n";
}
}
if( nsampled < na[0] && currentSum >= sa[0] ) {
std::cout << "Too few elements but at or above the sum cutoff. Dropping a random element and trying again.\\n";
dropRandomIndex = 0 + (rand() % (int)(i - 0 + 1));
res[dropRandomIndex] = 0;
}
}
return res;
'
sumto <- cxxfunction( signature(x="integer", n="integer", s="integer"), body=cpp_src, plugin="Rcpp", verbose=TRUE )
testresult <- sumto(x=x, n=20L, s=1000L)
testresult <- testresult[!is.na(testresult)]
testresult <- testresult[testresult!=0]
testresult
cumsum(testresult)
length(testresult)
いくつかの異なる値で試してみましたが、逃げない限り有効な答えが得られます。ここで注意点があります。それは、必要な数の要素から1つ離れていて、「打撃距離」内にある場合にチートすることです。たとえば、最後の値を描画するだけでなく、その数が有効かどうかを計算します。
ベンチマーク
比較コードについては、要点を参照してください。
別のアプローチですが、浮動小数点数を使用しているため、探しているものとは異なります。申し訳ありません。
randomsum <- function(nb, sum) {
tmp <- sort(runif(nb-1))
tmp <- c(min(tmp), diff(tmp), 1-max(tmp))
as.vector(quantile(0:sum, probs=tmp))
}
たとえば、次のようになります。
R> result <- randomsum(10, 1000)
R> result
[1] 35.282191 66.537308 17.263761 182.837409 120.064363 210.752735
[7] 143.201079 6.164731 34.936359 182.960064
R> sum(result)
[1] 1000
結果を使用roundして整数を取得することもできますが、もちろん、合計は取得したいものとわずかに異なる可能性があります。手っ取り早い回避策は、ランダムな値の1つを変更して、ベクトルの合計を必要な数にすることです。
randomsumint <- function(nb, sum) {
tmp <- sort(runif(nb-1))
tmp <- c(min(tmp), diff(tmp), 1-max(tmp))
res <- as.vector(quantile(0:sum, probs=tmp))
res <- round(res)
res[length(res)] <- res[length(res)]+(sum-sum(res))
res
}
どちらが与えるでしょう：
R> result <- randomsumint(10,1000)
R> result
[1] 42 152 0 11 74 138 9 138 172 264
R> sum(result)
[1] 1000
まれに、結果に負の値が生じる可能性があるため、これが完全とはほど遠いわけではありません。
これが別の試みです。は使用しませんsampleが、を使用しますrunif。合計を示すオプションの「メッセージ」を出力に追加しました。これは、showSum引数を使用してトリガーできます。Toleranceターゲットにどれだけ近いかを指定する引数もあります。
SampleToSum <- function(Target = 100, VecLen = 10,
InRange = 1:100, Tolerance = 2,
showSum = TRUE) {
Res <- vector()
while ( TRUE ) {
Res <- round(diff(c(0, sort(runif(VecLen - 1)), 1)) * Target)
if ( all(Res > 0) &
all(Res >= min(InRange)) &
all(Res <= max(InRange)) &
abs((sum(Res) - Target)) <= Tolerance ) { break }
}
if (isTRUE(showSum)) cat("Total = ", sum(Res), "\n")
Res
}
下記は用例です。
デフォルト設定と設定の違いに注意してくださいTolerance = 0
set.seed(1)
SampleToSum()
# Total = 101
# [1] 20 6 11 20 6 3 24 1 4 6
SampleToSum(Tolerance=0)
# Total = 100
# [1] 19 15 4 10 1 11 7 16 4 13
を使用して、この動作を確認できますreplicate。Tolerance = 0関数を5回設定して実行した結果は次のとおりです。
system.time(output <- replicate(5, SampleToSum(
Target = 1376,
VecLen = 13,
InRange = 10:200,
Tolerance = 0)))
# Total = 1376
# Total = 1376
# Total = 1376
# Total = 1376
# Total = 1376
# user system elapsed
# 0.144 0.000 0.145
output
# [,1] [,2] [,3] [,4] [,5]
# [1,] 29 46 11 43 171
# [2,] 103 161 113 195 197
# [3,] 145 134 91 131 147
# [4,] 154 173 138 19 17
# [5,] 197 62 173 11 87
# [6,] 101 142 87 173 99
# [7,] 168 61 97 40 121
# [8,] 140 121 99 135 117
# [9,] 46 78 31 200 79
# [10,] 140 168 146 17 56
# [11,] 21 146 117 182 85
# [12,] 63 30 180 179 78
# [13,] 69 54 93 51 122
Tolerance = 5また、機能を5回設定して実行する場合も同様です。
system.time(output <- replicate(5, SampleToSum(
Target = 1376,
VecLen = 13,
InRange = 10:200,
Tolerance = 5)))
# Total = 1375
# Total = 1376
# Total = 1374
# Total = 1374
# Total = 1376
# user system elapsed
# 0.060 0.000 0.058
output
# [,1] [,2] [,3] [,4] [,5]
# [1,] 65 190 103 15 47
# [2,] 160 95 98 196 183
# [3,] 178 169 134 15 26
# [4,] 49 53 186 48 41
# [5,] 104 81 161 171 180
# [6,] 54 126 67 130 182
# [7,] 34 131 49 113 76
# [8,] 17 21 107 62 95
# [9,] 151 136 132 195 169
# [10,] 194 187 91 163 22
# [11,] 23 69 54 97 30
# [12,] 190 14 134 43 150
# [13,] 156 104 58 126 175
当然のことながら、許容値を0に設定すると、関数が遅くなります。
これは「ランダムな」プロセスであるため、数字の正しい組み合わせを見つけるのにかかる時間を推測するのは難しいことに注意してください。たとえば、を使用してset.seed(123)、次のテストを3回続けて実行しました。
system.time(SampleToSum(Target = 1163,
VecLen = 15,
InRange = 50:150))
最初の実行には9秒強かかりました。1秒は7.5秒強かかりました。3番目は...381秒弱でした！それはたくさんのバリエーションです！
好奇心から、関数にカウンターを追加しました。最初の実行では、55026回試行して、すべての条件を満たすベクトルに到達しました。（私は2回目と3回目の試みをわざわざ試みませんでした。）
入力が妥当であることを確認するために、関数にエラーまたは健全性チェックを追加するとよい場合があります。たとえば、SampleToSum(Target = 100, VecLen = 10, InRange = 15:50)15から50の範囲では、100に到達する方法がなく、ベクトルに10の値があるため、入力できないはずです。
整数が必要であると仮定すると（そうでない場合はディリクレ分布を見てください）、これはボールと壺問題と考えることができます（数値間の関係にさらに制限はありません）。
20個の数字が必要な場合は、20個の壷で表すことができます。数字の合計を100にして、100個のボールにします。正確に20個の数字が必要なので（最大20個の数字が必要な場合はこの手順をスキップしますが、それより少なくてもかまいません）、まず各壷に1つのボールを置き、残りのボールを壷にランダムに分配します。各壷のボールの数を数えると、合計で100になる20個の数字が得られます。
Rコードとして：
as.vector(table( c( 1:20, sample(1:20, 80, replace=TRUE) ) ))
as.vectorテーブルクラスとラベルを削除するだけです。
迅速、シンプル、正確、ループなし、再帰など。
その他の合計または値の数については、上記の適切な部分を変更してください。
私は組み合わせ論における星とバーとパーティションについて考えました：
foo <- function(n,total) {
while(!exists("x",inherits=FALSE) || 1 %in% diff(x)) {
x <- sort(c(0,sample.int(n+total,n-1,replace=FALSE),n+total))
}
print(x)
sort(diff(x)-1)
}
もう1つの方法は、partitionsパッケージを使用することです。これは、すべてのパーティションを列挙するのに適していますが、今のところは問題ありません。総数が少ない限り動作します。
require(partitions)
foo <- function(n,total) {
x <- restrictedparts(total,n,include.zero=FALSE)
return(x[,sample.int(ncol(x),1)])
}