数学編第32回乱数生成 ― 逆関数法と Box-Muller 法

当ページのリンクには広告も含まれています。

2026年9月25日

位置: 数学編第32回 / 前提: M15 / 数学編: なし / この回でできるようになること: 一様乱数から、逆関数法と Box-Muller 法で正規乱数を作ることができる

0. 物語 ― 0が出ると対数が発散する

くじ引きで「はずれ」を引くのは、よくあることです。しかし「ゼロ」と書かれた特別な札を引いた瞬間に、ゲームそのものが壊れてしまうとしたらどうでしょうか。コンピュータで乱数を作るとき、まさにそういうことが起こります。

コンピュータがまず用意できるのは、ゼロから1までの数をでたらめに出す仕組みです。サイコロの目のように偏りのない乱数で、これを一様乱数と呼びます。しかし、金融の計算で大量に必要になるのは、つりがね型をした別の種類の乱数です。コンピュータはつりがね型の乱数を直接は作れません。一様乱数を数学的に加工して、つりがね型に変換する必要があります。

変換の方法の一つに、2つの一様乱数を組み合わせて、2つのつりがね型の乱数を同時に作る手法があります。この手法では、片方の一様乱数の対数を途中で計算します。ところが、もしその乱数がちょうどゼロだったら、ゼロの対数は計算できません。値が際限なくマイナスに発散して、変換後の値も際限なく大きくなり、計算全体が使い物にならなくなります。

だから、一様乱数を作る範囲にゼロを含めてはいけません。「端を含む」範囲ではなく「端を含まない」範囲で乱数を取る必要があります。たった1つの値を除外するだけのことですが、それを忘れると結果が丸ごと壊れます。

この回では、一様乱数からつりがね型の乱数を作る2つの方法を証明します。そして、乱数の作り方の質が悪いと何が起こるかを、数値で確かめます。

1. 動機 ― コンピュータで正規分布をどう作るか

この回を読み終えると、一様乱数から、逆関数法と Box-Muller 法で正規乱数を作ることができます。

金融の計算では、正規分布に従う乱数が大量に必要になります。たとえば株価の将来の動きを何万回もシミュレーションするとき（モンテカルロ法と呼ばれ、本編第96回で扱います）、そのたびに正規乱数を生成しなければなりません。

しかし、コンピュータが直接作れる乱数は $[0, 1]$ 上の一様乱数だけです。正規分布のような複雑な分布は、一様乱数を数学的に変換して作ります。今回は、その変換の理論を2つ証明します。

この記事で証明すること

THM-M32-1 一様乱数に分布関数の逆関数を適用すると、目的の分布に従う乱数が得られます（逆関数法）。
THM-M32-2 2つの独立な一様乱数から、独立な2つの標準正規乱数を同時に作れます（Box-Muller 法）。
THM-M32-3 決定的な漸化式で作る擬似乱数には周期があり、生成法によっては高次元で規則的な構造が現れます。

2. 設定と仮定 ― 一様乱数から出発する

コンピュータの乱数生成は、次の2段階で構成されます。

第1段階: 決定的な漸化式（擬似乱数生成器）で $[0, 1]$ 上の一様乱数を作る。
第2段階: 一様乱数を数学的に変換して、目的の分布に従う乱数を作る。

第1段階の最も単純な例は線形合同法です。

x_{n + 1} = (a x_{n} + c) mod m (1.1)

ここで $a$ , $c$ , $m$ は整数の定数、 $x_{0}$ は初期値（シード）です。 $x_{n} / m$ を $[0, 1)$ 上の一様乱数として使います。

この方法は決定的なので、同じシードからは同じ列が出ます。周期は最大でも $m$ です。今回の主題は第2段階の数学ですが、第1段階の質が結果にどう影響するかも後半で確認します。

以下では次の仮定を置きます。

$U$ , $U_{1}$ , $U_{2}$ は $[0, 1]$ 上の一様分布に従う確率変数とします。
$U_{1}$ と $U_{2}$ は独立とします。
分布関数 $F$ は右連続で単調非減少、 $lim_{x \to - \infty} F (x) = 0$ , $lim_{x \to \infty} F (x) = 1$ を満たします。

一様乱数さえあれば、変換でほとんどの分布が作れます。問題はその変換をどう構成するかです。

3. 主張 ― 逆関数と変数変換で分布を変えられる

定理 THM-M32-1（逆関数法） $U$ が $[0, 1]$ 上の一様分布に従い、 $F$ が分布関数であるとする。一般化逆関数を

$F^{- 1} (u) = in f {x \in R : F (x) \geq u}$
と定義すると、 $X = F^{- 1} (U)$ は分布関数 $F$ を持つ確率変数である。

逆関数法は「一様乱数を分布関数の逆関数に通すだけで、目的の分布が得られる」という強力な結果です。ただし、すべての分布で逆関数が閉じた形で書けるとは限りません。

定理 THM-M32-2（Box-Muller 法） $U_{1}$ , $U_{2}$ が独立に $[0, 1]$ 上の一様分布に従うとする。

$Z_{1} Z_{2} = - 2 lo g U_{1} cos (2 π U_{2}), = - 2 lo g U_{1} sin (2 π U_{2})$
とおくと、 $Z_{1}$ と $Z_{2}$ は独立に標準正規分布 $N (0, 1)$ に従う。

定理 THM-M32-3（擬似乱数の周期と相関） 線形合同法 (1.1) のように決定的な漸化式で作る擬似乱数には、 $m$ 以下の周期がある。さらに、3次元以上の空間に連続する出力を並べると、点は有限個の超平面上にしか存在しない。

定理 THM-M32-3 は Marsaglia の定理として知られています。周期の短い生成器ほど、平面の枚数が少なくなり、規則的な構造が目立ちます。

証明を読む（読まなくても記事は完結します）

4. 証明 ― 変数変換の公式でヤコビアンを計算する

この節は読まなくても記事は完結します。

この証明を追いたい人は、数学編第15回「変数変換とヤコビアン」へ。

THM-M32-1 の証明

$X = F^{- 1} (U)$ の分布関数が $F$ であること、すなわち任意の $x \in R$ に対して $P (X \leq x) = F (x)$ を示します。

一般化逆関数の定義から、次の同値関係が成り立ちます。

F^{- 1} (u) \leq x ⟺ u \leq F (x) (2.1)

これを確認します。 $F^{- 1} (u) = in f {y : F (y) \geq u}$ なので、 $F^{- 1} (u) \leq x$ ということは、 $F (y) \geq u$ を満たす $y$ の下限が $x$ 以下だということです。 $F$ は右連続なので、 $F (x) \geq u$ が成り立ちます。逆に $F (x) \geq u$ なら、 $x$ は ${y : F (y) \geq u}$ の元なので、その下限は $x$ 以下です。

(2.1) を使って $X = F^{- 1} (U)$ の分布関数を計算します。

P (X \leq x) = P (F^{- 1} (U) \leq x) = P (U \leq F (x)) (2.2)

2行目は (2.1) の同値関係を適用しました。 $U$ は $[0, 1]$ 上の一様分布に従うので、 $0 \leq F (x) \leq 1$ の範囲では

P (U \leq F (x)) = F (x)

が成り立ちます。一様分布の分布関数が $[0, 1]$ 上で恒等関数であることを使いました。

したがって $P (X \leq x) = F (x)$ が示されました。 ∎

THM-M32-2 の証明

Box-Muller 法の証明は、極座標への変数変換とヤコビアンの計算で行います。

ステップ1: 変数変換の設定。

$U_{1}$ , $U_{2}$ は独立に $(0, 1)$ 上の一様分布に従います（ $U_{1} = 0$ は $lo g 0$ で破綻するため除外します）。変数を

R = - 2 lo g U_{1}, Θ = 2 π U_{2} (2.3)

と置きます。すると $Z_{1} = R cos Θ$ , $Z_{2} = R sin Θ$ です。

ステップ2: $R$ と $Θ$ の分布を求める。

$R$ の分布を導きます。 $U_{1}$ は $(0, 1)$ 上の一様分布なので、 $r > 0$ に対して

P (R \leq r) = P (- 2 lo g U_{1} \leq r) = P (- 2 lo g U_{1} \leq r^{2}) = P (U_{1} \geq e^{- r^{2} /2}) = 1 - e^{- r^{2} /2} (2.4)

3行目では、 $lo g$ の単調性により不等号の向きが変わることを使いました。4行目は一様分布の性質 $P (U_{1} \geq a) = 1 - a$ です。

したがって $R$ の密度関数は

f_{R} (r) = r e^{- r^{2} /2}, r > 0 (2.5)

です。これはレイリー分布です。

$Θ = 2 π U_{2}$ ですから、 $Θ$ は $[0, 2 π)$ 上の一様分布に従います。密度は $f_{Θ} (θ) = 1/ (2 π)$ です。 $U_{1}$ と $U_{2}$ が独立なので、 $R$ と $Θ$ も独立です。

ステップ3: 極座標から直交座標への変数変換。

$Z_{1} = R cos Θ$ , $Z_{2} = R sin Θ$ の同時密度を求めます。逆変換は

R = Z_{1}^{2} + Z_{2}^{2}, Θ = arctan \frac{Z _{2}}{Z _{1}}

ヤコビアンを計算します。 $(Z_{1}, Z_{2}) \to (R, Θ)$ の変換のヤコビ行列は

\frac{\partial ( Z _{1} , Z _{2} )}{\partial ( R , Θ )} = (cos Θ sin Θ - R sin Θ R cos Θ)

行列式は

\frac{\partial ( Z _{1} , Z _{2} )}{\partial ( R , Θ )} = R cos^{2} Θ + R sin^{2} Θ = R (2.6)

三角関数の基本恒等式 $cos^{2} Θ + sin^{2} Θ = 1$ を使いました。

ステップ4: 同時密度を計算する。

$(R, Θ)$ の同時密度に変数変換の公式を適用します。

f_{Z_{1}, Z_{2}} (z_{1}, z_{2}) = \frac{f _{R} ( r ) f _{Θ} ( θ )}{∣ J ∣} = \frac{r e ^{- r^{2} /2}}{2 π} \cdot \frac{1}{r} = \frac{1}{2 π} e^{- r^{2} /2}

ここで $r^{2} = z_{1}^{2} + z_{2}^{2}$ を代入して

f_{Z_{1}, Z_{2}} (z_{1}, z_{2}) = \frac{1}{2 π} e^{- (z_{1}^{2} + z_{2}^{2}) /2} (2.7)

を得ます。これは

f_{Z_{1}, Z_{2}} (z_{1}, z_{2}) = \frac{1}{2 π} e^{- z_{1}^{2} /2} \cdot \frac{1}{2 π} e^{- z_{2}^{2} /2}

と因数分解できます。右辺は標準正規密度 $φ (z_{1})$ と $φ (z_{2})$ の積です。同時密度が周辺密度の積に等しいので、 $Z_{1}$ と $Z_{2}$ は独立であり、それぞれが標準正規分布 $N (0, 1)$ に従います。 ∎

THM-M32-3 の証明（要点）

線形合同法 (1.1) の状態空間は ${0, 1, \dots, m - 1}$ の $m$ 個の値しかとりません。漸化式は決定的なので、同じ状態に戻った時点で列は周期的になります。したがって周期は $m$ 以下です。

Marsaglia（1968）は、線形合同法の連続 $d$ 個の出力を $d$ 次元空間にプロットすると、全ての点が高々 $(d! \cdot m)^{1/ d}$ 枚の超平面上に乗ることを示しました。 $m$ が小さいと平面の枚数が極めて少なくなり、乱数の「均一さ」が $d$ 次元では大きく損なわれます。 ∎

5. 反例 ― 質の悪い乱数は平面に並ぶ

CEX-M32-1: 正規分布には逆関数法が直接使えない

逆関数法は強力ですが、分布関数の逆関数が閉じた形で書ける場合にしか直接使えません。指数分布の場合は $F (x) = 1 - e^{- λ x}$ なので $F^{- 1} (u) = - lo g (1 - u) / λ$ と明示的に書けます。

しかし、標準正規分布の分布関数

N (x) = \frac{1}{2 π} \int_{- \infty}^{x} e^{- t^{2} /2} d t

は初等関数で表すことができません。したがって $N^{- 1} (u)$ も閉じた形では書けず、逆関数法を直接適用して正規乱数を作ることはできません。これが Box-Muller 法のような別の手法が必要になる理由です。

CEX-M32-2: 周期の短い線形合同法は3次元で破綻する

$a = 21$ , $c = 1$ , $m = 256$ の線形合同法を使います。周期は最大でも $256$ です。連続する3つの出力 $(x_{n}, x_{n + 1}, x_{n + 2})$ を3次元空間にプロットすると、254 個の点は少数の平面上に規則的に並びます。

先頭の5組を見ると

$n$	$x_{n}$	$x_{n + 1}$	$x_{n + 2}$
1	22	207	252
2	207	252	173
3	252	173	50
4	173	50	27
5	50	27	56

これらの点は $[0, 255]^{3}$ の立方体の中で一様に散らばるのではなく、特定の平面上に集中します。Marsaglia の定理から、3次元では高々 $(3! \times 256)^{1/3} \approx 11.54$ 枚、すなわち11枚以下の平行な平面にしか乗りません。実際、 $x_{n} + 2 x_{n + 1} - 3 x_{n + 2}$ の値は 6 種類しかなく、254 個の点はすべて 6 枚の平行な平面に乗っています。このような乱数で金融のシミュレーションを行うと、結果が系統的に偏ります。

質の悪い乱数を使っても結果は変わらない、という理解は誤りです。高次元の構造が計算結果を歪めます。

CEX-M32-3: 一様乱数に 0 が含まれると計算が壊れる

Box-Muller 法では $- 2 lo g U_{1}$ を計算します。もし $U_{1} = 0$ が生成されると

lo g 0 = - \infty

となり、 $- 2 \times (- \infty) = + \infty$ で計算が破綻します。したがって $U_{1}$ の生成範囲は閉区間 $[0, 1]$ ではなく、開区間 $(0, 1)$ でなければなりません。

実装では np.clip(U1, 1e-300, 1.0) のように下限を設けるか、 $U_{1} = 0$ を棄却して再生成するのが標準的な対処です。 $U_{1} = 0$ が出ても問題ない、という理解は誤りです。

6. 数値例 ― 正規乱数を作って検定する

compute.py で実際に乱数を生成し、理論値と比較します。

逆関数法による指数分布

$λ = 1$ の指数分布に対して $X = - lo g (1 - U)$ で 100,000 個の乱数を生成しました。

統計量	標本値	理論値
平均	$0.996$	$1.0$
分散	$0.986$	$1.0$

標本値は理論値に十分近く、逆関数法が正しく機能していることが確認できます。

Box-Muller 法による標準正規乱数

独立な一様乱数の組 $(U_{1}, U_{2})$ を 100,000 組生成し、Box-Muller 法で $Z_{1}$ , $Z_{2}$ を作りました。

統計量	$Z_{1}$ の標本値	$Z_{2}$ の標本値	理論値
平均	$0.0007$	$- 0.0019$	$0$
分散	$0.9931$	$0.9975$	$1$
歪度	$0.0066$	―	$0$
尖度	$2.9893$	―	$3$

$Z_{1}$ と $Z_{2}$ の標本相関係数は $0.0045$ であり、理論値 $0$ に近い値です。これは定理 THM-M32-2 の「 $Z_{1}$ と $Z_{2}$ は独立」という結論と整合しています。

すべての統計量が理論値の近傍にあり、Box-Muller 法が正しく標準正規乱数を生成していることが確認できます。

7. 次回へ ― 本編に戻って読む

数学編はこの第32回で一巡します。ここまでに扱った道具――測度、積分、確率変数、確率微分方程式、そして今回の乱数生成――は、本編の各回で前提として参照されます。本編の議論で数学的な基礎に立ち返りたくなったとき、対応する数学編に戻って読めるように設計されています。

本編第96回では、今回の乱数生成を使って金融商品の価格を計算する方法を扱います。

よくある質問

コンピュータの乱数は本当にランダムですか

ランダムではありません。決まった規則で作られますが、統計的にランダムに見えるように設計されています。

なぜ一様乱数から始めるのですか

一様乱数さえあれば、変換でほとんどの分布が作れるからです。

Box-Muller 法は今も使われていますか

使われています。より速い方法もありますが、考え方の基本として重要です。

参考文献

Glasserman P.、Monte Carlo Methods in Financial Engineering、Springer、2004年
Knuth D. E.、The Art of Computer Programming Volume 2 ― Seminumerical Algorithms、第3版、Addison-Wesley、1997年

数値例は compute.py によりシード固定で検算済み。