擬似乱数（線形合同法・メルセンヌツイスタ）

← シミュレーション一覧

🎓 レベル：標準　|　重要度：A（必須）

📎 関連：逆関数法　|　準モンテカルロ法

要点（BLUF）

擬似乱数は、決定的な漸化式から「乱数らしく見える」数列を作ったもの。真の乱数ではないが、再現性があり高速。
線形合同法（LCG） $x_{n+1} = (a x_n + c) \bmod m$ は最も単純。必ず有限の周期を持ち（最大 $m$ ）、パラメータが悪いと格子構造が出ます。
実用はメルセンヌツイスタ（周期 $2^{19937}-1$ ）や PCG（numpy の既定）。小さな LCG の周期 16 をコードで確認します。

1. 擬似乱数とは

計算機は本質的に決定的なので、「真の乱数」は作れません。代わりに、決定的な漸化式から、統計的検定をパスする程度に乱数らしい数列を生成します。これが擬似乱数（pseudo-random number）です。

擬似乱数の良し悪しは次で測ります。

一様性： $[0,1)$ に均等に散らばるか。
無相関性：連続する値や $k$ 個飛ばしの値に相関がないか。
長い周期：同じ列が繰り返すまでの長さ。短いと使い尽くしてしまう。
再現性・速度：シードで同じ列を再現でき、高速に出せるか。

シードで再現できることは「欠点」ではなく長所です。デバッグ・検証・並列実行の管理に不可欠で、本テキストが全コードでシードを固定するのもこのためです。

2. 線形合同法（LCG）

最も古典的な生成器が線形合同法です。

x_{n+1} = (a\,x_n + c) \bmod m

整数の状態 $x_n$ を、乗数 $a$ ・増分 $c$ ・法 $m$ で更新し、 $u_n = x_n/m$ で $[0,1)$ に正規化します。状態は $0,\dots,m-1$ の $m$ 通りしかないので、周期は最大でも $m$ 。一度ある状態に戻れば、そこから先は同じ列の繰り返しです。

「最大周期 $m$ （フルピリオド）」を達成する条件は Hull–Dobell の定理で与えられ、(1) $c$ と $m$ が互いに素、(2) $m$ の全素因数が $a-1$ を割る、(3) $m$ が4の倍数なら $a-1$ も4の倍数、を満たすときです。

import numpy as np

def lcg(seed, a, c, m, n):
    """線形合同法で n 個の整数列を生成する"""
    out = np.empty(n, dtype=np.int64)
    x = seed
    for i in range(n):
        x = (a*x + c) % m        # 漸化式
        out[i] = x
    return out

# 小さな LCG（a=5, c=3, m=16）。Hull-Dobell条件を満たしフル周期になる
seq = lcg(seed=1, a=5, c=3, m=16, n=40)
first = seq[0]
period = next(i for i in range(1, len(seq)) if seq[i] == first)  # 戻るまでの歩数

print(f"最初の16個: {seq[:16]}")
print(f"周期 = {period}")

出力：

最初の16個: [ 8 11 10  5 12 15 14  9  0  3  2 13  4  7  6  1]
周期 = 16

出力の意味：最初の16個が 0〜15 の並べ替えになっていて、16歩で最初の状態に戻ります（周期 16 ＝フル周期）。 $m=16$ なのでこれ以上は不可能。実用の LCG は $m = 2^{32}$ や $2^{48}$ を使い周期を伸ばしますが、それでも周期は有限で、しかも下位ビットの規則性や格子構造（連続する値を多次元プロットすると平行平面に乗る）という弱点が残ります。

3. 一様性の確認と良いパラメータ

それなりのパラメータ（Numerical Recipes の $a=1664525,\ c=1013904223,\ m=2^{32}$ ）なら一様性は良好です。

import numpy as np

def lcg(seed, a, c, m, n):
    out = np.empty(n, dtype=np.int64); x = seed
    for i in range(n):
        x = (a*x + c) % m; out[i] = x
    return out

u = lcg(12345, a=1664525, c=1013904223, m=2**32, n=200_000) / 2**32
print(f"平均 = {u.mean():.4f}  (期待 0.5)")
print(f"分散 = {u.var():.4f}  (期待 {1/12:.4f})")

出力：

平均 = 0.4996  (期待 0.5)
分散 = 0.0836  (期待 0.0833)

出力の意味：平均 0.4996・分散 0.0836 が一様分布の理論値 0.5・ $1/12 = 0.0833$ とよく一致。1次元の一様性は合格です。ただし LCG は多次元での格子構造が弱点で、高次元のモンテカルロには不向き。実用では次の生成器を使います。

4. 実用的な生成器

メルセンヌツイスタ（MT19937）：周期 $2^{19937}-1$ 、623次元まで均等分布。長年の標準で numpy の旧 RandomState の既定。
PCG（Permuted Congruential Generator）：LCG に出力置換を足して統計的品質を高めた現代的生成器。numpy の default_rng() の既定で、本テキストもこれを使用。
暗号用途には予測不可能性が必要で、上記の擬似乱数は使えません（暗号論的 PRNG を使う）。

本テキストが np.random.default_rng(seed) を使うのは、品質・速度・再現性のバランスが良く、シードで完全再現できるためです。

数式の直観的意味

$x_{n+1} = (a x_n + c) \bmod m$ の $\bmod m$ は「 $m$ 個の引き出しの中をぐるぐる回る」操作です。状態が有限なので必ずいつか戻り、戻ったら同じ列が繰り返す——これが周期の正体。良い擬似乱数とは「この決定的な回り方が、統計検定では乱数と区別できないほど複雑」なものです。準モンテカルロ（準モンテカルロ法）は逆に、この規則性をわざと使って一様に敷き詰めます。

⚠️ よくある誤解・落とし穴

「擬似乱数は本物の乱数」ではない：決定的な数列です。だからシードで完全再現できます。
「周期が長ければ安心」ではない：周期だけでなく多次元一様性・無相関性も重要。LCG は周期を伸ばしても格子構造が残ります。
「適当な $a,c,m$ でいい」ではない：歴史的な RANDU（ $a=65539$ ）は深刻な格子構造で誤った結論を生みました。パラメータ選びは理論的条件に従う必要があります。
「シードを毎回時刻で変えるべき」ではない：研究・検証では固定が原則。再現できないシミュレーションはデバッグも査読もできません。
「並列でシードをずらせば独立」ではない：単純なシードずらしは部分列が重なる危険あり。並列には default_rng().spawn() など専用の仕組みを使います。

対応シミュレーション参照

本文の LCG 実装（周期16の確認・一様性チェック）。