この質問が多く寄せられるのを見てきましたが、それに対する真の具体的な答えを見たことはありません. rand()したがって、C++のように乱数ジェネレーターを使用するときに「モジュロ バイアス」が正確に存在する理由を人々が理解するのに役立つことを願って、ここに投稿します。
質問する
60746 次
10 に答える
438
で定義された定数であるrand()0 から の間の自然数を選択する疑似乱数ジェネレータも同様です (の一般的な概要については、この記事を参照してください)。RAND_MAXcstdlibrand()
0 と 2 の間の乱数を生成したい場合はどうなるでしょうか。説明のために、RAND_MAXが 10 で、 を呼び出して 0 から 2 の間の乱数を生成するとしrand()%3ます。ただし、rand()%30 から 2 までの数字は同じ確率で生成されません。
Whenrand()は 0、3、6、または 9 を返します rand()%3 == 0。したがって、P(0) = 4/11
Whenrand()は 1、4、7、または 10 を返します rand()%3 == 1。したがって、P(1) = 4/11
Whenrand()は 2、5、または 8 を返します rand()%3 == 2。したがって、P(2) = 3/11
これは、0 と 2 の間の数値を等しい確率で生成しません。もちろん、範囲が小さい場合、これは最大の問題ではないかもしれませんが、範囲が大きい場合、これは分布をゆがめ、小さい数値にバイアスをかける可能性があります。
では、rand()%n0 から n-1 までの範囲の数値が等しい確率で返されるのはいつでしょうか? いつRAND_MAX%n == n - 1。この場合、以前の仮定rand()が 0 と等しい確率で数値を返すことに加えRAND_MAXて、n のモジュロ クラスも均等に分散されます。
では、この問題をどのように解決すればよいでしょうか。大まかな方法は、目的の範囲内の数値が得られるまで乱数を生成し続けることです。
int x;
do {
x = rand();
} while (x >= n);
しかし、これは の値が小さい場合には非効率的です。これは、範囲内の値を取得nする機会しかないため、平均して の呼び出しをn/RAND_MAX実行する必要があるためです。RAND_MAX/nrand()
nより効率的な数式のアプローチは、 で割り切れる長さの大きな範囲を取り、範囲内にあるRAND_MAX - RAND_MAX % n乱数を取得するまで乱数を生成し続け、モジュラスを取得することです。
int x;
do {
x = rand();
} while (x >= (RAND_MAX - RAND_MAX % n));
x %= n;
の値が小さい場合、 をn複数回呼び出す必要はめったにありませんrand()。
引用された作品とさらに読む：
于 2012-06-11T17:44:04.150 に答える
37
ランダムに選択し続けることは、バイアスを取り除く良い方法です。
アップデート
で割り切れる範囲で x を検索すれば、コードを高速化できますn。
// Assumptions
// rand() in [0, RAND_MAX]
// n in (0, RAND_MAX]
int x;
// Keep searching for an x in a range divisible by n
do {
x = rand();
} while (x >= RAND_MAX - (RAND_MAX % n))
x %= n;
上記のループは非常に高速である必要があり、平均で 1 回の繰り返しとなります。
于 2012-06-12T00:10:56.777 に答える
22
@ user1413793 は問題について正しいです。1 つの点を指摘することを除いて、これ以上説明するつもりはありません。はい、 の値が小さい場合nと の値が大きい場合RAND_MAX、モジュロ バイアスは非常に小さい可能性があります。ただし、バイアスを誘発するパターンを使用するということは、乱数を計算するたびにバイアスを考慮し、さまざまなケースに対してさまざまなパターンを選択する必要があることを意味します。また、選択を誤ると、導入されるバグは微妙であり、単体テストはほとんど不可能です。適切なツール ( などarc4random_uniform) を使用する場合と比較すると、作業が増えるだけでなく、作業が減ることはありません。より多くの作業を行い、より悪いソリューションを取得することは、特にほとんどのプラットフォームで毎回正しく行うことが簡単な場合は、ひどいエンジニアリングです。
残念ながら、ソリューションの実装はすべて正しくないか、本来よりも効率的ではありません。(各解決策には、問題を説明するさまざまなコメントがありますが、それらに対処するために修正された解決策はありません。) これは、答えを探し求めている人を混乱させる可能性があるため、ここでは既知の適切な実装を提供します。
繰り返しになりますが、最善の解決策はarc4random_uniform、それを提供するプラットフォームで使用するか、プラットフォーム用の同様の範囲のソリューション ( Random.nextIntJava など) を使用することです。それはあなたにコードの費用をかけずに正しいことをします。ほとんどの場合、これが正しい呼び出しです。
を持っていない場合はarc4random_uniform、オープンソースの力を利用して、より広い範囲の RNG の上にそれがどのように実装されているかを正確に確認できます (ar4randomこの場合、同様のアプローチが他の RNG 上でも機能する可能性があります)。
OpenBSD の実装は次のとおりです。
/*
* Calculate a uniformly distributed random number less than upper_bound
* avoiding "modulo bias".
*
* Uniformity is achieved by generating new random numbers until the one
* returned is outside the range [0, 2**32 % upper_bound). This
* guarantees the selected random number will be inside
* [2**32 % upper_bound, 2**32) which maps back to [0, upper_bound)
* after reduction modulo upper_bound.
*/
u_int32_t
arc4random_uniform(u_int32_t upper_bound)
{
u_int32_t r, min;
if (upper_bound < 2)
return 0;
/* 2**32 % x == (2**32 - x) % x */
min = -upper_bound % upper_bound;
/*
* This could theoretically loop forever but each retry has
* p > 0.5 (worst case, usually far better) of selecting a
* number inside the range we need, so it should rarely need
* to re-roll.
*/
for (;;) {
r = arc4random();
if (r >= min)
break;
}
return r % upper_bound;
}
同様のことを実装する必要がある人のために、このコードに関する最新のコミット コメントに注目する価値があります。
2**32 % upper_boundとして 計算するように arc4random_uniform() を変更し-upper_bound % upper_boundます。コードを簡素化し、ILP32 アーキテクチャと LP64 アーキテクチャの両方で同じにします。LP64 アーキテクチャでは、64 ビットの剰余の代わりに 32 ビットの剰余を使用することでわずかに高速化します。tech@ ok deraadt で Jorden Verwer によって指摘されました。djm や otto からの異議なし
Java 実装も簡単に見つけることができます (前のリンクを参照)。
public int nextInt(int n) {
if (n <= 0)
throw new IllegalArgumentException("n must be positive");
if ((n & -n) == n) // i.e., n is a power of 2
return (int)((n * (long)next(31)) >> 31);
int bits, val;
do {
bits = next(31);
val = bits % n;
} while (bits - val + (n-1) < 0);
return val;
}
于 2013-11-18T15:32:40.500 に答える
17
意味
モジュロ バイアスは、モジュロ演算を使用して出力セットを入力セットのサブセットに縮小する際の固有のバイアスです。一般に、出力セットのサイズが入力セットのサイズの約数でないときにモジュロ演算を使用する場合のように、入力セットと出力セットの間のマッピングが均等に分散されていない場合は常にバイアスが存在します。
このバイアスは、数値がビットの文字列 (0 と 1) として表されるコンピューティングでは特に回避が困難です。ランダム性の真にランダムなソースを見つけることも非常に困難ですが、この説明の範囲を超えています。この回答の残りの部分では、真にランダムなビットの無制限のソースが存在すると仮定します。
問題の例
これらのランダム ビットを使用してサイコロ (0 ～ 5) をシミュレートすることを考えてみましょう。6 つの可能性があるため、数値 6 を表すのに十分なビット (3 ビット) が必要です。残念ながら、ランダムな 3 ビットから 8 つの結果が得られます。
000 = 0, 001 = 1, 010 = 2, 011 = 3
100 = 4, 101 = 5, 110 = 6, 111 = 7
モジュロ 6 の値を取ることで、結果セットのサイズを正確に 6 に減らすことができますが、これはモジュロ バイアスの問題を提示します: 1100 が111生成され、1 が生成されます。このサイコロはロードされます。
考えられる解決策
アプローチ 0:
ランダムなビットに頼るのではなく、理論的には小さな軍隊を雇って 1 日中サイコロを振って結果をデータベースに記録し、各結果を 1 回だけ使用することができます。これは見た目と同じくらい実用的であり、とにかく真にランダムな結果が得られることはほとんどありません (しゃれが意図されています)。
アプローチ 1:
モジュラスを使用する代わりに、素朴だが数学的に正しい解決策は、得られた結果を破棄110し111、単純に 3 つの新しいビットで再試行することです。残念ながら、これは各ロールに 25% の確率で再ロールが必要になることを意味します。これには、各再ロール自体も含まれます。これは、ごく些細な用途を除いて、明らかに実用的ではありません。
アプローチ 2:
より多くのビットを使用する: 3 ビットの代わりに 4 ビットを使用します。これにより、16 の可能な結果が得られます。もちろん、結果が 5 より大きい場合はいつでも再ロールすると状況が悪化する (10/16 = 62.5%) ため、それだけでは役に立ちません。
2 * 6 = 12 < 16 であるため、12 未満の結果を安全に取得し、モジュロ 6 を減らして結果を均等に分配できます。他の 4 つの結果は破棄する必要があり、前のアプローチと同様に再ロールします。
最初は良さそうに聞こえますが、数学を確認してみましょう。
4 discarded results / 16 possibilities = 25%
この場合、1 ビット追加してもまったく役に立ちませんでした。
その結果は残念ですが、5 ビットでもう一度試してみましょう。
32 % 6 = 2 discarded results; and
2 discarded results / 32 possibilities = 6.25%
確実な改善ですが、多くの実際のケースでは十分ではありません。良いニュースは、より多くのビットを追加しても、破棄して再ロールする必要がある可能性が高まることは決してないということです. これはサイコロだけでなく、すべての場合に当てはまります。
ただし、示されているように、1 ビット余分に追加しても何も変わらない場合があります。実際、ロールを 6 ビットに増やしても、確率は 6.25% のままです。
これにより、さらに 2 つの質問が生じます。
- 十分なビットを追加した場合、廃棄の確率が減少するという保証はありますか?
- 一般的な場合、何ビットあれば十分ですか?
一般的な解決策
ありがたいことに、最初の質問に対する答えはイエスです。6 の問題は、2^x mod 6 が 2 と 4 の間で反転し、偶然にも互いに 2 の倍数であるため、偶数 x > 1 の場合、
[2^x mod 6] / 2^x == [2^(x+1) mod 6] / 2^(x+1)
したがって、6 はルールではなく例外です。同じ方法で 2 の累乗を連続して生成するより大きな係数を見つけることは可能ですが、最終的にはこれをラップアラウンドする必要があり、破棄される可能性は低くなります。
これ以上の証拠を提供しなくても、一般に、必要なビット数の 2 倍を使用すると、破棄される可能性が小さくなり、通常は取るに足らないものになります。
コンセプトの証明
以下は、OpenSSL の libcrypo を使用してランダムなバイトを提供するサンプル プログラムです。-lcryptoコンパイルするときは、ほとんどの人が利用できるはずのライブラリにリンクしてください。
#include <iostream>
#include <assert.h>
#include <limits>
#include <openssl/rand.h>
volatile uint32_t dummy;
uint64_t discardCount;
uint32_t uniformRandomUint32(uint32_t upperBound)
{
assert(RAND_status() == 1);
uint64_t discard = (std::numeric_limits<uint64_t>::max() - upperBound) % upperBound;
uint64_t randomPool = RAND_bytes((uint8_t*)(&randomPool), sizeof(randomPool));
while(randomPool > (std::numeric_limits<uint64_t>::max() - discard)) {
RAND_bytes((uint8_t*)(&randomPool), sizeof(randomPool));
++discardCount;
}
return randomPool % upperBound;
}
int main() {
discardCount = 0;
const uint32_t MODULUS = (1ul << 31)-1;
const uint32_t ROLLS = 10000000;
for(uint32_t i = 0; i < ROLLS; ++i) {
dummy = uniformRandomUint32(MODULUS);
}
std::cout << "Discard count = " << discardCount << std::endl;
}
MODULUSとの値をいじって、ROLLSほとんどの条件下で実際に何回リロールが発生するかを確認することをお勧めします。懐疑的な人は、計算された値をファイルに保存して、分布が正常に見えるかどうかを確認することもできます。
于 2015-04-23T00:10:33.927 に答える
10
modulo の使用には 2 つの通常の不満があります。
1 つはすべてのジェネレーターに有効です。極限の場合の方が見やすい。ジェネレーターの RAND_MAX が 2 (C 標準に準拠していない) で、値として 0 または 1 のみが必要な場合、モジュロを使用すると (ジェネレーターが 0 と 2 を生成する場合) 2 倍の頻度で 0 が生成されます。 1 を生成します (ジェネレーターが 1 を生成する場合)。これは、値をドロップしないとすぐに当てはまることに注意してください。ジェネレーターの値から必要な値へのマッピングに関係なく、一方が他方の 2 倍の頻度で発生します。
ある種のジェネレーターは、少なくともいくつかのパラメーターについて、重要度の低いビットが他のものよりもランダムではありませんが、悲しいことに、それらのパラメーターには他の興味深い特性があります (RAND_MAX を 2 のべき乗よりも 1 小さくすることができるなど)。この問題はよく知られており、長い間、ライブラリの実装はおそらく問題を回避してきました (たとえば、C 標準のサンプル rand() 実装ではこの種のジェネレータを使用していますが、重要度の低い 16 ビットを削除しています)。それとあなたは運が悪いかもしれません
のようなものを使用して
int alea(int n){
assert (0 < n && n <= RAND_MAX);
int partSize =
n == RAND_MAX ? 1 : 1 + (RAND_MAX-n)/(n+1);
int maxUsefull = partSize * n + (partSize-1);
int draw;
do {
draw = rand();
} while (draw > maxUsefull);
return draw/partSize;
}
0 から n の間の乱数を生成すると、両方の問題が回避されます (また、RAND_MAX == INT_MAX によるオーバーフローも回避されます)。
ところで、C++11 では、リダクションと rand() 以外のジェネレーターに標準的な方法が導入されました。
于 2012-06-13T09:54:11.537 に答える
-1
受け入れられた答えが示すように、「モジュロバイアス」は の低い値に根ざしていますRAND_MAX。彼は非常に小さな値RAND_MAX(10) を使用して、RAND_MAX が 10 の場合、% を使用して 0 から 2 の間の数値を生成しようとすると、次の結果になることを示しています。
rand() % 3 // if RAND_MAX were only 10, gives
output of rand() | rand()%3
0 | 0
1 | 1
2 | 2
3 | 0
4 | 1
5 | 2
6 | 0
7 | 1
8 | 2
9 | 0
したがって、0 の出力は 4 つ (4/10 の確率) あり、1 と 2 の出力は 3 つだけです (それぞれ 3/10 の確率)。
だから偏見です。数値が低いほど、出てくる可能性が高くなります。
RAND_MAXしかし、それは が小さい場合にのみ明らかに現れます。またはより具体的には、改造している数が に比べて大きい場合RAND_MAX。
ループ(これは非常に非効率的であり、提案する必要さえありません)よりもはるかに優れた解決策は、出力範囲がはるかに大きい PRNG を使用することです。Mersenne Twisterアルゴリズムの最大出力は 4,294,967,295 です。すべての意図と目的のためにこのようMersenneTwister::genrand_int32() % 10に行うと、均等に分散され、モジュロ バイアス効果はほとんどなくなります。
于 2013-04-15T02:19:03.210 に答える