乱数の生成アルゴリズムとその応用|問題演習で理解する統計学【8】

下記などで取り扱った、乱数の生成に関する問題演習を通した理解ができるように問題・解答・解説をそれぞれ作成しました。
https://www.hello-statisticians.com/explain-terms-cat/random_sampling1.html

基本問題

線形合同法・乗算型合同法

・問題
線形合同法(linear congruential method)は$1948$年頃にレーマー(Lehmer)が提案した手法で、仕組みがシンプルであるので乱数生成のイメージを掴むにあたって適した手法である。また、乗算型合同法は線形合同法の特殊な場合と考えられるので、両者は同時に理解することができる。
さて、線形合同法を理解するにあたっては、「整数に関する合同」と「数列の漸化式」の理解があれば十分であると思われる。以下の問題に答えよ。
i) $x_{n+1} = x_{n} + 2, x_1 = 1$のように表せるとき、$x_2, x_3$を求めよ。
ⅱ) $x_{n+1} = 2x_{n} + 1, x_1 = 1$のように表せるとき、$x_2, x_3$を求めよ。
ⅲ) 合同は「整数$x$を整数$M$で割った時の余りが等しいか」を表す概念である。このときの$M$を法(modulus)という。たとえば$2$と$5$は$3$で割った際にどちらも余りが$2$であるので、$3$を法とするとき合同であるといえる。このことは「$2 = 3 \quad (\mathrm{mod} \, 3)$」のように表すことができる。
上記に基づいて、法を$5$とするときに$2$と合同な自然数を$5$つ答えよ。
iv) $9$を法とするとき、$x_{n+1}=4x_{n}+2 \quad (\mathrm{mod} \, 9), x_1=1$に従って$x_2, x_3$を求めよ。
v) iv)の例の$x_4$以降を計算すると下記のようになる。
$$
\begin{align}
x_4=7, x_5=3, x_6=5, x_7=4, x_8=0, x_9=2, x_{10}=1, x_{11}=6, x_{12}=8
\end{align}
$$
上記には周期があると考えることができるが、どのような周期か答えよ。

・解答
i)
$x_2, x_3$は下記のように求めることができる。
$$
\begin{align}
x_2 &= x_{1} + 2 \\
&= 1 + 2 \\
&= 3 \\
x_3 &= x_{2} + 2 \\
&= 3 + 2 \\
&= 5
\end{align}
$$

ⅱ)
$x_2, x_3$は下記のように求めることができる。
$$
\begin{align}
x_2 &= 2x_{1} + 1 \\
&= 2 + 1 \\
&= 3 \\
x_3 &= 2x_{2} + 1 \\
&= 6 + 1 \\
&= 7
\end{align}
$$

ⅲ)
$$
\begin{align}
7, 12, 17, 22, 27
\end{align}
$$

iv)
$$
\begin{align}
x_2 &= 4x_{1} + 2 \\
&= 4 + 2 \\
&= 6 \\
x_3 &= 4x_{2} + 2 \\
&= 24 + 2 \\
&= 26 \\
&= 8 \quad (\mathrm{mod} \, 9)
\end{align}
$$

v)
$x_1$〜$x_9$と$x_{10}$〜$x_{18}$は一致し、これを周期とみなすことができる。線形合同法の計算の仕組み上、$1$度でも同じ値が出てくると以後の値は周期的に繰り返す。

・解説
線形合同法は「整数の合同」と「数列の漸化式」について理解していればそれほど難しくない手法だと思います。乗算合同法は線形合同法の定数項がないパターンであり、重複した値が得られやすい一方で計算コスト自体は低いため、合わせて抑えておくと良いと思います。
線形合同法、乗算型合同法の数式は下記でまとめたので、詳しくは下記も参照してみてください。
https://www.hello-statisticians.com/explain-terms-cat/random_sampling1.html

M系列

正規分布に従う乱数の生成

・問題
一様分布に従う乱数が得られる前提で正規分布に従う乱数の生成を行うにあたっては主に「中心極限定理」を用いる手法と「ボックスミュラー法」を用いる手法がある。中心極限定理に基づく手法は複数の乱数の平均が中心極限定理に基づいて正規分布に従うことに基づく手法であるが、$1$つのサンプルを作成するにあたって$n=12$個がよく用いられるなどあまり効率的ではない。よって以下では「正規分布に従う乱数の生成」にあたって実用的に用いられることの多い「ボックスミュラー法」の原理に関して確認を行う。

一様乱数$U_1 \sim \mathrm{Uniform}(0,1), U_2 \sim \mathrm{Uniform}(0,1)$に対し、下記の変換によって得られる$X, Y$は互いに独立に標準正規分布$\mathcal{N}(0,1)$に従う。
$$
\begin{align}
X &= \sqrt{-2 \log{U_{1}}} \cdot \cos{(2 \pi U_{2})} \quad (1) \\
Y &= \sqrt{-2 \log{U_{1}}} \cdot \sin{(2 \pi U_{2})} \quad (2)
\end{align}
$$

上記に基づいて乱数を生成する手法を「ボックスミュラー法」という。ここまでの内容を元に下記の問いに答えよ。
i) 関数電卓やPythonなどを用いて$(U_1,U_2) = (0.1,0.3), (0.8,0.6), (0.2,0.9)$に対応する$(X,Y)$の数値の計算を行え。

・解答
i)
下記を実行することでそれぞれ計算を行うことができる。

import numpy as np

U_1, U_2 = np.array([0.1, 0.8, 0.2]), np.array([0.3, 0.6, 0.9])
X = np.sqrt(-2*np.log(U_1)) * np.cos(2*np.pi*U_2)
Y = np.sqrt(-2*np.log(U_1)) * np.sin(2*np.pi*U_2)

print(X)
print(Y)

・実行結果

> print(X)
[-0.66313997 -0.54046156  1.45147566]
> print(Y)
[ 2.04093497 -0.39266831 -1.05455879]

発展問題

モンテカルロ積分

MCMC

参考書籍

・自然科学の統計学(東京大学出版)