three.js + cannon-es でブラウザ上の 3D サイコロを作っています。物理演算で転がしているので「乱数をそのまま出目にする」より自然に見えますが、そのぶん出目が一様である保証はありません。そこで出目の偏りを χ²（カイ二乗）適合度検定で確認したところ、5 個同時に振ったときだけ 3 と 4 が有意に出にくいことがわかりました。原因の特定と修正、再検定までをまとめます。
検定の方法
- サイトの物理パラメータ（重力、摩擦、反発係数、減衰、スリープ判定、ステップ幅 1/120 s など）をそのまま Node.js に移し、描画なしのヘッドレスで大量に振る
- 出目の読み取りもアプリと同じ（+Y に最も近い面法線）
- 6 面が等確率という帰無仮説に対して Pearson の χ² 検定、自由度 5
- 臨界値: 11.07（p = 0.05）、15.09（p = 0.01）、20.52（p = 0.001）
注意点として、アリーナ形状は実機（カメラの視錐台から決まる）の近似です。また 5 個同時の場合、同じ投擲内のサイコロはぶつかり合うので独立ではなく、χ² はやや楽観的（有意に出やすい）になります。
χ² と p 値を計算するコード
外部ライブラリなしで動きます。p 値は χ² 分布の上側確率で、正則化不完全ガンマ関数 Q(df/2, χ²/2) として求めます。
// 適合度検定（一様分布）: χ² と p 値（上側確率）
function chiSquare(counts) {
const n = counts.reduce((a, b) => a + b, 0);
const e = n / counts.length; // 期待度数
return counts.reduce((s, o) => s + (o - e) ** 2 / e, 0);
}
// 正則化不完全ガンマ関数 Q(a, x) = 1 - P(a, x)
function gammaQ(a, x) {
const lnG = lnGamma(a);
if (x < a + 1) { // 級数展開
let sum = 1 / a, term = sum;
for (let k = 1; k < 500; k++) { term *= x / (a + k); sum += term; if (term < sum * 1e-15) break; }
return 1 - sum * Math.exp(-x + a * Math.log(x) - lnG);
}
// 連分数（Lentz 法）
let b = x + 1 - a, c = 1e300, d = 1 / b, h = d;
for (let i = 1; i < 500; i++) {
const an = -i * (i - a); b += 2;
d = an * d + b; if (Math.abs(d) < 1e-300) d = 1e-300;
c = b + an / c; if (Math.abs(c) < 1e-300) c = 1e-300;
d = 1 / d; const del = d * c; h *= del;
if (Math.abs(del - 1) < 1e-15) break;
}
return Math.exp(-x + a * Math.log(x) - lnG) * h;
}
function lnGamma(z) { // Lanczos 近似
const g = 7, p = [0.99999999999980993, 676.5203681218851, -1259.1392167224028,
771.32342877765313, -176.61502916214059, 12.507343278686905,
-0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7];
z -= 1; let x = p[0];
for (let i = 1; i < g + 2; i++) x += p[i] / (z + i);
const t = z + g + 0.5;
return 0.5 * Math.log(2 * Math.PI) + (z + 0.5) * Math.log(t) - t + Math.log(x);
}
const pValue = (chi2, df) => gammaQ(df / 2, chi2 / 2);
// 出目 1〜6 の度数
const runs = {
'v12 5個 (C+D)': [9923, 10222, 9505, 9845, 10383, 10122],
'v13 5個 (F+G)': [10018, 10027, 10049, 9951, 10078, 9877],
'v13 1個 (E)': [5040, 4996, 4975, 5022, 4998, 4969],
};
for (const [name, c] of Object.entries(runs)) {
const x = chiSquare(c);
console.log(name, 'χ²=' + x.toFixed(2), 'p=' + pValue(x, c.length - 1).toPrecision(3));
}
console.log('検算 χ²=11.07 ->', pValue(11.07, 5).toFixed(4)); // 0.0500
実行結果:
v12 5個 (C+D) χ²=48.58 p=2.70e-9
v13 5個 (F+G) χ²=2.71 p=0.745
v13 1個 (E) χ²=0.74 p=0.981
検算 χ²=11.07 -> 0.0500
修正前（v12）の結果
| ラン | 投擲数 | 1回の個数 | 出目 1〜6 | χ² | p |
|---|---|---|---|---|---|
| A | 12,000 | 1 | 2062 / 2014 / 1947 / 1895 / 2000 / 2082 | 12.30 | 0.031 |
| B | 30,000 | 1 | 5042 / 5048 / 4889 / 4954 / 5022 / 5045 | 4.20 | 0.52 |
| C | 6,000 | 5 | 5008 / 5045 / 4790 / 4895 / 5196 / 5066 | 20.00 | 0.0013 |
| D | 6,000 | 5 | 4915 / 5177 / 4715 / 4950 / 5187 / 5056 | 32.08 | <0.0001 |
| C+D | 12,000 | 5（60,000 個） | 9923 / 10222 / 9505 / 9845 / 10383 / 10122 | 48.58 | ≈3×10⁻⁹ |
1 個振りは微妙なところです。ラン B は余裕で通り、ラン A は p = 0.031 でぎりぎり棄却。弱い兆候はあるが決め手にはならない、という程度です。
5 個同時は独立した 2 回のランが両方とも強く棄却しました。合算すると 3 の目が 15.84%、4 の目が 16.41%（理想は 16.67%）。差は 1 ポイント未満なので遊んでいて気づくことはまずありませんが、6 万個振れば検出できます。
3 と 4 はこのモデルではちょうど ±X 軸の面です。特定の軸だけ出にくいということは、乱数ではなく向きに依存する処理が怪しいと考えられます。
原因：傾いて止まったサイコロを「平らにする」補正
物理演算のサイコロは、壁や他のサイコロに寄りかかって傾いたまま止まることがあります（面法線と +Y の内積 < 0.988）。v12 ではこれを次の順で補正していました。
- tip：小さく蹴って倒す
- nudge：それでもダメなら押し出す
- flatten：最後に強制的に平らにする
問題は 1 と 2 です。蹴るときの回転軸がサイコロ自身の向きに紐づいていたため、補正で倒れる先の面が軸によって偏っていました。5 個同時だとサイコロ同士で寄りかかるケースが多く、v12 の 5 個ランでは tip が 3,296 回、nudge が 392 回発生していました。1 個振りより 5 個振りで偏りが強く出たのはこのためと見ています。
v13 では次のように変えました。
- 傾いたサイコロは蹴らず、最寄りの平らな面へ 0.3 秒かけて slerp で回しながら、近くの空いた場所へスライドさせる
- 空きがないときだけ従来の tip / nudge にフォールバック
- 傾いたサイコロは高い位置のものから順に補正し、全アニメーションの完了を待ってから判定
- d6 のコライダーを角丸の ConvexPolyhedron から
CANNON.Boxに変更
結果、5 個 6,000 投擲あたりの tip は約 1,650 回から約 12 回に減りました。なお「補正経路が原因」というのは変更箇所からの推論で、各変更を個別に切り分けて検証したわけではありません。
修正後（v13）の結果
| テスト | v12 χ² (p) | v13 χ² (p) |
|---|---|---|
| 1 個, 30,000 投擲 | 4.20 (0.52) / 12,000 投擲で 12.30 (0.031) | 0.74 (0.98) |
| 5 個, ラン 1 | 20.00 (0.0013) | 1.54 (0.91) |
| 5 個, ラン 2 | 32.08 (<0.0001) | 5.98 (0.31) |
| 5 個 合算 | 48.58 (≈3×10⁻⁹)、3 の目 15.84% | 2.71 (0.75)、最小 16.46% |
すべてのランで一様性は棄却されなくなり、±X 面の不足も消えました。
結果の読み方で気をつけること
「有意差なし」は「公平だと証明された」ではありません。 p = 0.75 が意味するのは「一様だと仮定すれば、この程度のばらつきは普通に起こる」ということだけです。検出力が足りずに小さな偏りを見逃している可能性は常に残ります。言えるのは「この条件・このサンプル数では検出できる偏りはなかった」までです。
ほかにもハマりどころがいくつかあります。
- p は「偏りの大きさ」ではない。 サンプルを増やせば実用上無意味な差でも p は小さくなります。今回の v12 も差は 1 ポイント未満でした。比率を一緒に見るべきです。
- 何度も検定すれば偶然の棄却は起きる。 有意水準 5% で 20 回検定すれば 1 回は棄却されて当然です。v12 のラン A（p = 0.031）単独で「バグだ」と結論しなかったのはこのためで、独立したランで再現するかを確認しました。
- 独立性の仮定。 同じ投擲内のサイコロはぶつかり合うので厳密には独立ではありません。χ² は楽観的になります。
- シミュレーションは実機ではない。 アリーナは近似、タイミングは 60 fps 固定です。画面サイズによる差までは保証できません。
まとめ
- 物理演算のサイコロでも、偏りは物理そのものより「止まった後の補正処理」に潜んでいることがある
- 軸ごとに偏りが出たら、向きに依存する処理を疑う
- 1 回の有意な結果で判断せず、独立したランで再現を確認する
- 「有意差なし」は「偏りが検出されなかった」であって「公平の証明」ではない
みなさんは乱数や物理演算の出目の偏りをどう検証していますか？ 回帰テストとして CI に組み込んでいる例や、検出力の見積もり方などがあれば知りたいです。
検証に使ったのは自作のサイコロサイトです：https://3droll.com/ja/