乱数生成の基本アルゴリズムまとめ（線形合同法、M系列、中心極限定理の応用 etc）

乱数はゲームのプログラミングやMCMCなどを用いた近似解の計算など、様々な分野で用いられています。パッケージを用いた乱数の発生についてはよくまとめを見かける一方で、乱数生成の原理についてまとめてある記事は少ないです。そこで当稿では「自然科学の統計学」などを主に参考にしつつ、乱数生成のアルゴリズムについてまとめます。
乱数を用いた計算は「自然科学の統計学」の出版以降に大きく発展したので、基本的な原理にとどまらず昨今用いられるアルゴリズムについても確認します。

乱数生成手法の概要

乱数の生成にあたってはサイコロなどを用いて乱数表(table of random number)を生成しても良いが、円周率の計算にあたっての数値積分のようにコンピュータを用いた計算にあたって大量の乱数が必要とされる場合には不向きである。

よって数値積分などに用いる際は、プログラムに従って乱数を生成する擬似乱数(算術乱数)が用いられることが基本的に多い。この擬似乱数(算術乱数)の基本的なアルゴリズムを確認することを当稿の目的とする。

一様乱数の生成

線形合同法(LCG)

$1948$年にレーマー(Lehmer)が考案した線形合同法(LCG; linear congruential method)について確認する。線形合同法は、漸化式に基づいて乱数列を出力する手法である。
$$
\large
\begin{align}
X_{n+1} = aX_{n} + c \quad (\mathrm{mod} \; M), \quad n \geq 1
\end{align}
$$
線形合同法では上記の漸化式を用いて、乱数列の$X_0, X_1, X_2, …$を出力する。式に記載のある「$\mathrm{mod} \; M$」は知らないと難しく見えるが、単に「$M$で割った余り」を意味する。「$2=5 (\mathrm{mod} \; 3)$」のイメージで理解しておけば良い。

数式だけではイメージがつかないと思われるので、$X_0=2$、$a=3$、$c=2$、$M=7$の簡単な例で確認を行う。
$$
\large
\begin{align}
X_{0+1} &= aX_{0} + c \quad (\mathrm{mod} \; M) \\
X_1 &= 3X_{0} + c \quad (\mathrm{mod} \; M) \\
&= 3 \times 2 + 2 \quad (\mathrm{mod} \; 7) \\
&= 8 \quad (\mathrm{mod} \; 7) \\
&= 1 \quad (\mathrm{mod} \; 7)
\end{align}
$$

上記のように$X_1=1$が計算できる。以下同様に$X_7$まで計算する。
$$
\large
\begin{align}
X_{2} &= 3X_{1} + 2 \quad (\mathrm{mod} \; 7) \\
&= 3 \times 1 + 2 \quad (\mathrm{mod} \; 7) \\
&= 5 \quad (\mathrm{mod} \; 7)
\end{align}
$$
$$
\large
\begin{align}
X_{3} &= 3X_{2} + 2 \quad (\mathrm{mod} \; 7) \\
&= 3 \times 5 + 2 \quad (\mathrm{mod} \; 7) \\
&= 17 \quad (\mathrm{mod} \; 7) \\
&= 3 \quad (\mathrm{mod} \; 7)
\end{align}
$$
$$
\large
\begin{align}
X_{4} &= 3X_{3} + 2 \quad (\mathrm{mod} \; 7) \\
&= 3 \times 3 + 2 \quad (\mathrm{mod} \; 7) \\
&= 11 \quad (\mathrm{mod} \; 7) \\
&= 4 \quad (\mathrm{mod} \; 7)
\end{align}
$$
$$
\large
\begin{align}
X_{5} &= 3X_{5} + 2 \quad (\mathrm{mod} \; 7) \\
&= 3 \times 4 + 2 \quad (\mathrm{mod} \; 7) \\
&= 14 \quad (\mathrm{mod} \; 7)\\
&= 0 \quad (\mathrm{mod} \; 7)
\end{align}
$$
$$
\large
\begin{align}
X_{6} &= 3X_{5} + 2 \quad (\mathrm{mod} \; 7) \\
&= 3 \times 0 + 2 \quad (\mathrm{mod} \; 7) \\
&= 2 \quad (\mathrm{mod} \; 7)
\end{align}
$$
$$
\large
\begin{align}
X_{7} &= 3X_{6} + 2 \quad (\mathrm{mod} \; 7) \\
&= 3 \times 2 + 2 \quad (\mathrm{mod} \; 7) \\
&= 8 \quad (\mathrm{mod} \; 7) \\
&= 1 \quad (\mathrm{mod} \; 7)
\end{align}
$$

上記のように計算することができる。ここで、$X_0=X_6=2$、$X_1=X_7=1$になっていることは着目すべきで、これは線形合同法によって生成される乱数の周期性を意味する。また、線形合同法によって生成される乱数は現在の値にのみ基づいて次の値が決まるので、「$\mathrm{mod} \; 7$」では最大$7$通りの値しか持たないことから、周期は$7$以下となることがわかる。ここでの周期は$6$であり、$7$以下である。

「$\mathrm{mod} \; M$」とした際に乱数の周期が$M$以下となるのは同様に考えることで一般的に成立するとできる。また、乱数列${X_n}$の周期が最大値の$M$となる必要十分条件は下記であることも知られている。

i) cとMが互いに素
ⅱ) b=a-1がMを割り切る全ての素数の倍数
ⅲ) Mが4の倍数であればbも4の倍数

最大周期を実現するパラメータの設定に関する証明などは下記で取り扱いを行った。

一様乱数の生成に用いる線形合同法(LCG)の最大周期の法則の証明

また、線形合同法では多次元乱数を生成する際に「多次元疎結晶構造」を生じやすいことに注意が必要である。多次元疎結晶構造に関しては下記で詳しく取りまとめた。

線形合同法(linear congruential method)における多次元疎結晶構造

乗算型合同法

乗算型合同法の基本的な考え方は線形合同法と同じであるが、乗算型合同法においては下記の漸化式を用いる。
$$
\large
\begin{align}
X_{n+1} = aX_{n} \quad (\mathrm{mod} \; M), \quad n \geq 1
\end{align}
$$
線形合同法において定数項の$c$を消去すると乗算型合同法になると抑えておくとよい。

また、乗算型合同法の周期に関しては下記が成立する。
i) $M$が素数の場合は可能な最長周期は$M-1$である。
ⅱ) $M=2^e (e \geq 4)$を用いる場合は可能な最長の周期は$M/4$となる。

M系列

前述の線形合同法や乗算型合同法は簡易的な手法という意味では優れているが、生成する乱数が$1$つの乱数から決まることから大規模なケースで用いるのが難しいとされている。この話は確率過程を考える際のマルコフ性がシンプルで取り扱いやすい一方で、取り扱う題材によってはシンプル過ぎることで表現力上の問題が生じるのと同様に理解することができると思われる。

そのため、乱数の生成にあたって$2$つ以上の数字を用いる手法の方が望ましいと考えられる。
$$
\large
\begin{align}
a_{n} = a_{n-q}+a_{n-p} \quad (\mathrm{mod} \; 2, \quad q<p)
\end{align}
$$
初期値$a_0, a_1, a_2, …, a_{p-1}$をランダムに定めたのちに（全てが$0$や$1$にならないように気をつける必要あり）、上記のように数列${a_n}$の作成を考える。この数列${a_n}$はM系列と呼ばれる。このときに$X_n$のビット長を$l$とすると、$a_0$から$l$個の数字を取り出して並べたものを$2$進数表記とみなし、これを$10$進数表記に直すことで乱数を得ることができる。$X_1$以降も同様の手順を繰り返し乱数を得る。

解説よりも具体例を元に確認する方がわかりやすいので、「自然科学の統計学」の$11.2.2$節の例に基づいて具体的に確認する。$l=2$、$p=5$、$q=3$、$a_0=1, a_1=1, a_2=0, a_3=0, a_4=1$とする。このとき、$a_5$は下記のように計算できる。
$$
\large
\begin{align}
a_{n} &= a_{n-q}+a_{n-p} \quad (\mathrm{mod} \; 2) \\
a_{5} &= a_{5-3}+a_{5-5} \quad (\mathrm{mod} \; 2) \\
&= a_{2}+a_{0} \quad (\mathrm{mod} \; 2) \\
&= 0+1 \quad (\mathrm{mod} \; 2) \\
&= 1
\end{align}
$$

同様に$a_6, a_7, a_8, a_9$は下記のように計算できる。
$$
\large
\begin{align}
a_{6} &= a_{3}+a_{1} \quad (\mathrm{mod} \; 2) \\
&= 0+1 \quad (\mathrm{mod} \; 2) \\
&= 1
\end{align}
$$
$$
\large
\begin{align}
a_{7} &= a_{4}+a_{2} \quad (\mathrm{mod} \; 2) \\
&= 1+0 \quad (\mathrm{mod} \; 2) \\
&= 1
\end{align}
$$
$$
\large
\begin{align}
a_{8} &= a_{5}+a_{3} \quad (\mathrm{mod} \; 2) \\
&= 1+0 \quad (\mathrm{mod} \; 2) \\
&= 1
\end{align}
$$
$$
\large
\begin{align}
a_{9} &= a_{6}+a_{4} \quad (\mathrm{mod} \; 2) \\
&= 1+1 \quad (\mathrm{mod} \; 2) \\
&= 0
\end{align}
$$

上記より、$2$進数表記で$X_0=11$、$X_1=00$、$X_2=11$、$X_3=11$、$X_4=10$とできる。これを$10$進数に直すと、$X_0=3$、$X_1=0$、$X_2=3$、$X_3=3$、$X_4=2$となる。このように計算することで$2$つの値に基づいて一様乱数を生成することができる。

またこの数列の周期は$2^p-1=2^5-1=31$であり、線形合同法や乗算型合同法に比較して周期が長い乱数となることも着目しておくと良い。

メルセンヌツイスタ法

下記で詳しく取りまとめた。

メルセンヌ・ツイスタ(Mersenne Twister)法まとめ

LCGと多次元疎結晶構造

一様乱数以外の生成

正規乱数の生成（中心極限定理の応用）

中心極限定理を応用することで互いに独立に一様分布$[0,n]$に従う乱数から正規乱数の生成を行うことができる。$n$個の独立な一様乱数の$U_1, U_2, …, U_n$の和$U=U_1+U_2+…+U_n$を下記のように標準化することで正規乱数を得ることができる。
$$
\large
\begin{align}
X = \frac{U-E[U]}{\sqrt{V[U]/n}}
\end{align}
$$

ここで一様分布$[0,n]$に対して、$E[U]$と$V[U]$は下記のように計算できる。
$$
\large
\begin{align}
E[U] &= \int_{0}^{n} x \cdot \frac{1}{n} dx \\
&= \left[ \frac{x^2}{2n} \right]_{0}^{n} \\
&= \frac{n^2}{2n} – \frac{0^2}{2n} \\
&= \frac{n}{2}
\end{align}
$$
$$
\large
\begin{align}
E[U] &= \int_{0}^{n} (x-E[X])^2 \cdot \frac{1}{n} dx \\
&= \left[ \frac{(x-n/2)^2}{3n} \right]_{0}^{n} \\
&= \frac{n^3}{3 \cdot 8n} – \frac{-n^3}{3 \cdot 8n} \\
&= \frac{n^2}{12}
\end{align}
$$

このとき、$X$は下記のようになる。
$$
\large
\begin{align}
X &= \frac{U-E[U]}{\sqrt{V[U]/n}} \\
&= \frac{U-n/2}{\sqrt{n^2/(12n)}} \\
&= \frac{U-n/2}{\sqrt{n/12}}
\end{align}
$$
上記は$X$の確率分布が$\displaystyle n \to \infty$のとき標準正規分布に近づくという事実を利用し(中心極限定理)、正規乱数の作成を行っている。$n=12$がよく用いられる。

正規乱数を作成する他の手法よりもパフォーマンスはよくないとされるが、推測統計学を考える上での中核的な考え方である中心極限定理を実感することもできるので考え方は抑えておくと良いと思われる。

逆関数法

下記で詳しく取り扱った。

・逆関数法の原理と指数分布に基づく乱数の生成
https://www.hello-statisticians.com/explain-terms-cat/inverse-transformation-method.html

ボックス・ミュラー法

下記で詳しく取り扱った。

ボックス・ミュラー法(Box-Muller’s method)を用いた正規乱数の生成

別名法

下記で詳しく取り扱った。

別名法を用いた離散型確率分布に基づく乱数生成の手順とPythonプログラムまとめ