Math:rnd(min, max)等の分布が偏っている The Math:rnd(min, max) and its seeded variant have slight biases in distribution

[MineCake147E]

概要

aiscript/src/interpreter/lib/std.ts

Lines 408 to 427 in 537fa96

	'Math:rnd': FN_NATIVE(([min, max]) => {
	if (min && min.type === 'num' && max && max.type === 'num') {
	return NUM(Math.floor(Math.random() * (Math.floor(max.value) - Math.ceil(min.value) + 1) + Math.ceil(min.value)));
	}
	return NUM(Math.random());
	}),
	'Math:gen_rng': FN_NATIVE(([seed]) => {
	expectAny(seed);
	if (seed.type !== 'num' && seed.type !== 'str') return NULL;
	const rng = seedrandom(seed.value.toString());
	return FN_NATIVE(([min, max]) => {
	if (min && min.type === 'num' && max && max.type === 'num') {
	return NUM(Math.floor(rng() * (Math.floor(max.value) - Math.ceil(min.value) + 1) + Math.ceil(min.value)));
	}
	return NUM(rng());
	});
	}),

現行の実装では、 $[0, 1)$ の乱数を生成し、範囲を $[min, max+1)$ に整形してからfloorするという方法で生成されています。
しかし、この方法では偏りが生じます。

現行実装の問題点

浮動小数点数演算の丸め誤差による偏り

まず、 $[0, 1)$ 内の有効な倍精度浮動小数点数は $4,607,182,418,800,017,408$ ( $2^{52} * 3 * 11 * 31$ )パターン存在します。
Math.random()はこれらの数の中から何らかの分布(選定されたパターンを倍精度浮動小数点数として解釈し観測した場合に一様分布となるような分布)に従ってランダムに値を選定して返します。
これを範囲変換する際に浮動小数点数の乗算や加算を行っています。
浮動小数点数では、四則演算を行う度に丸めが発生します。ここでの丸めは切り捨てや切り上げではなく、Rounding half to even(日本語版)と呼ばれる方法に基づいて行われることが多いです。
また、浮動小数点数はそもそも無限桁の演算を行うことが出来ない為、四則演算は殆どの場合単射ではありません。
ここで偏りが発生します。

内部的に生成されうる全パターンを等分出来ないことによる偏り

まず、 $[0, 1)$ 内の有効な倍精度浮動小数点数 $4,607,182,418,800,017,408$ パターンの内、その約99.8%は $[0, 0.5)$ の範囲内にあり、 $[0.5, 1)$ の範囲内にあるパターン数は $4,503,599,627,370,496$ ( $2^{52}$ )パターンしか存在しません。
Math.random()では $[0, 0.5)$ と $[0.5, 1)$ の範囲がそれぞれ50%の確率で生成されます。つまり、50%もの確率で $2^{52}$ パターンの中から抽選されることになります。
$2^{52}$ パターンの中から選ばれた $[0.5, 1)$ の範囲内にある数を2の冪ではない整数(range)で乗算し、Math.floor()により丸めを行うと、 $[\lfloor {\frac{range}{2}} \rfloor, range)$ 内の整数のいずれかにたどり着きますが、この範囲の整数はそもそも $2^{52}$ パターンを等分できず、偏りが発生します。
他の範囲についても同様のことが言えます。

提案手法

偏りを除去するためのアルゴリズムは複数存在します。

Rejection Sampling系アルゴリズム

最もシンプルな方法として、"Rejection Sampling"と呼ばれる方法が広く知られています。
以下に示すアルゴリズムはそのうちの一種類です。

1. $\lceil log_2(range + 1)\rceil$ bitの乱数列を得ることにより、 $[0, 2^{\lceil log_2(range + 1)\rceil})$ の範囲の整数を一様分布で生成する。
2. 1で生成した数がrangeを
   • 超える場合、1に戻る。
   • 超えない場合、その数をそのまま返す。

TypeScriptでの実装例:

min以上max以下の整数を生成する例(max-minは2^32 以下である必要があります):

const trueMin = Math.ceil(min);
const scale = Math.floor(max) - trueMin + 1;
const array = new Uint32Array(1);
const shift = Math.clz32(scale);
let result: number;
do {
    result = (crypto.getRandomValues(array)[0] ?? -1) >> shift;
} while (result > scale);
return Math.floor(Math.random() * scale + trueMin);

この方法は最も簡潔ですが、maxInclusiveが二の冪(=パターン数が二の冪+1)であるときに約50%もの確率で生成をやり直すことになり、殆どの乱数列が無駄になります。

Lemire's algorithm 及びその変種

.NET RuntimeやOpenSSL等で採用された方法です。

詳細な解説は以下にあるのでそれをご確認ください。

• .NET Runtime及びNode.jsで採用されたLemire's algorithm
• OpenSSLで採用されたLemire's algorithmの変種

どちらも固定小数点数をスケールして範囲に収めていますが、偏りの原因となるパターンが生成されたことを検知して精度を向上した上で乗算する処理が入っており、ループ回数を制限しなければ偏りを取り除くことが出来ます。

[FineArchs]

ユースケースによってはシード値の保存とかをしている可能性がある以上、乱数生成・調整アルゴリズムの変更は破壊的変更ということになりそうですかね？
ユーザー側で手法を選べるようにするか、そうでなければnext（破壊的変更用ブランチ）に入れるのがよさそう？

[FineArchs]

論文を読むのがかなり苦手なのであまり参加できないかも…

[syuilo]

Misskey Pages/Playとかでランダム生成なのに何故かパターンがあるような生成がされるなと3年くらい前から思ってるけどそういうことだったのかしら

[MineCake147E]

Rejection Samplingの方は既に手元環境で試してみました。実装例をほぼそのまま書き込みましたがうまく動いてくれます。

Misskey Pages/Playとかでランダム生成なのに何故かパターンがあるような生成がされるなと3年くらい前から思ってるけどそういうことだったのかしら

欲しいパターン数nに対して $2^{53}\mod n$がとても大きくなるような数ならそう感じるかもしれませんが...
このIssueではあくまでも得られた乱数列の処理方法に関するものです。
「パターンがあるような生成がされる」と感じる原因は、人間が(サイコロで生成したような)真の乱数列からでも(無意識に)何らかのパターンを見出そうとする特性によるものであるか、内部の擬似乱数生成器が生成する擬似乱数列自体の品質に致命的な問題があるかのどちらかである可能性が高いです。

内部の乱数生成器の品質(Off-topic?)

また、昨夜に解析していたのですが、seedrandomの内部のシード生成には予測が比較的容易なシード生成アルゴリズムが使われており、生成したシードを使うアルゴリズムもRC4という現在は(RFC 7465等により)非推奨とされている擬似乱数生成器が用いられているようです。(これについては別Issueで出すつもりでしたが...)

シード生成(暗号目的には決して使えないアルゴリズム)

https://github.com/davidbau/seedrandom/blob/4460ad325a0a15273a211e509f03ae0beb99511a/seedrandom.js#L179-L186

seedrandom

https://github.com/davidbau/seedrandom/blob/4460ad325a0a15273a211e509f03ae0beb99511a/seedrandom.js#L44-L54

RC4実装

https://github.com/davidbau/seedrandom/blob/4460ad325a0a15273a211e509f03ae0beb99511a/seedrandom.js#L106-L147

[FineArchs]

うーん、話を聞く限りよほど大きい数を扱わなければ分布の偏りは無視できる程度に収まりそうですかね？
大きな数向けのより高精度な乱数や、より予測しづらく暗号に使えるような乱数は別の関数として実装するのがいい気がします

[FineArchs]

この先更に精度が良かったり安全だったりなアルゴリズムを実装する時の対応を考えたいですね…

例えば、Rnd:rejection_samplingやRnd:lemiresのようにアルゴリズムが識別できるような命名の関数をいくつか用意して、
それらのうち、例えば

• 動作の軽い乱数の最新版をMath:rnd
• 精度のよい乱数の最新版をMath:rnd_fine
• 予測しづらい乱数の最新版をMath:rnd_safe

というようにエイリアスを用意することで、多くのユーザーにはより改善されたアルゴリズムを提供でき、アルゴリズム変更が気になるユーザーにも以前のアルゴリズムを使い続ける選択肢がある、というのはどうでしょうか？

[MineCake147E]

よほど大きい数を扱わなければ分布の偏りは無視できる程度に収まりそうですかね？

残念ながら僅かな偏りに噛みつく人はどうしても存在します。
特に麻雀やカードゲーム等のように運が強く作用するゲームにおいて乱数の品質や"偏り"はしばしば疑惑の目で見られてしまいます。
運試し系のPlayを作っている身としては「乱数の偏りのせいで◯◯が出ない」などと言われてしまうと結構傷付きます。
(このような発想の殆どはギャンブラーの誤謬等確率論に対する理解不足に起因するものであるとも言えますが...)
それでも、仮に内部の乱数生成にハードウェア乱数生成器を用いたところで、乱数列の不適切な取り扱いを原因とする偏りが生じてしまっては結局疑いの目で見られることになるので、(乱数生成器自体の特性はともかく)乱数生成器を適切に扱えるアルゴリズムを提供することは重要であると思います。

[FineArchs]

（私は乱数生成器とそれを扱うアルゴリズムを一括りにして語っていましたが、）

乱数生成器を適切に扱えるアルゴリズムを提供することは重要であると思います。

概ね同意です。
それとして、多少ばらつきはあるが軽量なアルゴリズムにも需要はあると思うので、やはり両方提供するのがよさそうですね。

[marihachi]

発生している乱数の偏りがどの程度問題になるのか分かりません。
修正しなければならない程度の偏りなんでしょうか？

例えば、その偏りが実用上あまり影響がないのであれば、修正しても良いとは思いますが必ずしも修正する必要はありません。AiScriptとしてその方法を採用するかを選択できると思います。

[FineArchs]

AiScriptとしてその方法を採用するかを選択できると思います。

AiScriptはサードパーティライブラリのような仕組みがないので、やるなら公式で実装するしか無い、という点で私は採用したい寄りですね

[marihachi]

仮に、実用上の影響があまりなくてできるだけ本来の乱数に近づけたいという目的であれば、方法があるならどこまでも近づけるのかという話になりますし、実装の複雑さに影響を与えます。AiScriptとしてこれを採用するかどうかの判断ができるという意味です。

[FineArchs]

もちろんその通りですし、その上で

AiScriptはサードパーティライブラリのような仕組みがないので、やるなら公式で実装するしか無い、という点で私は採用したい寄りですね

です。それが際限のない探求になるとしても、例えば

例えば、Rnd:rejection_samplingやRnd:lemiresのようにアルゴリズムが識別できるような命名の関数をいくつか用意して、 それらのうち、例えば

• 動作の軽い乱数の最新版をMath:rnd
• 精度のよい乱数の最新版をMath:rnd_fine
• 予測しづらい乱数の最新版をMath:rnd_safe

というようにエイリアスを用意することで、多くのユーザーにはより改善されたアルゴリズムを提供でき、アルゴリズム変更が気になるユーザーにも以前のアルゴリズムを使い続ける選択肢がある、というのはどうでしょうか？

のような仕組みを用意することで対処できると考えています。

[marihachi]

メンテナンスできる人は限られているのでそこも考える必要があります。
アルゴリズムの良し悪しを判断できる人はあまり多くありません。

[marihachi]

現在のアルゴリズムによる偏りの話です。

発生している乱数の偏りがどの程度問題になるのか分かりません。
修正しなければならない程度の偏りなんでしょうか？

[MineCake147E]

発生している乱数の偏りがどの程度問題になるのか分かりません。
修正しなければならない程度の偏りなんでしょうか？

用途にもよりますが、過去に以下のようにご回答した通り、修正は必要であると考えています。

残念ながら僅かな偏りに噛みつく人はどうしても存在します。
特に麻雀やカードゲーム等のように運が強く作用するゲームにおいて乱数の品質や"偏り"はしばしば疑惑の目で見られてしまいます。
運試し系のPlayを作っている身としては「乱数の偏りのせいで◯◯が出ない」などと言われてしまうと結構傷付きます。

また、.NET Runtimeでは「偏りが計測不可能なほど小さい」ことは「僅かに偏っている」のと同義であると結論付け、Lemire's algorithmを実装する際にループの実行回数に制限を設けませんでした。

[MineCake147E]

実装に取り掛かっているのでどなたかAssignして頂けませんか?

[MineCake147E]

後で気づきましたが、現行の実装では範囲外の値が帰ってきた結果Play等が意図せずクラッシュする可能性がありますね。

極稀に範囲外の値が生成される可能性がある(off-topic?)

Math.random()や、seedrandom("seed")()は、 $[0, 1)$ の範囲の数値を返すと記されています。
1より小さい最大の倍精度浮動小数点数は $0.99999999999999988897769753748434595763683319091796875$ ( $= 1.0 - 2^{-53}$ )(JavaScriptでの表示は0.9999999999999999)です。
この値が実際に返された場合、minが非ゼロの場合に範囲外の値が生成される恐れがあります。

aiscript/src/interpreter/lib/std.ts

Line 410 in 537fa96

return NUM(Math.floor(Math.random() * (Math.floor(max.value) - Math.ceil(min.value) + 1) + Math.ceil(min.value)));

maxが333、minが253、乱数生成器の生成結果が1 - (2 ** -53)と仮定して、これを上記の方法で変換すると、

$$\displaylines{\lfloor (1 - 2^{-53}) \times (333 - 253 + 1) + 253 \rfloor \newline = \lfloor (1 - 2^{-53}) \times (81) + 253 \rfloor \newline = RoundFraction_{52}(80.99999999999999100719350053623202256858348846435546875) + 253 \rfloor \newline = \lfloor 80.9999999999999857891452847979962825775146484375 + 253 \rfloor \newline = \lfloor RoundFraction_{52}(333.9999999999999857891452847979962825775146484375) \rfloor = \lfloor 334 \rfloor = 334}$$

となり、指定したmaxを超えてしまいます。
これは上でもご説明したRounding half to even(日本語版)によるものです。
この挙動そのものは

return NUM(Math.floor(Math.random() * (Math.floor(max.value) - Math.ceil(min.value) + 1)) + Math.ceil(min.value)); 

と書き直せば解決しますが、乱数生成器が0.98765432098765415513952348192106001079082489013671875(JavaScriptでの表示は0.9876543209876542)を返した場合等、一部の値が違う整数にマップされることになるため、破壊的変更は避けられません。

$$\displaylines{\lfloor 0.98765432098765415513952348192106001079082489013671875 \times (333 - 253 + 1) + 253 \rfloor \newline = \lfloor 0.98765432098765415513952348192106001079082489013671875 \times (81) + 253 \rfloor \newline = \lfloor RoundFraction_{52}(79.99999999999998656630140203560586087405681610107421875) + 253 \rfloor \newline = \lfloor 79.9999999999999857891452847979962825775146484375 + 253 \rfloor \newline = \lfloor RoundFraction_{52}(332.9999999999999857891452847979962825775146484375) \rfloor = \lfloor 333 \rfloor = 333}$$

だったのが、この変更により、

$$\displaylines{\lfloor 0.98765432098765415513952348192106001079082489013671875 \times (333 - 253 + 1) \rfloor + 253 \newline = \lfloor 0.98765432098765415513952348192106001079082489013671875 \times (81) \rfloor + 253 \newline = \lfloor RoundFraction_{52}(79.99999999999998656630140203560586087405681610107421875) \rfloor + 253 \newline = \lfloor 79.9999999999999857891452847979962825775146484375 \rfloor + 253 \newline = 79 + 253 = 332}$$

へと変わってしまいます。

[FineArchs]

乱数の最大値付近が出た時、丸め誤差が掛け算により増幅され、結果が最大値を超す場合がある、みたいな感じですかね？
普通に使っていればそうそう発生する事象ではないとは思いますが、もし修正するのであれば破壊的変更ついでにアルゴリズムごと変えてしまうのは確かにアリですね。バグ修正の体で次のリリースに入れてもいいと思います

[MineCake147E]

乱数の最大値付近が出た時、丸め誤差が掛け算により増幅され、結果が最大値を超す場合がある、みたいな感じですかね？

どちらかというと、掛け算まではぎりぎり切り上げられなかった結果を足し算後の丸めが切り上げさせてしまう感じですね。

普通に使っていればそうそう発生する事象ではないとは思いますが、もし修正するのであれば破壊的変更ついでにアルゴリズムごと変えてしまうのは確かにアリですね。

私もそう思います。(とりあえず現在は新しい関数として実装してますが...)