J. Plasma Fusion Res. Vol.96, No.7 (2020)３
７２‐３７８講座 シミュレーションのための乱数入門

４．乱数の品質と検定法
4. Quality and Verification Methods of Random Numbers
佐竹真介
SATAKE Shinsuke
自然科学研究機構 核融合科学研究所（原稿受付：２
０
２
０年３月１
７日）

疑似乱数には様々な発生法がありますが，良い疑似乱数，悪い疑似乱数の判断はどのように考えられているでしょうか．本章ではまず，結晶構造や多次元均等分布性について説明した後，乱数の満たすべき統計的な性質に基づくいくつかの検定テスト法と，それを２．
１節と３章で紹介した疑似乱数の発生法に対して試した結果を紹介します．
Keywords:
pseudo random number, quality of random number, verification of pseudo random number

本章では，「良い乱数」「悪い乱数」とはどういうものません．２．
１節で線形合同法（混合合同法）の問題点は，周か，疑似乱数の品質の判断基準と検定の仕方について概説

期長の短さと，多次元分布における結晶構造であることをします．第２章でも繰り返し述べてきたことですが，疑似触れました．この結晶構造とは，乱数 列 %#)&を用いて，乱数をシミュレーションに利用するにあたって最も気にするべきことの１つが，この乱数の「品質」です．疑似乱数の"次元空間に点列#!
!
#' ）を並べた時
#
#
#
$!""
'!
!"!
!"（
!!
!
!
に，その点列が "次元空間を一様に埋め尽くすのではな

発達の歴史は，いかに品質の良い疑似乱数を効率よく決定

く，ある決まった大きさの "次元格子の頂点にのみ点列が論的なアルゴリズムから生み出すかの探求であると言って存在することを意味します．図１に２．
１節で紹介した線形も過言ではないでしょう．また，疑似乱数に対する品質の不合同法による乱数列 ##"!$&"&#%&#$%&#'*(#$#$!"""!
!!"

満や不信が物理乱数発生器の開発を促したとも言えます．

を３次元空間に並べた場合の例を示します．見る角度を選とは言え，何をもって「良い」疑似乱数であるとするか

ぶと，図のように点列#がある平面上に整列し
#
#
#
!"!
!#$
!!
!
!

というと，これは非常に難しい問題でもあります．例えばていることがわかります．このようなパターンは，真に一偏微分方程式や数値積分の数値解法であれば，より高精度様な乱数であれば出ないため，線形合同法は品質が劣る疑で誤差の蓄積が少ない，保存すべき量を保存する，数値不安定性を起こしにくい，といった明確な良し悪しの基準をその数値アルゴリズムの用途に応じて定められます．しかし乱数の場合，確率論的な統計量に対しての判定の議論が主となり，「この試験をパスした疑似乱数は真の乱数に近い」と言い切れるものがある訳ではありません．むしろ，発生した数列が，真の乱数であるならば満たす，考えうる
限り多くの統計的な性質 A，B，C，…を満たすならばそれは良い乱数だ，という少し漠然とした判断基準になります．また，数学的に厳密な乱数の品質の議論をするには高度な数学（統計学だけでなく整数論，Galois（ガロア）体等）の知識を必要とするため，この講座で全てを説明することはできません．そこで，ここではわかりやすい検定法の実例を挙げながら解説したいと思います．

・疑似乱数の結晶構造と均等分布性

図１「とても良い」乱数を定義することは難しいですが，「とても悪い」疑似乱数を見つけることはそれほど難しくあり

線形合同法 ri+1 ＝ 1103515245 ri + 12345 (Mod 232)で発生した乱数から，３次元空間に点列(ri, ri+1, ri+2)を並べた場合に見られる結晶構造の例．

National Institute for Fusion Science, Toki, GIFU 509-5292, Japancorresponding author’s e-mail:satake@nifs.ac.jp

３７２!2020 The Japan Society of Plasma
Science and Nuclear Fusion ResearchLecture Note

4. Quality and Verification Methods of Random Numbers

S. Satake

"#'
#$"
!
!
(!"(の 各 ビ ッ ト の 選 び 方 に 依 存 し ま#0!

似乱数だと判断されるわけです．
線形合同法については，次元 'をどのように取っても必す．２．
１節 で も 述 べ ま し た が，世 の 中 に 出 回 っ て い る

ず結晶構造を持つことが数学的に証明されています［１］．

Tausworthe 法のプログラムは適切な初期化ルーチンと通対 照 的 に，２．
１節 で 紹 介 し た Tausworthe（ト ー ズ ワ ー

常セットになっているのでユーザがあまり気に病まずともス）法や Mersenne-Twister（メルセンヌ‐ツイスター）法

よいですが，周期長や均等分布性のよい初期値の設定の仕（MT 法）は，長周期性だけでなく，多次元空間における分方について詳しく知りたい方は［３，４］を参考にして下さい．

布の均等性が良いという特徴も持ちます．例えば，最もよところで-ビットベクトルの組+"#""!
"#"(!",が取"#!
0!

!"の３２ビット整数バー

りうるビットの組み合わせは #(- !"$#(!"通りあるの

ジョンの MT 法の場合，６２３次元空間に均等分布するとさで，これらの組み合わせが一周期中に１回ずつ出てくるよ

れています．しかし，実際にこのような超多次元空間の分うな数列を作ることができれば，より長周期でより高次元

布を，長い疑似乱数の一周期に渡って調べることは不可能に均等分布する疑似乱数列を作ることができそうです．そです．そのため，これらの疑似乱数における多次元均等分

こで，Tausworthe 法では各ビットが（１）式から独立に決布性は，乱数発生法の原理から以下のように説明されます

まっていたものを，ビット間でのシャッフルを取り入れる［２］．

ことによってそれを実現したものが MT 法です．MT 法の"(($'く使われている，周期長が #まず Tausworthe 法の場合，漸化式

漸化式は"#"( $"#"(! "#-!* *"#"' $"#"& !'
"#
"#""(
!
-

（１）

（"#は+
!
!
",を元とする-ビットの行ベクトル）に対応する(（２）

でした（式の詳細は２．
１節を参照）．この場合，"#ベクトル)特性多項式 . ". ""が Galois 体 GF（２）上で定義され，こを-ビット乱数として先ほどと同様に漸化式の階数'だけ

れ が +$#(!"の 時 は 多 項 式 .+!"を 割 り 切 る が，並べた組 +"#!
"#""!
"#"'!",を考えると，取りうるビット0!

+"# !"に対しては割り切れない場合に，この特性多項式の組み合わせは #'-!*!"通りになります．MT 法ではこれを (次の原始多項式と呼び，生成される数列の各ビットのらが一周期の中で１通りずつ出るようにパラメータ'
!(
&!
'!
*
!

周期がそれぞれ #(!"となると２．
１節で説明しました．まが上手く選択されています．例えば，２．
１節で紹介したサン((

た，# !"が（１）式の形の漸化式が取りうる最大周期長でプルコード mt１
９９３７ar.f では，'$&#%，&$$('，*$$"，もあります．このように最大周期長を取っている場合，"#-$$#であり，周期長 !$#&#%#$#!$"!"$#"(($'!"になっの $番目のビット .
#!
$に着目して，それを (個ずつの組にし（#$"て並べたもの+!
!
.
.
.
##(!"）を
0!
0!
""!
"(!"!
#
!
$!
#
$!
#
$,

ており，-ビット精度で)
'-!*
#
-*
'
(
$&#$次元均等分布することになります．なお，MT 法では - ビット精度の均等考えると，これらの組の取りうるパターンは，全てのビッ

分布性だけでなく，各上位 ,ビットに対する均等性につい

トが０というパターンを除いた #(!"通りです．実は，特ても工夫されています．上位 ,ビットに関する均等分布性'-!*#,*以下にしかなり得ませんが，MT'(は，理論上 )性多項式が (次の原始多項式となる場合には，この .
#!
$を(個並べた組+が取りうる#(!"通りの全パ
.
.
0!
"(!"!
#
!
$!
#
$,

!
!
-!"のそれぞれについても極力最大値#
0!
法では ,$"

ターンが，乱数列の一周期の中にちょうど１回ずつ含まれに近づくように，"#ベクトルに追加の"Tempering"と呼ばることが証明されています．Tausworthe 乱数の各ビットれる操作を付け加えています（２．
１節の（４）式を参照）．各はこのような意味で (次元に均等分布している，と言われ,ビット精度での均等分布の次元%',(が具体的にどの様にます．

なっているかは，［５］に示されています．このように，MT

ところで Tausworthe 法では - ビットベクトル "#の各

法は単に周期長が長いだけでなく，多次元分布の均等性に

ビットは独立に決まるため，"#の初期値を上手く選ばない

ついても Tausworthe 法より良好になるように作られていと-ビットベクトル全体としての周期を各ビットの最大周ます．

１節で説明しました．で期長#(!"に一致させられないと２．
は - ビットベクトル "#に対する均等分布性に関してはど

・統計的性質による疑似乱数の検定

うでしょうか？その説明をするために，まず定義を明確に

ここまで多次元分布の均等性について見てきましたが，しましょう．疑似乱数の周期を !とした時，- ビットベク

それさえ良ければ品質のよい疑似乱数であるという訳では

トル "#を '個並べた組 +
!
!
!）
#
0!
"#!
"#""!
"#"'!"（,#$"
0!
が ,ビット（,%-）の精度で '次元均等分布するとは，ありません．実際に我々が疑似乱数を使う時には，- ビッ
ト整数乱数 *!
!
!
#- !"/を #- !"で割って .
!
"/区間の
#&.

',

"#!
"#""!
"#"'!",の上位 ,ビットが取りうる # !"個+0!

一 様 乱 数 と し て 使 う こ と が 一 般 的 で す．以 下，（全てゼロは除く）のパターンが，一周期中に同じ回数だけ

の品質につ
.
.
#- !"
(としてこの実数疑似乱数列+
'
#$*
##
#,

現れることを意味します．Tausworthe 法では，"#ベクト

いて調べていきましょう．+.
#,が「良い一様乱数」として

ルを (個（ここで (は周期長ではなく，（１）式の階数）並満たすべき統計的な性質は，いくらでも挙げることができ

べた組 +
"#!
"#""!
"#"(!",が，全 -ビットの精度で均等分0!

ますが，例として以下のようなテストを考えてみましょ布しうる次元は最大で )(#-*であることが知られています

う．

（)"*は "以下の最大整数の意味）．しかし，実際に何次元（!）連検定：例えば，.
!"$.
""".
"#"0".
"& $
#
#".
#
#
#

で 均 等 分 布 す る か は，周 期 長 と 同 様 に 初 期 条 件.
"&""となるような区間のことを，長さ&の上昇連と呼び#
３７３

Journal of Plasma and Fusion Research Vol.96, No.7 July 2020

ます（逆は下降連）．このような連が全体の長さ # の数列

"
"
$
'
0"
今回の検定では，（i）の連検定につい て は /###の中に出現する数を .
/ と表します．あるいは，平均値 1/2と /&(の計６ケースの出現率，（ii）の近接値の出現率に#以上（未満）の値が / 回続く連を考え，その出現数を +/ついても &- #)
!
"
!
"
!
"
!
"
!
"
!
"
!
"
""
"
"#
"
"$
"
"'
"
#"
"
$"
#
""*としたとします．)
4,*が一様乱数であるなら，それらの期待値は場合の計６区間について調べるので，そ れ ぞ れ 自 由 度

以下のように求まります［１，６］．

＝６−１＝５の ($分布に従うかどうかをチェックします．

#
.
*,
#
/ -

一方，（iii）の平均値，二乗平均値の検定は Z 検定によっ

$
#'
/$"%
/"#(
,
/"%(
!
'

て行います．ある確率変数 &の母集団の平均値と分散がの標本平均
4+$(と仮定される時，大きさ #3の標本 )
)""
'
,*

，
/%"%
/$!/!&(!'-

を 4とし，次式で定義される値 Z によって「母集団の平均#"%!/ ．
#
+
*,
#
/ /"#

値が )"である」という仮定が正しいか否かを判定します．

$

十分な長さ #の乱数列を発生させ，その中の長さ/の連の'#出現回数と，上記の期待値との間に統計的に有意な差が生4!)"(#'.3+（５）

じていないかを確認します．

仮 説 が 正 し け れ ば，-は 平 均０，分 散１の 正 規 分 布（!）近接値の出現率：２つの連続した乱数の組 )
4
4
"#*
,"
,

"
'+が 1.96,'"#(に従うはずです．もし標本から評価した +

（0##
"
"
$
/ ）が，区間 &/!##+
4%&/ に現れる回数
"#!4
,
,+

以上の場合，仮説が正しければその確率は５％なので，こ#*/ の期待値は，の場合「有意水準５％で仮説（母平均が )"である）は棄却
された」ことになります．同様に，+%$!
-+')となる場合は

$
$
#
*
#!#(
&/ !&/&/!#!&/!#
*,
$
$
#'
,
'
(
!'
(
/-

有意水準１％で棄却されます．ここでは，疑似乱数 )4,*となります．&/ $#に対して調べることで，近接した値がの平均値，二乗平均値について，母集団が ,""#-の一様乱不自然に連続して出ていないかを確認します．

数であるとの仮説から得られる平均と分散の期待値（"）平均値と二乗平均値の分布：これは最も単純なテス+$(を用いて -検定します．
)""'

トの一つと言えます．一様乱数 )4,*の平均，二乗平均の期

なお通常，-検定は #3個の標本のサンプリング１回につ$待値はそれぞれ *,
4
$
$，*,
##
$
4
%ですが，大きさ##,,

いて -を評価し，仮説が正しいか否かを判定しますが，こ#3の標本を多数取った場合，)
4,，,*が一様乱数であれば 4

こではそのようなサンプリングを %回行い，評価した Z$，の 平 均 値 は そ れ ぞ れ，正 規 分 布 ,'
$
"
$
#3(
#
$
#
#$
(
'
4
,

が実際に,'""#(の正規分布をしているか，５％，１％の棄

$
"
$
#3(
,'
#
%
&
&'
'
(に従って分布するはずです．生成された却域に入るサンプル数が実際に "
!
!
"'%，"
"#%に近いかを

疑似乱数列の標本平均が，この分布に従うかを確認します．

調べます．

なお，連検定と近接値の出現率の検定については，（!）

乱数検定の対象としては，２．
１節で紹介した線形合同法の連の長さ，あるいは（"）の近接値の区間 &/ のインデック，("###"%'#'$&'，，，$#$%$）（4)##$%&'
+1.$(
"4
"##(
,
,'

#
"
"
$
0 に対する出現数分布の期待値が上の *,
，
ス /##
.
/ -

Tausworthe 法 （ !,"1 #!,"1! !,，1#$'"，2##"%），

#
#
，*,
+
*
*,
/ / -のように求められているので，長さ # の標
本中の長さ / の連の出現数の分布が，期待値の分布とどの

Mersenne-Twister 法（mt19937ar.f）と，第３章で紹介した程度よく一致しているかをチェックします．このような「分SFMT については Jump Method によって 2100 飛ばしで作っ布の当てはまりのよさ」のチェックには ($検定が使われます．

た初期値から開始した１０本の乱数列それぞれについてテス

ある確率的事象が -個の互いに背反なクラス !#"
!$"
!- に
0"

０万の標
トしました．それぞれの疑似乱数から大きさ ##１

分別される時，%個の標本中の ,番目の事象 !,の観測度数
0,，その期待度数 ",（" ,#" ,##",#%）に対して次,##0

本を %#１万回サンプリングし，上記の検定を行いました．

KMATH̲RANDOM 版の SFMT の４つとし ま す．なお，

まず，２種類の連と，近接値の出現頻度の分布を見てみ式から求められる値を ($値と呼びます．
$
0,!",('($#" ,##",ましょう．結果はどの疑似乱数もほぼ同じものになったた
め，代表として mt19937ar の例を図２に示しました．この
図では３つのテストそれぞれにおける，連や近接値の実際（３）の出現数（図中点で表示）と期待値（線で表示）を比較しもし観測度数の分布（例えば，疑似乱数列に含まれる長さていますが，一様乱数を仮定して得られる期待値とよい一/ の上昇・下降連の個数の分布）が，疑似乱数が一様乱数

致をしていることがわかります．この出現率の分布の期待であると仮定して求めた期待値の分布に実際に従っている

値との適合度を見るために，（３）式で評価した ($値の度なら，（３）式に従って長さ#の標本を%個用意して評価し数分布をグラフにしてみると，図３のようになり，確かにた ($値は，自由度 *#-!#の ($分布に従うと期待されま

自由度５の ($分布に従っていることがわかります．紙面のす．ここで自由度 *の ($分布（確率密度関数）は次式で定都合上全ての結果について図は載せませんが，SFMT の１
０

義されます．

本の乱数列だけでなく，Tausworthe 法や，乱数としての#$!# !($$
*
$
+
("
($(
*(
/ $
'
'
# *$$
*
$
$ ''
$(
$

品質が劣るとされている線形合同法についても図２や図３（４）をプロットしても有意な差は見られませんでした．
次に疑似乱数の平均値，２乗平均値の分布について見て
３７４

Lecture Note

4. Quality and Verification Methods of Random Numbers

図２

mt19937ar の疑似乱数に対する２種類の連と近接値出現回
数 の 分 布．“up-down”が 上 昇‐下 降 連，“＞ ＜ 0.5”が 0.5以上又は以下の連，が #m−1 ＜ |xi+1 - xi|! #M の近接値の出現数であり，点がそれぞれの観測値，線が期待値を表します．長さ１０万の乱数列を１万回サンプリングして評価しました．

図３

mt19937ar の疑似乱数に対する２種類の連と近接値出現回数の出現数分布に対する $2 値の度数分布（サンプル数＝１
万）．点がそれぞれの観測値，曲線が（４）式による自由度
５の $2 分布を表します．

みましょう．図４に線形合同法の結果を代表として示します．こちらも，一様乱数を仮定して求めた正規分布とよい

S. Satake

図４

線形合同法で発生させた疑似乱数（標本の大きさ１０万）の平均値と２乗平均値の分布（サンプル数１万）．なお２乗
平均値は 0.16 右にオフセットしてプロットしてあります．
“Gaussian”は一様乱数の場合に期待される正規分布を表
します．

図５

平均値，２乗平均値に対する Z 検定と，２種類の連と近接
値 の 出 現 率 に 対 す る $2 検 定 で，そ れ ぞ れ 有 意 水 準
５％，１％で棄却された標本数．横軸は乱数の発生法を表
０個 の SFMT，１
１が 線 形 合 同 法，１
２がし，ID ＝ 1〜10 が１
Tausworthe 法，１３が mt19937ar の結果です．また，ID ＝ 0
は１０個の SFMT の平均値を表します．

一致を示していることがわかります．この場合も，乱数のきる疑似乱数の発生法が「より乱数らしい」と言えるだけ発生法による有意な差は現れませんでした．（５）式によるです．また，様々な検定を組み合わせても，「実際上この(検定の結果についても見てみましょう．図５では，平均種の検定で棄却できるのは異常に不都合なアルゴリズム

値と二乗平均値について，有意水準５％，１％で棄却されで，非常によいものを積極的に見出すことは，はなはだ難るサンプルが１万個のうち実際にいくつあったかをプロッしい」（［１］の p．
３２より抜粋）ことが以前から指摘されてお#トしてあります．また，図３で見た $ 値の分布についてり，それゆえにより複雑な乱数検定法が多く考案され続けも，$#値が大きくて連などの分布が期待値によく一致してています．例えば，NIST SP 800-22［７］と呼ばれる乱数検定いると言えない標本数（自由度５の場合 $#"""!
!%が５％，

用のツール群には１５種類もの検定法が含まれています．し

#

$ ""$
!
!'が１％の棄却域）も合わせて示しました．どのかしながら，その検定法の中には問題点を指摘されている

乱数発生法のどの検定結果も，実際にほぼ５％，１％の標ものも存在しているため［８］，利用には注意が必要です．

本がそれぞれの棄却域に入っています．

なお，NIST の検定法では疑似乱数に対する各テストの合

これらの結果から，ここで取り上げた４つの疑似乱数

否判定は，棄却域に入るサンプル数に対する Z 検定によっては，どれも今回紹介した検定をパスしたと言えそうです．

行っており，具体的には各テストにおいて１％棄却域に入った

驚くべきことに，線形合同法は明らかに結晶構造という乱数らしからぬ性質を持っているのに，これらの検定はパス

!
!
!"$(
標本数が全体数 ! に対して !"!
#できました．２．
１節で整数の線形合同法乱数が偶・奇，偶・（(##!
$&）の範囲に入っていれば，疑似乱数が「期待通り

奇，のパターンを繰り返すことを紹介しましたが，このよに有意水準１％の棄却域でテストに不合格した」ため「乱うな乱数らしからぬ性質も今回のテストでは見逃されてし数検定のテストに合格した」という判断をします．ここでまいます．本章の冒頭でも述べた通り，「このテストさえも，乱数の品質の検定はあくまで統計学的なスタンスでな

合格すれば OK」と言い切れるような，万能な乱数の検定

されることが現れています．

法がある訳ではなく，考えうる限り多くのテストを合格で

!
!
!
!"%
"!!
!"&
，
"
!

NIST の合否判定を今回参考として行ったテストに当て
３７５

Journal of Plasma and Fusion Research Vol.96, No.7 July 2020

はめると，#"１万なので合格の範囲は100±25ということ

化 MT 法として紹介した SFMT において改善されていまになります．テスト結果の内訳を見ると，１０個の SFMT

す．もし，BigCrush テストの結果を重視するので あ れ#の内の２つと線形合同法の上昇・下降連テストの " 検定ば，MT 法を改良して省メモリ化し，かつ BigCrush テスト

において１２６あるいは１２７個の標本が１％の棄却域に入ったも全て通るように出力を工夫した TinyMT［１０］が MT 法

他は，全て 100±25 の範囲に収まりました．上昇連，下降の作者らによって発表されているので，それを使うという連のテストについては筆者自身の経験として，一標本当た

選択肢もあります．ただし均等分布性の次元は元の MTりの長さ "をかなり大きく取らないと，乱数の発生法に依法より低くなっています．TinyMT はまだ新しく一般的に##らずテストから評価した " 値が期待される " 分布より若普及していないため本講座では紹介しませんでしたが，干上振れする傾向があります．これは長さ "の数列に含まDynamic Creator や Jump Method も備えているので，最新れる連の総数自体が確率変数であることや，図２からわかの疑似乱数を試したい方は利用してみてはいかがでしょう

るように平均値以上・以下の連の分布と比較して，上昇・か．

下降連の分布は長さ５や６以上の連の出現率が長さ１，２の連に比べて極端に少ないため，"#値の誤差（分散）が大・並列乱数の品質について

きく出やすいためではないかと予想しています．いずれにここまで，１つの疑似乱数発生法に対する品質の検定にしても，棄却されたサンプル数の１，２程度の差をもってついて説明してきましたが，最後に，第３章で紹介したよ

直ちに SFMT の内の２つと線形合同法がテストに不合格

うな並列疑似乱数についての品質について少し述べたいとだとみなすのは早計でしょう．実際，１サンプルの長さを思います．乱数を並列化する上で重要な性質は，それらの１０万から１００万に延ばすと，全ての &検定，"#検定において乱数列同士が互いに独立で相関関係を持たないことです．

１％の棄却域に入った標本数が 100±25 の合格圏に収まり第３章では MT 法の並列化の方法として，過去に本学会誌ました．このように，検定法によって適切な１サンプルの［１１］で紹介した Dynamic Creator の他に，１本の MT 法乱

大きさが異なったり，"#値の出方に癖があったりするの数列のそれぞれ遠く離れた点を出発点として，並列計算にで，NIST やその他，世の中に出回っている乱数の検定法

おいて異なる部分区間の乱数列を利用するという Jump

を使う際，テストのサイズ設定や合否判定の解釈には注意

Method を新たに取り上げました．並列計算によるシミュが必要です．

レーションで，これら並列疑似乱数を利用する時には，そ検定プログラム群としては他にも，TESTU01［９］と呼ばれらの乱数列が互いに独立であるという前提で利用するわれる乱数の品質検定プログラム群が知られています．こちけですが，その保証はどこにあるのでしょうか？

らはテストの規模に応じて SmallCrush，Crush，BigCrush

Dynamic Creator の場合は，MT 法を特徴づける漸化式

いう３つの検定プログラム群が用意されており，最も詳細の形が異なるものを必要な数だけ用意するというものでしな検定をする BigCrush では１０６種類もの検定法による１
６０

た．その独立性については，作者自身も厳密な数学的証明

通りのテストが行われます．文献［９］の Table.I には多くのはなく，漸化式を特徴づける特性方程式が互いに素な疑似乱数発生法に対する検定テストの結果が載っており，棄却

乱数同士は独立である，という仮説に基づいていると説明

域 10−8％ で不合格だったテストの数が示されています．残

しています［１２］．一方，Jump Method では１本の MT 法の念ながら本節で試したものと同じ線形合同法はありません

疑似乱数列を例えば 2100 飛びに使うという方法であり，そが，パラメータの異なる様々な線形合同法の結果が載ってれぞれの部分区間が統計学的な意味で独立とみなせるかど

おり，おしなべて多くのテストで不合格になっています．

うかについて，筆者が調べた限り具体的な説明は見当たり

Tausworthe 法（本節と同じものは文献［９］の Table.I 中のませんでした．しかし，MT 法が持つ一般的な意味での疑GSFR（２５０，１０３）），MT 法（同 Table 中の MT19937）の結

似乱数の品質の良さや高次元均等分布性を踏まえると，部果も載っています．Tausworthe 法に比べて MT 法の方が分区間 A と B から取った数列 %%$&&と %&&に対して，例え不合格になったテストの数が少ないことがわかりますが，

ば $&%
（'"!
!
!
"
' $に対する相関係数を具体的に評価し!'&MT 法でも BigCrush の１
６０個のテスト中２つで不合格にたところで，何かおかしな相関が現れることはないと期待なっています．これを「たった２つ」とみなすか，「２つもして良さそうです．このように，MT 法による疑似乱数並不合格」とみなすかは，先ほども述べたように万能な乱数列化に関しては，まず間違いなく大丈夫だろう，という感のテストというものがないため難しいところです．なおじで実応用に供されているというのが実情のように思われMT 法については初期のバージョンには，初期化の seed

ます．

の値を変えても初期ベクトルの状態があまり変わらない，

筆者自身は１
５年くらい前から MT 法（当時は Dynamic

一旦漸化式中のベクトルのビットが０が過多な状態に陥る

Creator，現在は Jump Method）を並列シミュレーションとそこからなかなか抜け出せない，といった問題点が指摘に応用してきましたが，当時から並列乱数列の独立性につされていました．これらが２つのテストで不合格になった

いて議論や検証の話が見当たらないことに不安を感じつつことと関係があるかについては，筆者が調べた限りでは

使 い 始 め ま し た．最 終 的 に は，１６０本（後 に１０２４本）のはっきり述べられたものはありませんでした．初期化の問Dynamic Creator の MT 法乱数列について， "%!!##"!#$!#$

題については本講座で紹介した mt19937ar において既に改通りの全ての組み合わせに対して，相関チェックを自ら試善されており，また０超過状態からの回復についても並列すという実力行使に出ました．この際に考えた Checker３７６

Lecture Note

4. Quality and Verification Methods of Random Numbers

S. Satakeboard test と 名 付 け た 検 定 法［１３］で は，２つ の 乱 数 列

す．空間セルの数が例えば１０万個ある場合，統計的には

$
!
#&,
!
#(.に対し，)#(##(*の２次元のマス目を考+(+&,

１０００個のセルで疑似乱数が１％の棄却域に落ちる可能性が

え，)
#&!
$
&*を順次打っていき，i ステップ目まで進んだ時ある，ということです．

に，全部で "$##( 個のマス目のうちまだ１度も点を打た

この例のように，初期分布を作るのに乱数を使う場合はが無相関な一様
れていない数')$%#&,+*を観測します．+&,

１回きりのことなので，検定を掛けて棄却域に入った場合

乱数だと仮定した場合の期待値とその分散はは乱数を作りなおしたとしても，総計算時間に対する計算時間の増加は無視できるでしょう．また，テスト粒子の平

%

"!"
%
')
%
&"!!
*
.
$"#
"$'"'*)!! "，""均初期速度がなるべく仮定した初期条件通りの状態から計算をスタートさせたいといった要請があるかも知れませ

#

"!#
"')
%
"")
"!"*
*
.
$"#
!
"
"
#

ん．しかしそれならば，平均値が棄却域に入るような乱数でもそのまま "個のテスト粒子の初速の決定に使い，後で

%

"!"
"!"*
#
!)
!
"$
"

粒子の速度を全体的にずらして仮定した初期条件に粒子の平均速度を合わせてもよいわけです．また，時間発展シと 与 え ら れ ま す．こ れ に 対 し，観 測 し た ')%*がミュレーションの中で，ランダムウォークや衝突計算を乱

!%
%
()
"&$*
*
.
%)
#")
*の範囲に現れる割合が，上式の期数で模擬する場合，１％の棄却域に入る乱数を一々排除し待値と分散を持つ正規分布に従うかどうかをチェックしまていると計算時間の増大につながりかねませんし，何よりした．もしテストの結果がそうなっていなければ +
$
#&,
+
&,

本当に１％の確率で起こるはずの「レアイベント」を人為が無相関な一様乱数だという仮定が棄却される，即ち相関的に排除してしまう可能性があります．

があると考えられるわけです．MT 法は３２ビットの精度を

筆者自身は初期化においても時間発展計算の最中におい持っていますが，その精度でチェックすることは時間的にても，疑似乱数の検定はしてもいいが，それは乱数の品質不可能だったので，+
$
!
#&,
!
#"%.にビットを落として+(&,の確認のためだけであり，棄却域に入る乱数もそのまま使%$#!
" まで追跡しました．その結果は［１３］に示すよううというスタンスを取りますが，これはシミュレーションに，独立な乱数として期待される通りに ')%*が%と共に減コードが扱う問題の性質や研究者個々人の乱数に対する考少するということが，Dynamic Creator で作った全てのえ方によって変わりうるものだと思います．ここで強調しMT 法の乱数列の組み合わせについて確認できました．

ておきたいことは，乱数の検定法が統計的な性質に基づく

このようなテストをわざわざするほど並列乱数の独立性ものである以上，真の乱数でも疑似乱数でも，ある確率でについて心配する必要は今にして思えばなかったかも知れ

検定の棄却域に入る乱数の標本は必ず生じ，むしろそれはません．また，このテストに限らず乱数の検定全般につい乱数として自然な振る舞いであるということです．

てよく指摘される点ですが，いくら乱数のテストを行ったとしてもそれはテストを行った部分区間における品質を

本章では筆者自身の過去の経験を交えながら，疑似乱数

チェックしたに過ぎないということは否めません．しかしの品質について概説してみました．乱数を使うシミュレー

ながら，並列乱数の独立性については自らテストをしなけ

ションを行う研究者の方々が，ふと「自分が使っているこ

れば何も分からないという程に，検定法や検定結果が発表の疑似乱数は大丈夫なのだろうか？」という疑問を感じた

されていないというのが実情です．並列計算に疑似乱数を時にこの記事で紹介した内容が参考になれば幸いです．

応用するプログラマは，その事実を頭の片隅に留めておいてもよいのではないか，と私は考えます．よりシビアに乱

参考文献

数の独立性を求めたいユーザは，第２．
２節で紹介したよう［１］津田孝夫：モンテカルロ法とシミュレーション（三訂版）（培風館，１９９５）．
［２］松本 眞，栗田良春：Twisted GFSR：新 し い 乱 数 発 生
法，京都大学数理解析研究所講究録 85, 86-95 (1993).
［３］伏見正則：乱数（東京大学出版会，１
９８９）．
［４］B.J. Collins, G. Barry Hembree, J. Association for Computing Machinery 33, 706-711 (1986).
［５］松 本 眞：有 限 体 の 擬 似 乱 数 へ の 応 用，http://www.
math.sci. hiroshima-u. ac. jp / 〜 m-mat / TOKYOHOMEPAGE/TEACH/1011.pdf［６］H. Levene and J. Wolfowitz, Ann. Math. Stat. 15, 58-69(1944).
［７］https://csrc.nist.gov/projects/random-bit-generation/
documentation-and-software［８］吉 田 等 明 等：電 子 情 報 通 信 学 会 信 学 技 報 112, 13(2012).
［９］P. L'ecuyer and R. Simard, ACM T. Math. Software, 33,
Article 22 (2007).

な物理乱数の利用を検討するのも一つの手段だと思いま
す．

・おわりに
最後に本章の補足として，実際にシミュレーションを行う際に乱数を検定してから使うべきか否か，という問題について少し述べたいと思います．例えば粒子シミュレー
ションで１空間セル当たり " 個のテスト粒子の初期速度を乱数で与える際に，大きさ "の疑似乱数に対して本節で例を示したような統計学的な検定テストを行え
ば，５％，１％など任意の有意水準で「一様乱数らしくない」疑似乱数のセットを排除することも可能です．しかし，疑似乱数が実際に真の乱数に近い挙動をしている場
合，１％の棄却域に入るような " 個の疑似乱数のセットは，実際に約１％の確率で生じて然るべきものでもありま
３７７

Journal of Plasma and Fusion Research Vol.96, No.7 July 2020

［１０］M. Saito et al., RFC 8682, DOI 10.17487/RFC8682 (https://www.rfc-editor.org/info/rfc8682).
［１１］講座「核融合プラズマシミュレーションの技法 ５．粒子シミュレーションのコーディング技法」プラズマ・核融合学会誌 89, 245 (2013).
［１２］M. Matsumoto and T. Nishimura, "Dynamic Creation of
Pseudorandom N umber Generators", Monte-Carlo andQuasi-Monte Carlo Methods pp. 56-69 (Springer, Berlin,
Heidelberg, 1998).
［１３］S. Satake et al., Lecture Notes in Computer Science bookseries Vol. 4759, pp. 344-357 (Springer-Verlag, Berlin,
Heidelberg, 2006).

なお、本講座の執筆のために統計数理研究所の研究者に意見を聞きに伺ったことがきっかけで、朝日新聞の記者から乱数のシミュレーション研究への応用について取材を受けました。その記事は２０２０年６月３日の朝日新聞夕刊および朝日新聞デジタルに『現場へ！「乱数の世界へようこそ」』というタイトルで掲載されました。

!!!!!!!!!!!!!!!!

さ

たけ

しん

すけ

佐 竹 真 介自然科学研究機構 核融合科学研究所 ヘリ
カル研究部 核融合理論シミュレーション
研究系 准教授，２００
３年 総合研究大学院大学 博士（学術）．
モンテカルロ法を使った３次元磁場配位中の新古典輸送現象，新古典粘性のシミュレーションや最適化配位の研究が主なテーマ．乱数については深い思い入れがありますが，特にギャンブル好きというわけではありません．

!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!

３７８