J. Plasma Fusion Res. Vol.96, No.6 (2020)２
９０‐２９１

講座

シミュレーションのための乱数入門
Introduction to Random Numbers for Simulations

１．はじめに
1. Introduction
宇佐見俊介
USAMI Shunsuke
自然科学研究機構 核融合科学研究所（原稿受付：２
０
２
０年３月１
７日）

非平衡系であるプラズマ中の現象を解析するには，コン（a）０，１，０，１，０，１，０，１，０，１，０，１，０，１，０，１

ピュータによるミュレーションが非常に有力な手段となっ（b）１，０，１，１，１，０，１，０，０，０，０，１，０，０，１，０ています．シミュレーションには様々なモデルや計算手法がありますが［１］，そのうち，いわゆるモンテカルロ法にの２種類を挙げてみます．
（a）は，明らかに０と１が順番基づいた輸送シミュレーションでは，乱数を用いた計算をに出現していて，その規則性ゆえに乱数には見えません．

繰り返すことから，大量の乱数データを必要とします．他

一方，（b）は，私自身がコインを投げて表を１，裏を０とにも粒子シミュレーションでは，粒子の初期位置や速度分して記録した結果ですので，この数列は乱数に近いと言え

布（Maxwell 分布など）を与える際に，乱数が用いられるるでしょう（＊）．一見，（b）では，同じ値が連続することのが一般的です．また，乱数は，これらのような科学技術が頻繁にあって作為的な感じがしますが，乱数にはこのよ

計算だけに留まらず，インターネットショップでの情報通

うな並び方が起こり得ます．０と１が規則性なく出現する

信の暗号化，二段階認証におけるワンタイムパスワード生ということから，同じ値が連続で出現する率はどのくらい成などにも用いられ［２］，現代社会で必要不可欠なツール

かという統計的性質が導けます．そのため，逆に言えば一となっています．

見乱数に見える数列も，乱数が満たすべき統計的性質から本講座では，そのような乱数について学んでいきます．

外れていると，乱数とは言えないでしょう．ただし，厳密

講座を読み始めるにあたって，最初に思い浮かぶ疑問は，には，無限の数列がないと，ある数列が乱数かどうかは判「そもそも乱数とは何か？」ではないでしょうか．乱数が，定できません．実は（a）の場合も，無限に続く乱数列中に含どのように定義されているのか調べてみると，大辞林（第まれうることに留意する必要があります．ならば，科学技三版）では，「乱数とは出現する値に規則性のない数」と術計算で（a）
も乱数として使ってよいのか，という疑問があります．この「規則性がない」という定義から，乱数が出てきますが，答えは NO です．科学技術計算では，有限満たすべき性質が導かれるように思いますが，数学的に厳の数列しか利用しないので，その有限の中で乱数としてふ

密に定義しようとすると簡単ではありません［３］．そこで，

さわしい統計的性質を満たしていることが望ましい，つままずは，以下のような例を用いて，乱数の性質の一端を覗り（b）が望ましいということになります．乱数発生法においてみましょう．０と１だけからなる数列として，

いて，この点は，議論のあるところではないでしょうか．

（＊）コインの動きは古典力学に従いますので，表裏どちらが出るのかは決定論的であり，ランダムではないという反論もあるかもしれません．一方，コインの動きは，条件のわずかな違い（コインを投げる速度，机の凸凹など）に非常に敏感なので，結果は事実上予測できません．このことから，表裏どちらが出るのかはランダムと言ってよいでしょう．
National Institute for Fusion Science, Toki, GIFU 509-5292, Japanauthor’s e-mail: usami.shunsuke@nifs.ac.jp

２９０!2020 The Japan Society of Plasma
Science and Nuclear Fusion ResearchLecture Note

1. Introduction

上記で述べたような性質は，コンピュータで乱数を用い

S. Usami（第２回目）

る際に，極めて重要な要素です．一般的なコンピュータは第３章では，第２章で紹介した様々な乱数発生法を高速

真の乱数を作れず，様々なアルゴリズムにより「疑似乱数」

化し，並列実行のプログラム中で用いる方法について説明を発生させています．コンピュータ（すなわち人間）が作します．

り出す疑似乱数の値は，一見ランダムに見えますが，決定第４章では，乱数の品質と検定法を解説します．乱数ら論的な演算（数式）によって算出されたものですので，原しさについて，周期長や多次元均等分布性といった数学的理的に周期性や，偏りを持っています．そのため，乱数のな観点から説明し，乱数が満たすべき統計的性質を利用し性質を満たしているかどうかが大きなチェック対象となりて，乱数の検定方法をいくつか紹介します．

ます．
コンピュータシミュレーション分野では，黎明期から現（第３回目）

在まで，より長周期で，偏りが少なく，かつ高速に疑似乱

最後に第５章では，乱数の変換例・利用例を示します．

数を生成できるアルゴリズムが開発されてきました［４］．

まずは，様々な乱数をいかに生成するか，つまり「乱数の現在では，一般のユーザはプログラミング言語に標準で備変換」について論じ，一様乱数から任意の確率分布に従う

えられている，あるいは数値ライブラリで提供される乱数乱数を作る方法を説明します．また，基礎的なモンテカル

発生ルーチンをブラックボックスのように利用することがロシミュレーションに加えて，これまであまり焦点の当できます．しかし，その疑似乱数がどのように生成されてたっていなかった，輸送方程式（非線形偏微分方程式）のいて，その品質がどのようなものか，また間違った利用方

解を求める過程を大域的最適化問題に帰着させるという解法をしていないか，ということを一度振り返って考えてみ法を紹介し，乱数がどのようにシミュレーションにおいてることは非常に重要なことです．

利用されているか具体的に解説します．

本講座では，実際に自ら，乱数を使ったシミュレーショ
ンモデルを開発して実行している研究者が，乱数が満たす

本講座が，これからシミュレーションで乱数を用いよう

べき条件，乱数発生方法のいくつかのバリエーション，検とする研究者，学生，さらには，乱数自体に興味のある定方法，高速化，応用例について，詳しく説明します．た人々にとって入門的な役割を果たせたなら幸いです．

だし，著者らは乱数そのものの専門家ではないことから，数学的に厳密であることを追求しすぎることは避け，シ

参考文献

ミュレーション研究者としての実践的な説明をすることを［１］プラズマ・核融合学会（編）：プラズマシミュレーショ
ン（京都大学学術出版会，２０１８）．
［２］「ランダムと乱数の奇妙な世界」Newton 2018 年 9 月号
pp.70-83.
［３］D.E. Knuth, The Art of Computer Programing Volume 2 Seminumerical Algorithms, Third Edition (Addison-Wesley,1997).
［４］伏見正則，逆瀬川浩孝（監訳）：モンテカルロ法ハンド
ブック（朝倉書店，２
０
１
４）
．

心がけました．
掲載内容は，以下のように予定しています．
（第１回目）
第２章では，コンピュータにおいて，疑似乱数列を生成する様々な手法を主に周期長という観点から紹介し，実際の使い方も説明します．さらに，コンピュータで真の乱数を発生させる「物理乱数」に関しても触れることにします．

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

う

さ

み

しゅん すけ宇 佐 見 俊 介
自然科学研究機構 核融合科学研究所 ヘリ
カル研究部基礎物理シミュレーション研究
系 准教授．
粒子シミュレーションを用いて，プラズマの様々な複雑な現象を調べていて，現在の主な研究テーマは
磁気リコネクションです．SF が好きなためか，現実世界はコンピュータ上のシミュレーションでは？と密かに妄想しています．今回の講座をとりまとめている際に，世界がシミュ
レーションなら，量子ゆらぎに基づく乱数も実は疑似乱数
で，何らかの規則性があったりしないか？とさらに妄想を膨らませています．

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

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

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

２９１

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

２．乱数発生の原理
2. Random Number Generation
佐 竹 真 介，菅 野 龍 太 郎
SATAKE Shinsuke and KANNO Ryutaro
自然科学研究機構 核融合科学研究所（原稿受付：２
０
２
０年３月１
７日）

乱数は，発生法により大きく２つに分けることができ，１つは，決定論的な演算による疑似乱数，もう１つは，ランダムな自然現象を乱数源として利用する物理乱数です．乱数を利用するシミュレーションは，当然のことながら利用した乱数に，その計算結果が影響されます．そのため，利用する乱数について，事前に理解しておくことが大切です．本章では，まずは，疑似乱数，物理乱数それぞれの発生法について概説します．
Keywords:
pseudo random number, physical random number, quantum random number generator

２．
１ 疑似乱数の発生法

では我々がコンピュータプログラムで利用している乱数

コンピュータシミュレーションでは，様々な用途で乱数とは何なのでしょうか？本節ではまず，我々が通常プログが必要となります．必要な数が１０個程度であれば，人の手ラムで利用している「疑似乱数」の発生原理について説明

でサイコロを振るなり乱数表を引くなりすればよいです

します．疑似乱数とは決定論的なアルゴリズムから生成さが，乱数を使うシミュレーションで要求される乱数の数はれる「乱数のように見える」数列のことです．本節ではい

一般的に数万から数億個といった膨大な量になります．こくつかの代表的な疑似乱数の生成方法と，それらの利点・

れは乱数を利用するシミュレーションの多くが，ある種の欠点を紹介します．なお本節の内容は［１‐３］を参考に書か

ランダムな事象をプログラムで模擬し，その統計的な性質れています．一方２．
２節では，まだあまり一般的に普及してに基づいて解を求めるという考え方，言い換えれば，あるいませんが，疑似乱数とは異なり現実の確率的事象の観測

種のアンサンブル平均を数値計算で評価することによって値から乱数列を取得する「物理乱数」について紹介します．

解を得るという考え方に基づいて作られているからです．

また，疑似乱数の発生法の高速化・並列化といった数値的

一般的に，確率変数のアンサンブル平均の統計誤差は標本な技法や，疑似乱数の品質，つまり「乱数らしさ」をどう

数 ! に対し "#!
( に比例して小さくなります．しかし，も定義し，評価するかについてはそれぞれ第３章，第４章でし乱数と思って使っている数列に予期せぬ偏りや相関があ

詳しく解説します．

る場合，それを使っていくら多くの標本を用意したとして!線形合同法も，得られる答え（統計量）の信頼性は上がりません．よって，性質の良い乱数を効率よく発生させることは，乱数をさて，決定論的な計算機の演算から乱数のように見える

使うシミュレーションの規模が大きくなるほど，またより数列を生み出す方法ですが，一般的に使われる方法は全て高い精度の計算が要求されるほど重要になると言えるで何らかの漸化式の計算によって生成します．また，コン

しょう．

ピュータの内部では数字は２進数で表現されますから，こところで，我々が現在利用できるコンピュータは，内部

こで説明する疑似乱数とは４バイトや８バイトの２進数で的にはデジタルに表現される数値の四則演算と論理演算表現された整数だと考えて下さい．単純で古くからよく知（AND，OR など），ビット操作（２進数で表現された数値

ら れ た 疑 似 乱 数 は「線 形 合 同 法」と 呼 ば れ る も の で，を全体に右または左にシフトする，特定のビットの値を変

: %#"
!
!
#
) を疑似乱数列とした時，一般的に以下のよ'&%'更するなど）を組み合わせた手続きによって動いており，

うな漸化式で表現されるものを指します．

完全に決定論的な演算しか行えません．一方で，乱数とい

'
$&%"%
$
!'
"'
""##
!""*"#
!&"$
%
%"#
%
&'
%

うのはその名の通りランダムで予測不可能な数列というこ（１）とですから，決定論的な計算機のアルゴリズムからこれをこ こ で パ ラ メ ー タ #&，$，"は い ず れ も 非 負 の 整 数 で生成するというのはおかしな話にも聞こえます．

#&""，$""であり，$&%"は（１）式右辺を "で割った

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

2. Random Number Generation

S. Satake and R. Kanno

余りを取る，という意味です．したがって疑似乱数 .
%はの 実 数 に 変 換 し た 場 合，(次 元 空 間 に"より小さい非負の整数になります．

，'#!
!
!
!
"
#
) と点を並べていく
+
+
+
*!
%
('""!
('"#!
(%
'""&&

これがどのような疑似乱数となるかはパラメータと，結晶のように規則的な格子状に点が並んでしまうこと$!
"&の選び方に強く依存します．疑似乱数のアルゴリ#&!
%を指します．例えば混合合同法の場合，２次元平面に並べズムを考える上で重要となる性質はその「統計的性質」，

るとこの結晶構造が現れます．これは係数%
#!!
$
!
"&をどの

「周期性」，「再現性」の３つです．統計的性質とは生成さように選んでも現れます．また，"は周期長が最大になるれた性質がいかに乱数っぽく振る舞うかのことで，詳しいように，疑似乱数のバイト長と同じ（４バイト整数なら

説明は第４章で行いますが，ここでは私たちが一般的な確"##$#とするなど）に取ることが一般的ですが，"が偶数

率的事象に対して期待するような振る舞いを，疑似乱数がの場合，一番下位のビットが必ず 0，1 が交互に現れる，つどれだけ上手く模擬できているかのことだと考えて下さまり偶数，奇数が交互に現れるといった性質もあります．

い．周期性についてはこの線形合同法に限らず，有限の有乱数を用いたシミュレーションの中には，このような規則

効桁数で演算されるコンピュータの漸化式による疑似乱数

性が計算結果に影響を与えてしまうものもあるため，疑似発生法には必ず固有の周期長がある，というのが重要な点

乱 数 と し て 望 ま し い 性 質 と は 言 え ま せ ん．周 期 長 がです．つまり，（１）式のような漸化式で数列を生成してい#$"$#!億というのも，１００万個の粒子を使ったモンテカル

くと，どこかでまた同じ数列の繰り返しになるということ

ロシミュレーションを２０００ステップ進めただけで１周期使です．また，決定論的なアルゴリズムなので，同じ初期条

い切ってしまうということで，現在の大規模計算において

件'
.
.
.
.
*!
!&!
!&
""!
!"!
%
%
%
%(を 与 え れ ば 同 じ 乱 数 列は十分に長い周期とは言えなくなってしまいました．

.
.
* が毎回得られる，つまり疑似乱数には再現性
""!
"#!
%
%

そこで，線形合同法とは異なる疑似乱数がい ろいろがあります．乱数なのに再現性があるというのもおかしくと 考 案 さ れ て き ま し た．こ こ で は そ の 代 表 例 と し て聞こえますが，同じ計算を同じ条件で走らせた場合に同じ

Tausworthe（トーズワース）の方法と，その発展形である

結果が得られるかどうかは，シミュレーションコードの開

Mersenne-Twister（メルセンヌ・ツイスター）法について

発時のバグ取りや検証のためには必要不可欠な性質です．

簡単に紹介します．

線形合同法のうち，特に &#!で，.
+-,"&
!.
""##
%
%%
!Tausworthe 法と い う 最 も 簡 単 な 形 を と っ た も の は 乗 積 合 同 法，.
+-,"&の形のものは混合合同法とも呼ば
%
!.
""##
%
%"$
れ，一つ前の値 .
%のみ記憶していればよいので，メモリや造を示さない疑似乱数として知られており，M 系列と呼ば計算速度に制限の多かった古い時代からよく使われた疑似れる疑似乱数の仲間の一つです．!%を -ビットの２進数を乱数です．なお，現代では混合合同法のことを線形合同法

表す横ベクトルとするとき，以下の漸化式

Tausworthe 法は上に述べた多次元分布における結晶構と呼ぶことが多いようです．

!%") #!%"*! !%C 言語や JAVA，Visual Basic，Excel や UNIX OS の標準（ただし )$*$!）（２）

的なライブラリーに入っている乱数ルーチンの中身は，現によって生成されます．ここで，! は各ビットごとに排他在でも混合合同法であることが多いようです．こうした標的論理和（XOR）を取ることを意味します．M 系列乱数に準化された乱数ルーチンでは，疑似乱数列の初期化のためついて詳しく説明するには，'!
!
"(である２つの要素に対に seed という値を最初にユーザが指定することになってするガロア体 GF（２）における議論が必要になりますの

いますが，これは乗積合同法や混合合同法における最初ので，詳しくは［１］を参照ください．ここではごく簡単な説明

.
#!!
$
!
"&は，発生する!を与えているに過ぎません．係数 %にとどめておきます．
（２）式に対して GF（２）上の多項式

乱数列の周期長がなるべく長く，乱数としての性質がよく.)".*""を特性多項式と呼びます．もし，.,!"という式なるように選ばれていますが，原理的にその周期長は "が ,##)!"の場合この多項式で割り切れ，,"#)!"に対を超えることはありません．例えば，C 言語の stdlib に実装

しては割り切れない場合，この特性多項式を特に )次の原されている RAND
（ ）関数は４バイト整数の疑似乱数ルー

始多項式と呼びます．このとき，（２）式から生成される数

チ ン で す が，そ の パ ラ メ ー タ は #!#""!$&"&#%&，

列の各ビットに着目すると，それぞれの '!
!
"(の並びのパ$#

ターンの周期が #)!"となることが知られています．原理$"$#"#$%&，"## であり，周期長は # です．
さて，古くから使われてきた線形合同法ですが，その乱的に，（２）式の漸化式が作る数列の周期は，#)!"より長数としての性質にはいくつかの問題があることが知られてくなることはないので，原始多項式から生成される数列はいます．「乱数らしさ」とは何かについては，また第４章

最大周期長を持つことになります．

で詳しく述べますが，線形合同法に関する問題点として

線形合同法と同様に Tausworthe 法でも，原始多項式とは，２進数で表現した乱数の下位ビットに，乱数列全体となるような %)!
*&の組み合わせを上手く選ぶ必要がありましての周期よりも短い周期性が現れること，「結晶構造」

す．代 表 的 な 例 と し て，%
，%
)
!
*&
!
!
'!(
#($&
"
#(*
%")&
#%

と呼ばれる分布の偏りが見られること，そして最近使われ等が知られています．さらに Tausworthe 法では )個のるようになった他の新しい疑似乱数ルーチンに比べて周期-ビットベクトル !%を初期条件として与える必要があり，長が短いことが挙げられます．
「結晶構造」とは，例えば漸化式を回すために )個分の過去の情報を保存しておく必! により%!".
!
!
"&区間%"! の整数疑似乱数列を +%#.
%#要があります．一つ留意しなければならないのは，原始多２９３Journal of Plasma and Fusion Research Vol.96, No.6 June 2020作を行ったベクトル 'を出力として取ります．

項式によって生成される数列が最大周期長を与えるのは.$"ビットの場合であり，.'#ビットになった場合に

，
&
%&-)
'(
"$%

は，%$ベクトルの初期条件を上手く設定しないと，. ビッ

，
&
&
)& #)
#$&
"'(
"%+

トベクトル %$全体のビットパターンの周期長が#(!"より

，
&
&
)& $)
$$&
#'(
#%,

短くなってしまう場合がある点です．幸い適切な初期ベク

（４）

，
'$&
&
)
$'(
$&%

トルの与え方は Tausworthe 法のサンプルプログラムに付随していますので，ユーザは自分で適当に初期ベクトルを

ここで -，+，,，%はある整数定数，(#!
$)は . ビットの定

定義せずにそれを利用しましょう．サンプルプログラムは数ベクトル，そして &-，%%はそれぞれ -ビット右シフ［３‐５］などで見つけることができます．

ト，%ビット左シフトを表します．
（３）式と（４）式は一見

Tausworthe 法の疑似乱数の性質については，多次元分

複雑な演算を行っているように見えますが，XOR，AND，

布における一様性が良好であることが特徴です．また，線ビットシフトなど簡素な２進数演算のみで構成されてお"%形合同法の周期長 "はプログラムが扱える最大の整数（#

り，コード化すると思いの外短く記述できます．

や #$#）を 超 え ら れ ま せ ん が，Tausworthe 法 の 場 合，MT 法の特徴の１つはその周期長の長さにあります．適(!
))の上手い組み合わせさえ見つけられれば，(&.と取(切にパラメータ &，'，*，"を選ぶことにより，生成されれるため格段に長い周期長 #(!"を実現できます．XOR

る数列 %$の周期長が #(!"$#'.!*!"になることが数学

演算もCPUで高速実行されますし，使用するメモリ量も線的に証明されています［７］．先に挙げた Tausworthe 法と形合同法より少し増える程度です．しかしながら，多次元

同様，その肝は漸化式（３）を決めるパラメータの選び方に分布の均等性以外の乱数としての性質については必ずしもあります．その詳細については数学の専門家ではない筆者

よいというわけではありませんでした［６］．また，（２）式の理解を超えるため，考案者の講義資料［３］を参照いただ

において各 .ビットが独立に計算されるため，ビット間でくことにします．ごく簡単に説明すると，数列の漸化式の情報のやり取りがないことが欠点として指摘されています．

形を表す特性多項式が /(!"の原始多項式になる，そのようなパラメータの組み合わせを見つけ出すことが一般的に!Mersenne-Twister 法は非常に難しいのですが，(が Mersenne 数と呼ばれる，

そこで，更なる改良として考案されたのが Mersenne-($#'!"の形を取る特別な素数の場合には，特性方程式

Twister 法［７（］以下 MT 法と略す）です．これは，日本のが原始多項式か否かの判定が容易になるという性質を使っ数学者，松本 眞氏によって考案された方法で，非常に長いているそうです．Mersenne-Twister の名前もこれが由来

周期性を持つ疑似乱数列が生成でき，乱数としての性質もになっています．また，（４）式のような操作は乱数の高次

良いことから，今では世界中の乱数を使うシミュレーショ

元分布の均等性をより高めるための工夫です．
「高次元分

ンで幅広く利用されるようになりました．

布の均等性」の説明は第４章で改めてしますが，例えば下に示す MT 法のサンプルでは周期長が #"''$&!"で，６２３次

MT 法の漸化式は形式的に以下のように記述されます．
.!* *%$"' $%$"& '(
%$
%$"")
!
,

元空間に均等分布する疑似乱数を生成します．

（３）

MT 法 に つ い て は，様 々 な 人 が Fortran，C，Java，

ここで '"&""です．また，%$は（２）式と同様に . ビッPython など様々な言語で実装しており，MT 法の考案者のトの２進数を表す横ベクトルであり，!は *!
!
"+を要素と松本氏のwebページにはそれらへのリンクがまとめてありする .#. の正方行列で，次のような形を取ります．
"
!
! " !
&
%
&
%! ! "&%&%&%! ! !
&
!$%
&
%
&
%
&
%
&
%
"
$
#
#! #" ## - #.!"

ます［８］．なお，初期のバージョンの MT 法のソースコー
ドには，乱数の初期化の seed の与え方があまり良くないものが使われており，紹介されているコードの中にも旧い初期化の仕方のままになっているものもあるので注意が必要
です．C言語については，MT法の作者自身が２００２年に公開した修正版（mt19937ar）［９］か，更にSIMD化して改良された SFMT 版（［９］中のリンク）をダウンロードして利用するのがよいでしょう．

ここで，"はある定数ベクトルで，演算 %!は以下のような

ここでは Fortran 版［１０］の利用法について紹介します．

演算に相当します．

なお，他の疑似乱数と同様，MT 法でも漸化式の形を決めるパラメータ("!
#!
$)を適切に選ぶ必要があ
&!
'!
*
!
-!
+
!
,
!
%
!

%!
!
%を左に１ビットシフト （%の最下位ビットが０）＝$"%を左に１ビットシフト ' "（%の最下位ビットが１）$#

りますが，これらは定数としてサンプルコード中に指定されているのでユーザが特に気にする必要はありません．
文献［１０］の mt19937ar.f は初期化と発生のルーチンのみ含

.!* *
また，%$
%$""は %$の上位 .!*ビットと %$""の下位,

まれています の で，これらを利用するためのメインプ*ビットをつなぎ合わせて１つの .ビットベクトルにする

ログラムを書きましょう．以下にサンプルプログラムという操作です．MT 法の漸化式は（３）式ですが，乱数とmtsample.f90を記載します（GitHubにも同じコードを公開して利用する際には，更に以下のような一連のかき混ぜ操してあります［１１］）．
２９４

Lecture Note

2. Random Number Generation

S. Satake and R. Kanno

input seed number (0 to continue)
program mtsample
implicit noneと聞かれるので，適当な自然数（４５６７など）を入れます．次

integer(kind=4), parameter :: N=624に，

integer(kind=4) :: seed,i,nrnd
input number to be generated
real(kind=8),allocatable :: rnd(:)

character(len=1) :: yesnoと聞かれるので，例えば１００と入れてみましょう．すると即!
!
"'区間疑似乱数が生成さ座に１００個の倍精度実数 の &real(kind=8) :: genrand_real1

れ，”mt_rand.txt”に格納されます．最後に，integer(kind=4) :: mti,initialized, mt(0:N-1)

common /mt_state1/ mti,initializedrecord the MT status? (y/n)
common /mt_state2/ mtと聞かれるので，y と入力すると乱数列を次回の計算に継!
print *,”input seed number (0 to continue)”

続するための情報が”mt_cont.dat”に記録されます．前回

read(5,*) seedの続きから疑似乱数を発生させたい場合は，seed に 0 を入if(seed>0) then

れて下さい．なお，mt_rand.txt と mt_cont.dat は実行する
度に上書きされるので注意して下さい．

call init_genrand(seed)

ここで，「疑似乱数を継続する」というのは，例えばあelsecall mt_initln

る seed からスタートして１万個を一気に生成したのと，同open(7,file=”mt_cont.dat”,action=”read”)

じ seed から始めて５千個の乱数を２回発生させたので同read(7,*) mti,initialized,mt

じものになるようにするということです．これは，実際にclose(7)疑似乱数を使って，１回のジョブに実行時間制限のあるス
パコンなどで，非常に長いシミュレーションを複数回繋いend ifで計算する場合に必要となります．サンプルの MT 法の場!
print *,”input number to be generated”

合は，項数６２４の漸化式になっているのでその状態を記録

read(5,*) nrnd

して，次に実行する時に読み直すことで継続できます．
０個ず
図１は seed=4567 で１００個一気に生成した場合と，５allocate(rnd(nrnd))

つを２回継続で生成した場合の乱数の分布をプロットした!

ものですが，両者が完全に一致していることが確認できま!
!
"#す．なお，mt19937ar.fにはgenrand̲real1の他にも，"，

do i=1,nrnd
rnd(i)=genrand_real1() ! [0,1]-real8!
!
"#など様々な区間の乱数を生成する function が定義さ&

end doれているので，用途に応じて rnd(i)=genrand_real1() の所!

を切り替えてください．

open(8,file=”mt_rand.txt”,action=”write”)

なお，サンプル中で適当に与えた seed の値の役割です

write(8,’(f18.15)’)rnd(:)

が，これは乱数の初期値を設定するものです．Tausworthe!
print *,”record the MT status? (y/n)”法や MT 法では漸化式をスタートさせるために，数百の漸read(5,*) yesno

化式の初期値 $!!%を指定する必要があります．MT 法の初if(yesno==”y”) then期化では１項だけの漸化式によるシンプルな疑似乱数ルー

open(7,file=”mt_cont.dat”,action=”write”)

チンを使ってその初期値を与えています（線形合同法より

write(7,*) mti,initialized,mt

は凝ったものですが）．seed の値を変えることは漸化式のclose(7)初期値のセットを変えることに相当し，より具体的に言えば周期長#"&&$%!"の数列のどこを開始点にするかを疑似乱

end if!

数で決めることになります．Tausworthe 法の場合，初期

end program mtsample

値を上手く設定しないと乱数列の周期長が期待通りの長さにならない場合があると述べましたが，MT 法ではビット

実 行 す る た め に は，［１０］か ら ダ ウ ン ロ ー ド し て き た

間で情報をやり取りするために，初期値の設定と実現されmt19937ar.f とこの mtsample.f90 を合わせてコンパイルしる周期長に依存性はありません［６］．

てください．例えば gfortran なら

ところで，このサンプルプログラムで生成する数を百万など大きめに取ると，それなりに計算時間が掛かることがgfortan mt19937ar.f mtsample.f90

体感できると思います．実はこのサンプルでは，１回につとすれば良いです．できあがった a.out を実行すると，

き１つの疑似乱数を返す function を必要な回数だけ call
する形になっています．実用上は，配列 rnd(nrnd)に一気に nrnd 個の疑似乱数を生成できた方がコードも書きやす
２９５

Journal of Plasma and Fusion Research Vol.96, No.6 June 2020

射線や，ダイオードの出力雑音，量子ゆらぎなどの物理現象を利用しようというアイデアは，古くからあり，日本では，遅くとも１
９５７年には物理乱数発生器を計算機に接続し，モンテカルロ計算に物理乱数が利用されていたようで
す［１２］．以下では，定電圧ダイオード（ツェナーダイオー
ド）の出力信号のゆらぎ，およびレーザー光の偏光の量子
ゆらぎを利用した乱数発生の原理について説明します．
!ダイオードの出力雑音による乱数発生
図１

MT 法で生成した乱数のサンプル．１０
０個を１度に生成した
場合と，５０個ずつ２回継続して生成した場合をプロットしてあります．

まず，定電圧ダイオードを利用した乱数発生器についてですが，出力信号に含まれるわずかな（数 μV 程度の）ゆらぎ（ノイズ波）をノイズ源として利用します［１３］．このゆいですし，効率もよさそうですよね？また，現在では大規らぎを ADC（Analog to Digital Converter）の変換範囲まで

模シミュレーションは MPI などの並列環境で実行するの増幅し，ある一定間隔の離散時刻でノイズ波のデジタル値が当たり前になっていますが，疑似乱数の並列化はどのよを取得します．各時刻に得られたデジタル値を００００〜１１１１

うにすればよいのでしょうか？こうした効率化，並列化にまでの１６個の２進数のうちの１つに対応させることで，各関するテクニックについては次章で解説したいと思いま

時刻における４ビットの乱数（０または１の値で表現されす．

る二値乱数）が生成されます．２回のデジタル値の取得で８ビット（つまり１バイト）の乱数となり，これを繰り返

２．
２ 物理乱数すことで，必要なバイト数の二値乱数を作ることができま

本節においては，乱数の生成にランダムな物理現象を利す．Fortran や C などのプログラミング言語における単精!
!
"!に含まれる一様乱数を生成す
度浮動小数点型の区間 "用した「物理乱数」について紹介します．乱数の生成に放!"""""""""""""""""""

乱 数 茶 話!"""""""""""""""""""!"""""""""""""""""""""""""
少し余談になりますが，疑似乱数は数値シミュレー

数調整」と称して一見無意味な行動をしてからイベント

ションの世界で利用されるだけではありません．一般のを起こすと（宝箱を開けるなど）決まった結果が起こる，

人々が疑似乱数と係わるのは，例えば，コンピュータというようなことをたまに聞きますが，これらは乱数のゲームで遊んでいる時ではないでしょうか．敵キャラク

代わりに使われている上に述べたようなデータの隠れた

タの出現位置や行動などを，プレイヤーが飽きないよう

規則性を見抜いて，都合のよい結果になるように「調整」にランダムにするなどの用途に疑似乱数が使われていましているわけです．

す．ところで筆者（S.S.）が子どもの頃のコンピュータと今回の記事の執筆のために調べていて興味深かったの

言えば，８ビットの CPU に 32 kB のメモリという，非常が，８ビット PC の時代から線形合同法や XOR を使ったに貧弱な計算資源の中でゲームを動かしていました．当

M 系列の乱数も利用されていたことを知ったこ とで然，乱数発生にメモリを割いたり複雑なアルゴリズムをす．ただし８ビットなので周期が非常に短いですが．線採用したりする余裕はありません．当時のパソコンのメ

形合同法を使ったゲームでも，偶奇が交互に出る性質をモリ（RAM）は放っておくと電気的に記録されたビット

理解せずにくじのイベントの当たり判定に使ってしまっ

情報が減衰してしまうので，時々データを書き直す必要たために，プレイヤーに癖を見抜かれてしまったというがありました．このタイミングをコントロールするため事例もあったようです．筆者のいとこは以前パチンコ機に CPU 内部にリフレッシュカウンタというレジスタが器のプログラマーをしていましたが，少ないプログラム

あり，タイマのようにCPUのクロックごとに値が増えて行数で統計性のよい疑似乱数を発生させる方法はないか

いきました．当時のCPUはせいぜい数kHzで動作していと相談されたことがありました．パチンコでは不正防止

ましたが，それでも人間から見れば高速にレジスタの値のために，当たりの確率が規定された通りになっているが変わるので，このレジスタの値が疑似乱数としてよくことが認可を受けるために必要な一方，プログラムに使使われていたようです．他にも，CPU には，命令を実行

えるメモリ量にも制約があるため疑似乱数ルーチンの選するたびにその結果に応じてレジスタの特定のビットが択で悩んでいたようです．

０になったり１になったりする，フラグレジスタという

このようにゲームやギャンブルの世界では疑似乱数はものがあります．これは本来，直前に行った計算がオー広く使われていますが，人工知能が発達すると，やがてバーフローしたとか，論理演算の真偽値を確認するためは疑似乱数に潜む隠れた規則性をあっさり見抜いて，囲のものですが，このフラグレジスタの値を乱数として使碁や将棋だけでなくコンピュータゲームでも人間のプレ

うこともありました．古いゲームの攻略法などで，「乱イヤーを凌駕してしまうかもしれませんね．

!"""""""""""""""""""""""""
２９６

Lecture Note

2. Random Number Generation

S. Satake and R. Kanno

る場合は，例えば，４バイトの乱数を用いて得た正整数をして，一様乱数表を利用する方法をここでは紹介しました

その最大値に１を足した数で割ることで得られます．以上

が，別の方法もあります．詳しくは，文献［１３］や関連するは，４ビ ッ ト の 乱 数 を 基準 に 説 明 し ま し た が，一 般 に

特許情報などを参照ください．

! ビットの乱数としても同様です．

このように乱数発生原理の概略がわかると，次に乱数発

ところで，上記の乱数発生法を実際の装置で実現化する
際には，得られる二値乱数の一様性が気になるのではない

生の速度が気になるのではないかと思います．この発生速度については，近年，640 MB/s（#5.1 Gbps）程度になってでしょうか．そこで，乱数の一様性を向上させる方法についるようです．複数台の乱数発生器を用意できるのであれいて紹介します．乱数生成のスピードを落とさずに一様性

ば，大規模並列計算機におけるシミュレーションでの利用を向上させる手法として，排他的論理和による補正手法がに（必要とする乱数の数や呼び出す頻度にも依りますが）考え出されました［１３］．文 献［１３］に 従 い，こ の 手 法 をほぼ問題の無いレベルの発生速度であると言えます．

２ビットの乱数を例として説明します．ADC が２ビット

定電圧ダイオードを用いた物理乱数発生器による乱数をの変換値として０〜３を出力するとし，変換値 "の発生確利用したい場合は，例えば，統計数理研究所の乱数ライブ

率を #" とします．ここで，""!
!
!
!
"
#
$で，"$
" ""
""!#

ラリー［１４］における乱数取得サービスから乱数をダウンです．あらかじめ与えた２ビットの一様乱数表から等しい

ロードする方法があります．この記事を書いている時点で確率"で補正値０〜３を生成し，ADCの変換値との排他的は，東京エレクトロンデバイス製の物理乱数発生器による論理和を取ります．表１に ADC の変換値と補正値との排他的論理和を行ったときの補正後の値（上段）とその生成乱数を取得できるようです．取得できる乱数のデータ型!
!
"$の単は，区間［−231，231］の３２ビット符号付整数，区間%確率（下段）を示しました．ここで，排他的論理和は，２．
１

精度浮動小数点数および倍精度浮動小数点数です．

節でも使いましたが，以下のような演算（! と表記）であ!レーザー光の偏光の量子ゆらぎによる乱数発生

り，表１では，２ビットの値００〜１１について，ビットごとに演算して補正後の値を得ています．

次に，光学的な装置を使った「量子乱数」と呼ばれる物
理乱数発生器の原理について紹介します．量子力学的な不! ! ! " !

確定性を持つ物理過程の観測値を乱数のソースとして利用! ! " " "しようというアイデアは１
９５０年代には既にありましたが" ! ! " "［１４］，初めに考えられたのは放射性物質から出てくる放射" ! " " !

線をガイガー計測器で計測した時の，２つのパルスの発生

例えば，変換値１と補正値３で考えると，!"! """ "!

間隔が Poisson 分布の確率分布をとるという性質を利用すですので，補正後の値は，２となります．表１を見ると，補るものでした．しかし放射性物質を線源として内蔵するこ

正後の値が２となるケースは全部で４つあり，その生成確との問題や，ガイガー計測器の計測速度に原理的な限界が率の和は，#!"!#""!##"!#$"""となるので，補正後のあったこと（数 kbps 程度），小型のレーザー発振器が実用

二値乱数１０は確率 "で生成されることがわかります．他の化されたことなどから，現在では光学デバイスによる量子補正後の値についても，すべて同じ確率 "で生成され，こ

乱数発生器が主流になってきています．文献［１５］の Tableのことから，生成された２ビットの乱数の一様性が向上しII にいくつかの異なる手法による光学的な量子乱数発生器ていることが分かります．上記の説明は，簡単のためにの一覧が載っていますが，乱数の生成速度は数 Mbps から

２ビットの場合で説明しましたが，２を ! に一般化して，

数 Gbps に達しています．

! ビットの乱数とした場合も同様な補正を行えば，生成すレーザー光を使った量子乱数の例としてここでは，［１６］る乱数の一様性が向上します．また，補正値の生成方法とのレーザーパルスの位相差を利用したものを紹介します．
図２のようなレーザー光の回路を組み，発振したレーザー

表１

パルスを一度２つの光ファイバーのケーブルに分岐させ，

ADC の変換値と補正値との排他的論理和により与えられる補正後の値（上段）とその発生確率（下段）．例えば，変換値１と補正値３の補正後の値は，２行‐４列目の上段に表示され，２です．この表は，文献［１３］の表１を参考に作成
しました．

０
０
#!"

補正値
１
２
１
２
#!"
#!"

３
３
#!"

１

１
#""

０
#""

３
#""

２
#""

２

２
##"

３
##"

０
##"

１
##"

３

３
#$"

２
#$"

１
#$"

０
#$"

０

変換値

片方の経路を長く取ることで，一方をちょうど１パルス分ずらします．その後，２つのパルスを再び合流させてから

図２

２９７

レーザー光を使った量子乱数発生器の概念図（［１６］の Fig.2を引用）．Adapted with permission from [16] ! The Optical Society.

Journal of Plasma and Fusion Research Vol.96, No.6 June 2020

その光の強度を計測します．レーザー光のパルスは量子力

NIST テストと呼ばれる一連の乱数検定テストもパスして学的揺らぎにより，一つ一つ異なる位相 #"を持った偏光にいます．このケースでの量子乱数の発生速度は理論的にはなっています．したがって２つのパルスを合成した後の40 Gbps を超えるとされていますが，実際には，Hash 関数

レーザー光の強度は，干渉によって ##!!#!!!$の依存$%#の演算処理による時間を含めると実用上のスピードは性を持ちます．#!と #!!!が量子力学的な揺らぎから決まる

数 Mbps〜Gbps 程度のようです［１９］．

ため，#!と #!!!は互いに独立で一様分布する確率変数であ今回紹介したレーザー光の偏光を利用した手法の他にると期待できるので，そこから観測光の強度の確率分布がも，レーザー光を２つにスプリットしたものを再度合わせ計算できます．仮に観測光の強度を #"%!
!&に規格化し!!

て計測することで，観測に混じる量子力学的な真空のゆらたとすると，その確率分布は !
"!
' !#"に比例すると期待さぎを増幅して観測するという原理に基づいた光学的量子乱れます．しかし，実際の観測値には観測のノイズやパルス数発生器も開発されています［２０］．この手法も信号をデジ幅に関係する補正項が入るので，単純に理論通りに分布す

タルサンプリングした生データの情報の一部を捨てる代わるわけではなく，図３のようなひずんだ凹型の分布になりりに統計的性質を向上させる手法を使って，計測データにます．この分布は，２つのパルスに何らかの相関関係があ含まれる量子的ゆらぎ以外の計測ノイズの影響を落とす後るため生じているわけではなく，上に述べたパルス間の位処理を使います．この原理に基づく乱数発生器がオースト

相差がランダムであることに起因しています．しかし，こラリア国立大学の研究グループによって公開運用されておの測定値には計測誤差やレーザー光源の揺らぎなどの誤差

り［２１］，乱数発生速度は 5.7 Gbps に達するそうです．インによる期待値からのずれも含まれており，これらの誤差はターネットを通じて生成された乱数を利用でき，様々なプ

必ずしもランダムであるとは言えません．また，非一様なログラム言語からこの乱数サーバのデータを取得するルー

強度分布のままでは乱数として使うには適さないので，デチンが公開されているので，関心のある方は［２１］のリンク

ジタル化した計測データをWhirlpool-Hash関数［１７］というから情報をたどってみてください．

ものを使って変換したものを乱数として使います．これところで比較参考のために MT 法の発生速度を見る

は，入力されたデータをある決まった規則に従って変換と，２．
１節のサンプルプログラムを Xeon のワークステー

し，入力データか ら は 予 測 で き な い 固 定 長 の デ ー タ列ション（クロック周波数 2.5 GHz）で実行した場合，３２ビッ（Hash 値）を返すもので，暗号化の技術の一つとして使わ

ト乱数を１千万個生成するのに約 0.2 秒，発生速度に直す

れているものです．この際，ある確率分布を持つ確率変数と約 1.6 Gbps でした．筆者（S.S）は数年前に［１６］の量子乱の最低エントロピー状態という理論を基に，デジタルサン

数発生器の試作機のテストに参加させてもらいましたが，プリングした元データのビット数よ り 低 い ビ ッ ト数の

当時はまだ乱数発生器からコンピュータのメインメモリにHash 値を乱数として採用することで，元のデータの情報

データを転送するところが律速となり，理論性能ほどの速量を一部落とす代わりに一様性の高い乱数列を得るという度は出ていませんでした．しかしそれでも当時のワークス手法［１８］を取っています．具体的には，１４ビットでサンプ

テーション用CPUでMT法を走らせた場合の１０分の１程度

8

リングした 1.2×10 個の観測データを Whirlpool-Hash 関数の速度が出ていました．物理乱数発生器は，まだあまり世で 1.25×108 個の７ビットのバイナリデータに変換していの中に出回っていませんが，［１４］や［２１］の乱数サーバの実

ます．こうして得られた数列は非常によい一様性を示し，績値を見ると，物理乱数の発生速度は疑似乱数に比較しう
る時代になっていると言え，今後の展開が期待されるとこ
ろです．

謝

辞

第２．
２節における定電圧ダイオードを利用した物理乱数の発生原理の詳細についてご教示いただいた，田村義保氏（情報・システム研究機構 統計数理研究所 特任教授・名誉 教 授）に 感 謝 い た し ま す．ま た，第２．
２節 に 対 し て，
河村学思氏（自然科学研究機構 核融合科学研究所 助教）より，貴重なコメントをいただきました．ここに感謝申し上
げます．

参考文献
図３

［１］津田孝夫：モンテカルロ法とシミュレーション（三訂版）（培風館，１９９５）．
［２］伏見正則，逆瀬川浩孝（監訳）：モンテカルロ法ハンド
ブック（朝倉書店，２
０１
４）．
［３］松本 眞：有限体の擬似乱数への応用，
http : / / www. math. sci. hiroshima-u. ac. jp / 〜 m-mat /短，長２つの経路を通ったそれぞれのパルス強度（〇と□の印，左の縦軸）と，２つを合わせて観測した場合の強度（＋印，右の縦軸）の分布．２つのレーザーパルスを合成す
ると干渉によって広がりを持った分布になります（［１６］のFig.3（a）を引用）．Adapted with permission from [16] !
The Optical Society.

２９８Lecture Note

2. Random Number Generation

S. Satake and R. Kanno

たけ

しん

すけ

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

かん

の

りゅう

た

ろう

菅 野 龍 太 郎
自然科学研究機構 核融合科学研究所 ヘリ
カル研究部 核融合理論シミュレーション
研究系 准教授．乱数を用いた確率論的な計算手法であるモンテカルロ法全般に興味があります．最近の研究では，プラズマの衝突輸送現象に対するモンテカルロ法に基づいたドリフト運動論シミュレー
ションを行っています．

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

さ!!!!!!!!!!!!!!

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

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

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

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

［１２］統計数理研究所 統計科学技術センター 計算機の歴史
https://www.ism.ac.jp/computer̲system/jpn/outline/
hist.html［１３］田村義保 他：日本統計学会誌 35, 201 (2006).
［１４］統計数理研究所 乱数ライブラリー
http://random.ism.ac.jp/［１５］Miguel Herrero-Collantes et al., Rev. Mod. Phys. 89,
015004 (2017).
［１６］C. Abellán et al., Optics Express 22, 1645 (2014).
［１７］P.S.L.M. Barreto and V. Rijmen, The Whirlpool hashing
function, https://www.researchgate.net/publication/
228610491［１８］N. Nisan and A. Ta-Shma, J. Comput. System Sci. 58, 148(1999).
［１９］Carlos Abellán et al., Optica 3, 989 (2016).
［２０］T. Symul et al., Appl. Phys. Lett. 98, 231103 (2011).
［２１］The ANU Quantum Random Numbers Server,
https://qrng.anu.edu.au/!!!!!!!!!!!!!!!!

TOKYOHOMEPAGE/TEACH/1011.pdf［４］統計数理研究所 乱数ライブラリー 乱数について
http://random.ism.ac.jp/info01.html［５］宮村 修，牧野 純：ベクトル化乱数発生プログラム，大阪大学大型計算機センターニュース 75, 39 (1989).
http://hdl.handle.net/11094/65855［６］松本 眞，栗田良春：Twisted GFSR：新しい乱数発生法，京都大学数理解析研究所講究録 85, 86 (1993).
http://www.kurims.kyoto-u.ac.jp/〜kyodo/kokyuroku/
contents/pdf/0850-08.pdf［７］M. Matsumoto and T. Nishimura, ACM T MODEL COMPUT S 8, 3 (1998).
［８］http://www.math.sci.hiroshima-u.ac.jp/ 〜 m-mat/MT/
VERSIONS/versions.html［９］http://www.math.sci.hiroshima-u.ac.jp/ 〜 m-mat/MT/
MT2002/mt19937ar.html［１０］http://www.math.sci.hiroshima-u.ac.jp/ 〜 m-mat/MT/
VERSIONS/FORTRAN/mt19937ar.f［１１］https://github.com/satakeshinsuke/random̲number!!!!!!!!!!!!!!!!!!!!!!!!!!

２９９

J. Plasma Fusion Res. Vol.96, No.7 (2020)３
６７‐３７１

講座 シミュレーションのための乱数入門

３．乱数発生の高速化と並列化技法
3. Techniques of Acceleration and Parallelizing the Random Number Generation
佐竹真介
SATAKE Shinsuke
自然科学研究機構 核融合科学研究所（原稿受付：２
０
２
０年３月１
７日）

現代のスーパーコンピュータを用いた大規模シミュレーションは，多数の計算ノードを用いた並列計算が標準的になっています．本章では，MPI 分散並列プログラムにおける疑似乱数の並列発生方法について，数値計算
ライブラリ KMATH̲RANDOM を例に解説します．また，乱数発生の高速化チューニングの実例や，物理乱数発生器を使う場合の高速化，並列化についてのアイデアについても紹介します．
Keywords:
pseudo random number, physical random number, parallel computing

・はじめに依存しますが，基本的に分散並列の各プロセスは異なる乱

本章では，大規模シミュレーションでの乱数の利用を想

数列を利用しなければならないということです．例えば，定し，乱数生成を効率的に行うためのプログラミング上の技

分子動力学（MD）シミュレーションでタンパク質の折り畳法について解説します．現代の大規模シミュレーションではみ構造を調べる計算を，１０プロセスの分散並列で行うとし

特に，MPI や OpenMP といった分散並列，スレッド並列をましょう．異なる初期配置を持ったタンパク質分子を，乱駆使したプログラムを利用することが一般的になっていま

数を利用して多数用意する際に，それぞれのプロセスで同

すから，そのような並列化プログラムの中で乱数を利用すじ疑似乱数発生法，例えば MT 法で，第２章で説明した

る上での留意点や工夫について紹介したいと思います．な

seed の値も同じに取ってしまったら，１０個のプロセスで全

お，本 章 で は 疑 似 乱 数 発 生 法 と し て 第２章 で 紹 介 し た

く同じ計算を開始してしまうことになり，並列化の意味がMersenne-Twister（メルセンヌ‐ツイスター）法（MT 法）を 対ありません．そんなことはわかりきったことであり，seed

象として説明しますが，物理乱数器の利用を想定した場合の値をプロセス毎に変えればよい，と考える人もいるかもの並列化や高速化に関するアイデアも紹介したいと思います．

知れませんが，それはお勧めできません．疑似乱数の seed

まず乱数の並列化の説明に先立って，MPI による分散並の値は，利用する乱数列の開始点を一周期のどこに取るか

列と，OpenMP（あるいはコンパイラの自動並列化機能にを決めるものだと２．
１節で説明しましたが，例えば seed のよる）スレッド並列の違いを明確にしておきましょう．分

値を１０個のプロセスででたらめに（人手で決める，時計を散並列は，計算を実行する CPU がそれぞれ物理的，あるい使う，他の疑似乱数を使うなど）与えても，それらの部分は論理的に分割された異なるメモリ領域に配列や変数を記

列の開始点が互いに十分離れているという保証は全くあり

憶し，シミュレーションに必要な計算を各CPUが分担してません．seed の選び方がたまたま悪くて２つのプロセスの行い，必要に応じて他のCPUと明示的な通信命令によって乱数列の開始点が数個しかずれてなかった，ということもデータをやり取りする並列計算です．一方，スレッド並列

考えづらいですが，起こりえない話ではありません．各プ

は，CPU の中にある複数のコアが，同一のメモリ上のデー

ロセスが参照する乱数の部分列が確実に重複しないようにタを使った演算を分担するやり方で，基本的には Fortran

選ぶことは，プロセス間で擬似乱数によるおかしな相関関の do loop や C 言語の for loop などの繰り返し計算を分割し係が生じる懸念を排除するために必要なことです．それでてそれぞれのコアが並列計算する方法です．この両者を組は，具体的にどのように乱数発生を並列していけばよいのみ合わせたハイブリッド並列も大規模シミュレーションでか見ていきましょう．

は一般的に行われています．以下，分散並列における並列・Mersenne-Twister 法による乱数生成の並列化

化の区切りを「プロセス」，スレッド並列におけるそれを分散並列計算における MT 法の乱数発生の並列化につい「スレッド」と呼んで区別します．
さて，並列シミュレーションにおいて乱数を利用する上ては，以前にもプラズマ・核融合学会誌の講座で紹介させで留意しなければならないことは，プログラムの詳細にもていただいたことがあります［１］．その時に紹介したのは，

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

３６７!2020 The Japan Society of Plasma
Science and Nuclear Fusion Research

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

Dynamic Creator と呼ばれる，MT 法の考案者の松本眞氏

コンパイラとコンパイルオプションの選択，NTL をインス

ら自身が開発した技法［２］で，ソースコードも公開されてトールしたディレクトリの指定です．random/arch/以下

い ま す［３］．Dynamic Creator で は，第２章 で 説 明 し た，にいくつかのコンピュータにおけるサンプルがあるので，

MT 法の漸化式の形を決めるパラメータのうち，ベクトル

それを参考にご自身の環境に合わせて設定します．なお，!!
"!
#$を変えることで，特性多項式が互いに素で，最大#コンパイルするには MPI 環境も必須です．また，MT 法の!

周期長がある Mersenne 素数 !に対し " !!となり，かつ

周期長 "!!!はここでコンパイルオプション “-DDSFMT_

高次元に均等分布する疑似乱数の漸化式の型を多数用意す

MEXP=p”で指定します．make コマンドを実行し，コンパイる，という方法を取ります．数学的に，特性多項式が互いルが完了すると，random ディレクトリ下の c/，c++/，f90/に素な漸化式同士は相関を持たないとされているため，

以下に libkm_random.a というそれぞれの言語用のライブ

Dynamic Creator は言わば互いに独立な乱数表を複数用意

ラリができるので，自分のプログラムで利用するにはこれすることに相当します．ただし，特性多項式が互いに素なをリンクしてコンパイルすることになります．

パラメータ #!!
"!
#$の探索はトライ＆エラー的になされる

次に，ジャンプファイルという N 個飛ばしの疑似乱数列

ため，数千並列の並列計算のための準備には非常に長い時の開始点の状態を記録したファイルを生成します．まず，間がかかります（筆者が十数年前に１０２４個用意するのに，

make ptool を実行すると，random/ptool/に km_rand_gen

ワークステーション１台を１か月以上稼働し続ける必要が_jump ができるので，これを実行します．

ありました）．

（例）km_rand_gen_jump -seed S1 S2 -max_ranks R -rand

そこで，超多並列計算のための実用的な MT 法の並列化_range M

として今回紹介したいのは，Jump Method と呼ばれるもの

ここで，-max_ranks の値Rは想定される並列計算の最大プです［４］．MT 法の特徴である長大な周期長を活かし，１

ロセス数，-rand_range の値 M はジャンプの間隔 N=2M

つの疑似乱数列の中の，遠く離れた複数の点を並列計算での指数部分を指定します．-seed の値 S1，S2 は，いわゆる乱数を利用する開始点として明示的に指定する方法です．

通常の意味での seed の値を S1，S1＋1，…，S2 と変えたもこれは，先に述べたような seed の値を単に変化させるのとのを（S2−S1＋1）個生成するということです．これは N 個は異なり，開始点間の距離を #"!個ずつ飛ばしに指定す飛ばしで R 個用意する開始点全体の起点が異なるジャンプ

ることができます．N 個飛ばしした後の漸化式の内部状態

ファイルを複数作成することを意味します．通常，１つのを並列プロセス数 P だけ事前に用意するには，原理的にはシミュレーションに使うジャンプファイルはどれか１つをN 個乱数を発生させては内部状態を記録する（第２章のサ

選べばよいですが，異なるジャンプファイルは，乱数をンプルプログラム に お け る”mt_rand.txt”の デ ー タに相使ったシミュレーションを複数回行って統計平均を取りた当）という操作を P 回繰り返せば可能です．ですが，P 個い場合などに利用します．試しに S1＝1，S2＝2，R＝200，の乱数列の部分区間が重ならないように N を十分大きく

M＝100で実行すると，jump/ディレクトリに file_00001，（１００兆など）取り，かつ P が数千となると，実用的な方法

00002 の２つのジャンプファイルのセットができます．ことは言えません．しかし，Jump Method に使われる SIMD-

れで，MT 法の疑似乱数に対し，2100 飛ばしの開始点の情報

oriented Fast Mersenne Twister（SFMT）という新しいタが２００個×２系列作られました．この計算は例えば最近のイプの MT 法［５］では，［６］に示される原理に基づいて N

Xeon CPU 上では数秒で終わってしまいます．これは逐次個飛ばし後の MT 法の漸化式の内部状態に一気にジャンプ

計算で 2100×200 個×２回，通常の MT 法で乱数を発生させできる計算ルーチンが用意されているので，それを利用す

るのとは比較にならないほど速いです．

ることで乱数列の遠く離れた先の開始点から始めるための

これで前準備は完了です．次にこのジャンプファイ

内部状態を高速に取得できます．SFMT は名前から推察で

ルを使って実際に並列計算で SFMT を利用する方法をきるように，高速化チューニングされた乱数発生ルーチン

見ていきましょう．第２章と同様に，サンプルプログラム

です．従って Jump Method の利用に限らず，逐次的に MT

rand_test.f90 と km_rand_wrapper.f90 を GitHub に挙げ

法で疑似乱数を発生させる用途としても，従来のバージョておきます［８］．ifort コンパイラを利用する場合のコンパンよりも SFMT を用いた方が高速実行が期待できます．

イル用 makefile のサンプルが用意してあるので，ご自分

さて，実際の利用法についてですが，ここでは［５］のオリの 環 境 に 合 わ せ て 修 正 し て make し，実 行 プ ロ グ ラ ム

ジナルのソースではなく，C，C++，Fortran９０での実装がkmrand_time を作成してください．実行方法ですが，この用意されている，理化学研究所の大規模並列数値計算技術

サンプルプログラムは MPI とスレッド並列のハイブリッ

研究チームが整備・公開しているライブラリ"KMATH̲

ドコードになっているので，環境変数 OMP_NUM_THREADS にRANDOM"［７］をベースに紹介したいと思います．このラ

スレッド並列数を設定し，また KMATH_RAND_JUMP_FILE_PATH

イブラリの利用許諾についてはマニュアルの記載事項を参に先ほど作成した file_00001 等のジャンプファイルが格照してください．まず前準備として，NTL という数値計算

納されているディレクトリを絶対 path で設定します．そしライブラリを別途ダウンロード，インストールする必要がて，”mpirun -np N kmrand_time”のような MPI 実行コマンありますが，［８］に説明があるので詳細は省きます．次に，ドで計算を開始してください．ここで N は試したい MPIKMATH̲RANDOM を 展 開 し て，random/Makefile.ma-

並列数（N はジャンプファイル作成時の-max_ranks の値 R

chine というファイルを編集します．必要なのは，利用する以下）です．
３６８

Lecture Note

3. Techniques of Acceleration and Parallelizing the Random Number Generation

S. Satake

このテストプログラムでは，rand_test.f90 の冒頭に指る場合，ある程度の量の乱数列に対する変換をスレッド並定された，１プロセス当たり ntest 個の乱数を，ntime 回発

列で一気に行った方が効率が良さそうです．更に，シミュ生させてその平均時間を計測します．そして，そのテスト

レーションの中で不定長の乱数を発生させるためにを nseed 個の異なるジャンプファイルに対して繰り返しまKMATH_Random_Vector を高頻度に呼ぶより，ある程度の大す．計測結果は"log.txt"に書き出されます．rand_test.f９０

きさの，固定長 nbuff の乱数列を一気に生成し，必要な分の主要部分を抜粋すると以下のようになっています．

だけ切り崩して使っていき，足りなくなったらまた一気に-----------------------------------------

作る，というやり方にします．井戸水を使う人は，使うた

call MPI_INIT(IERR)

びに必要な量を井戸まで汲みに行くのではなく，一定量水

call MPI_COMM_SIZE(MPI_COMM_WORLD,nmpi,ierr)

瓶にため込んでおいてそこから必要な量だけ柄杓ですくっcall MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ierr)て使うのと同じように，ある一定量の乱数を貯めこんでか

call KMATH_Random_Init(handle1,MPI_COMM_WORLD, ierr)ら使う方法を，私は「水瓶方式」と呼んでいます．そのような一連の作業をしているのが grnd_dsfmt になります．
ここで，最も効率が良くなる nbuff を選ぶことで高速化

do k=1,nseed
call KMATH_Random_Seed(handle1, k , ierr)につながりますが，その値は計算機の環境によります．

do j=1,ntime

図１は，IFERC-CSC の ス ー パ ー コ ン ピ ュ ー タ JFRS1 上
で，8 MPI×5 スレッド並列で kmrand_time を様々な nbuff

call grnd_dsfmt(rnd1, ntest)と ntest に対して実行した際の，乱数１００万個当たりの平均

end do

発生時間を表しています．ntest が大きくても小さくても，

end donbuff=64000 の時に最も乱数発生時間が短くなりました．
call KMATH_Random_Finalize(handle1, ierr)

この例では，チューニングの効果はあまり大きくありませcall MPI_Finalize(IERR)

んが，他の計算機環境では，特に ntest が小さい時に nbuff--------------------------------の選び方で数十％の速度の差が出ることもありました．な

また，上のコードの流れは，以下のようになっています．

おこの水瓶方式ですが，grnd_dsfmt は MPI の各プロセスで１．MPI_INIT で MPI 並列計算を起動

独立に動くので，それぞれのプロセスで乱数の発生数やタ２．KMATH_Random_Init で KMATH_RANDOM の利用準備

イミングがずれていても問題ありません（並列化乱数発生

３．KMATH_Random_Seed で k 番目のジャンプファイルを読法とは当然そうあるべきものですが）
．また，
２．
２節で紹介

み込む

したような物理乱数を利用する際も，物理乱数のデバイス

４．grnd_dsfmt(rnd1, ntest)で配列 rnd1 に ntest 個の乱にプログラムからアクセスして乱数データを取得するオー数を発生させる

バーヘッドタイムを考慮すると，必要に応じて乱数を取得

５．KMATH_Random_Finalize で KMATH_RANDOM の利用終了

するよりも，水瓶方式を用いてある程度の量の乱数を一気

６．MPI_Finalize で MPI 並列計算の終了に取得する方が効率がよいと考えられます．

このプログラムでは nseed 個のジャンプファイルについ・物理乱数発生器を利用した並列化て，各プロセスで発生した乱数の最初の ntest 個を，”rand_s.txt”に全プロセス分をまとめて書き出します（全部をさて物理乱数の話が出たので，次に分散並列計算を行い出力すると膨大になるので 1/100 に縮めてあります）．プ
ロセス毎，seed 毎に異なる乱数列が生成されることが確認できると思います．

・KMATH̲RANDOM の高速化チューニング
さて，このサンプルプログラムで発生速度が計測できる
ので，次に高速化について考えましょう．乱数の発生は
grnd_dsfmt というサブルーチンで行っていますが，これはkm_rand_wrapper.f90 の 中 に 定 義 さ れ て い ま す．実 は，KMATH_RANDOM が 用 意 す る 乱 数 発 生 ル ー チ ン KMATH_Random_Vector は，高速化のために区間 !
!
"#"の８バイト実数を，３８６個以上の偶数個発生するという少し変わった作りになっています．なので，１回の必要な発生数が少ない場
合や，奇数の場合の対処法が必要です．また，通常シミュ
レーションで使われる乱数は !
!
!
""区間に発生させた一様乱数を利用することが前提になっていることが多いですか

図１

ら，KMATH_Random_Vector で発生させた乱数を!
!
!
""区間に
変換する操作が必要です．スレッド並列が有効になってい３６９

水瓶方式でバッファの大きさ nbuff を変えた場合の，乱数の発生時間の変化（百万個当たり）．nbuff ＝ 64000 で極小値となります．

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

たいが物理乱数発生器が１つしかない場合，乱数の並列化をどうするか考えてみましょう．図２のように，計算の１ステップが Calc. A, B, C, …といったルーチンからなり，Calc.B において 乱 数 が 必 要 で あ る 場 合 を 想 定しま
す．この時，例えば量子乱数発生器（QRNG）が rank0 の
MPI プロセスが走るノードにのみあるとします．考えられる並列化の方法としては，rank0 は物理乱数の取得と，それを他のプロセスに分配することに専念し，他のプロセス
図２

が計算を担当するというやり方です．この際，「１ステッ

MPI 並列プログラムで量子乱数発生器（QRNG）が１つしかない場合の「源泉かけ流し方式」におけるプログラムの流れの模式図．rank0 のプロ セ ス は 乱 数 生 成 に 専 念 し，Calc.A が終わった他の MPI プロセスは乱数を受信し，乱数を利用する Calc.B に備えます．

プ当たり必要とされる乱数の最大個数（プロセス毎に異なる場合も含め）」×「計算に従事するプロセス数」の乱数を
rank0 のプロセスが乱数発生器から取得し，乱数が必要になる Calc.B の手前の Calc.A ルーチンの間に MPI_ISCATTER（あるいは ISEND/IRECV の組み合わせ）（＊）による非同期

integer :: n,nmpi

通信によって各プロセスに乱数を配信します．このやり方

real(kind=8) :: rand(n*nmpi), rand_local(n)では，乱数を利用しないルーチンにかかる時間が，必要な!$OMP PARALLEL

乱数の発生にかかる時間に比べて長ければ，ほぼ待ち時間!$OMP MASTER ! only master thread does the follow-

なしで並列計算が進められるため，計算の高速化にもつな

ing process

がります．各プロセスでは受け取った乱数を１ステップでif (myrank==0) then !assume MPI rank 0 has the

全部利用しない場合もありますが，コードの実装を単純に

QRNG board

するために，利用しなかった乱数は捨てて，次のステップ

call QRNG_GEN(rand,n*nmpi) !generate random

で新しい乱数を rank0 から受けとる形にします．この方法

numbers

を私は「源泉かけ流し方式」と名付けました．温泉で源泉

end ifから汲んだお湯を，下流の各浴槽の使用量にかかわらず常

call MPI_SCATTER(rand,n,mpi_real8,rand_local,n, &に溢れるぐらい十分な量を流し込み続け，使わなかったお&湯は再利用せずに捨ててしまう方式をイメージしたもので

mpi_real8,0,mpi_comm_world,ierr)!scatter random numbers ( n per each MPI rank)

す．Master-Slave 形式の並列計算（パラメータ並列計算!$ OMP END MASTER

で，Master プロセスが計算全体の管理を行い，手の空いた!$ OMP DO schedule (dynamic)

Slave プロセスに次の計算パラメータを順次伝える形式）!thread parallel

do j=1,mで初期条件に乱数が必要な場合にも，この源泉かけ流し方…. ! (Calculation A)

式で Slave プロセスのどれかが計算を終えるまでに，Mas-

end doter プロセスが次に使う乱数をため込んで準備しておく，!$ OMP END DOという形での利用が可能です．また，乱数発生法が物理乱!$ OMP END PARALLEL

数に限定されるわけでもなく，疑似乱数を使った並列計算!... Both communication and calculation A are com-

をしたいが，先に説明した Jump Method を準備するのが面

pleted here.

倒という場合にも，一つの MT 法で発生させた疑似乱数列

call sub_calc_B (n, a, b, c, •••, rand_local)

を源泉かけ流し方式で各プロセスに分配して使うことも可! Calc B which uses the random numbers

能です．

call sub_calc_C (.....)

次に MPI×OpenMP のハイブリッド並列の場合を考え•••ます．もし図２の Calc.A の計算時間が，QRNG で全 MPI

各 MPI プロセスには，0〜k-1 の k 個のスレッドがありまプロセスが利用する乱数を発生する時間より十分長い場

す．このうちの０番スレッドが Master スレッドであり，合，rank0 を乱数生成だけに専念させるのは CPU コアを無Master スレッドだけが，$OMP MASTER 〜 END MASTER の区間

駄に遊ばせていることになりますので，rank0 には計算にを実行します．Rank0 の MPI プロセスの Master スレッドはも従事してもらいたいです．そこで，Calc.A の並列実行の乱数を発生器から取得し，MPI_SCATTER で n 個ずつ各プロセ

仕方を図３のように変更してみましょう．この計算のプロ

ス に 分 配 し ま す．一 方，Rank0 以 外 の MPI プ ロ セ ス のグラムのイメージは以下のようになります（GitHub［８］にMaster スレッドは，乱数の受信に携わります．その間，も gensen_hybrid_img.f90 として上げてあります）．

全ての MPI プロセスにおいて，Master スレッド以外の（＊）図２では簡単のため MPI_ISCATTER を使っていますが，これでは rank0 自身にも乱数を分配することになります．この無駄を省くには，総プロセス数を nmpi，１プロセス当たりに分配する乱数の数を n とした場合，まず n×(nmpi−1)個の物理乱数を rank0
で生成し，それを配列 rand(n*nmpi)の n+1 番目以降に保存します．配列 rand の先頭の n 個の値は与えなくてよいです．そして，
この配列 rand を MPI_ISCATTER で配信すれば，rank0 自身には無意味なデータが渡り，rank1 以降に生成した乱数が n 個ずつ分配
されます．

３７０

Lecture Note

3. Techniques of Acceleration and Parallelizing the Random Number Generation

S. Satake

は，MT 法の内部状態と水瓶方式の バ ッ フ ァ の 中 身 を
save，load す る た め の サ ブ ル ー チ ン save(load)_grnd_
dsfmt が用意してあります．その利用法のサンプル cont_test.f90 が GitHuB に上げてありますのでご参照ください［９］．は じ め に，input_cont.txt 中 の njob を０に し て
kmrand_cont を MPI 並列で実行すると，インデックスがiseed のジャンプファイルを使って新規に ntest 個の疑似乱数を各 MPI プロセスが生成し，”out_r**R_s**I_i000”に書き出します．ここで**R は MPI のランク，**I は選択した
ジャンプファイルの seed 番号です．また，乱数の継続用のファイル ”rand_cont_s**I_i000” ，”grnd_buff_ s**I_i000”も作られます．これらはそれぞれ，MT 法の内部状態と水瓶方式のバッファの状態を全 MPI ランク分まとめて書き出したものになって い ま す．次 に，njob=1 と し て 再 度
kmrand_cont を実行すると，上の継続用ファイル２つを読

み込み，疑似乱数列の続きを ntest 個作ります．以降，njob=2,3,…と す る ご と に，”rand_cont_s**I_i**N”,”grnd_
buff_ s**I_i**N” (N=njob-1)のデータを読んで疑似乱数列を一つ前のジョブの続きから生成します．ntest の数をいろいろ変えて実行し，各 MPI ランクが作る疑似乱数がそ
図３

れぞれ異なる系列の乱数列で，かつ各々はジョブの区切り

MPI-OpenMP ハイブリッドコードにおける源泉かけ流し
方式の模式図．各プロセスの Master スレッドは，乱数の送受信が終わったところから Calc.A のスレッド並列計算に途
中参加します．

を跨いで連続した乱数列になることを確認してください．
駆け足になりましたが，これで KMATH_RANDOM ライブラリを利用した，並列化 Mersenne-Twister 法の利用と高速化

1〜k-1 番スレッドは，Calc. A の Do loop を並列実行し始と，物理乱数発生器を並列プログラムで利用する方法の解めています．$OMP DO に SCHEDULE (dynamic)が指定されて説はおしまいです．本章の内容とサンプルプログラムが，いるため，Master スレ ッ ド が 乱 数 の 生 成 と 分 配 を終え乱数を利用した多並列シミュレーションコードを書く方々て，まだ他のスレッドが Calc.A を終えていなければ，図３の一助になれば幸いです．

に示すように Master スレッドは途中から Calc. A の並列計算に参加することになります．このようにすれば，源泉か

参考文献

け流し方式でも乱数発生器を持つプロセスも計算に参加［１］講座「核融合プラズマシミュレーションの技法 ５．粒子シミュレーションのコーディング技法」プラズマ・核融合学会誌 89, 245 (2013).
［２］M. Matsumoto and T. Nishimura, "Dynamic Creation of
Pseudorandom Number Generators", Monte Carlo andQuasi-Monte Carlo Methods 1998, Springer, 2000, pp. 5669.
［３］http://www.math.sci.hiroshima-u.ac.jp/ 〜 m-mat/MT/
DC/dc.html［４］http://www.math.sci.hiroshima-u.ac.jp/ 〜 m-mat/MT/
SFMT/JUMP/index-jp.html［５］http://www.math.sci.hiroshima-u.ac.jp/ 〜 m-mat/MT/
SFMT/index-jp.html［６］H. Haramoto et al., "A Fast Jump Ahead Algorithm for
Linear Recurrences in a Polynomial Space", Lecture
Notes in Computer Science, vol. 5203, Springer, Berlin,
Heidelberg (2008).
［７］https : / / www. r-ccs. riken. jp / labs / lpnctrt / projects /kmath-random/［８］https://www.shoup.net/ntl/index.html［９］https://github.com/satakeshinsuke/random̲number/blob/master/Chapter3.tar

し，かつ乱数の分配のMPI通信時間を Calc.A の演算に隠ぺいすることが可能です．もちろん，この方式が効率よく稼働するには MPI 並列数とスレッド並列数を適切に選び，Calc.A にかかる時間より乱数生成時間が短くなるようにしなければならないので，どんなプログラム・計算環境にも使える万能な方法ではないことに留意してください．

・KMATH̲RANDOM における継続計算の方法
最後に，KMATH_RANDOM を使った計算を複数回継続する方法を説明します．第２章でもみたように，１回のジョブの最後に MT 法の内部状態をファイルに書き出し，継続ジョ
ブの最初にそれを読み直せばよいのですが，MPI 並列プロ
グラムではそれを各プロセスが行う必要があります．更に，grnd_dsfmt のように水瓶方式を使っている場合，「水瓶」の中身も保存し，継続計算に引き継がなければなりま
せん．そうでなければ，例えば１００ステップを３回つないだ乱数を使った計算と，１５０ステップを２回つないだ計算とで途中から利用する乱数がずれてしまいます．
それでは，サンプルプログラムを使って，継続の仕方を具体的に説明していきましょう．km_rand_wrapper.f90 に
３７１

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).

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

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

さ

たけ

しん

すけ

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

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

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

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

３７８

J. Plasma Fusion Res. Vol.96, No.8 (2020)４
０５‐４１７

講座 シミュレーションのための乱数入門

５．乱数の変換と利用例
5. Transformation of Random Numbers and its Applications
佐 竹 真 介，菅 野 龍 太 郎，本 多充１）

SATAKE Shinsuke, KANNO Ryutaro and HONDA Mitsuru１）
自然科学研究機構核融合科学研究所，１）量子科学技術研究開発機構（原稿受付：２
０
２
０年５月２
２日）

この章では，プラズマ・核融合分野，特に，磁場によって閉じ込められた核融合プラズマの理工学におけるシミュレーションにおいて，乱数がどのように利用されているのかを紹介します．シミュレーションにおいては，一様乱数だけでなく，正規乱数など様々な乱数が目的に応じて利用されています．そこで，まずは，乱数を用いたシミュレーションを行う際に必須である，様々な乱数を生成する「乱数の変換」について説明します．次に，シミュレーションにおいて，具体的にどのように乱数が用いられているのかを解説します．乱数を用いるシミュ
レーション手法は，一般に，モンテカルロ法と呼ばれます．ここでは，モンテカルロ法に基づくシミュレーショ
ンについて，基礎的な例を挙げて説明します．また，近年，様々な科学技術分野において注目を集めている遺伝的アルゴリズムと，そのアルゴリズムを用いた大域的最適解を求めるシミュレーションについても説明します．
Keywords:
random number transformation method, Monte-Carlo method, Dirichlet-Poisson solver, stochastic differential equation,
genetic algorithm, global optimizationに持つ単調増加関数になります（ 'が不連続な場合につい５．
１ 乱数の変換ては後述）．よって，,は !の逆関数 !!"を使って

本講座の第２章から４章で紹介してきた疑似乱数や物理乱数は基本的に +!
!
".の実数値一様乱数でした（利用((-+*,#!!")

するサブルーチンによって )
，!
!
!
!
"*
!
"*ないし )!
".の場
合もありますが，ここでは厳密な区別はしません）．しかのように表すことができます．さて，今 #が区間 !
!
".

し，様々なシミュレーションに乱数を利用する際には，一に一様分布する確率変数とすると，#$+となる確率は様乱数以外の様々な確率密度関数に従う乱数が必要になり#$+,,*で す．こ の 時，こ れ に 対 応 し て，
$&
+
#+#!)

ます．本節では，そのような任意の確率密度関数に従う疑,*で$#!!")#*+*と な る 確 率 も 同 じ く !)$,#!!")

似乱数を発生させる方法を紹介します．

す．すなわち，

任意の確率密度関数に従う乱数を生成する方法として

!!")
#*
#"!)
,*
,*
$&
+
$,,
#$&
+
,
#!)

は，大きく分けて２つの代表的な方法が知られています．
まず１つ目は，逆変換法（直接法）と呼ばれる方法で，確が成り立ちます．上式より，区間 !
+!
".の一様乱数 +(,率密度関数 'に対し累積分布関数!"
を用いて ,,+(,は累積分布関数(#! )(*として与えた +,!),*に従う分布をとる，つまり確率密度関数 'に従う分布%になることがわかります．

+#!)
,*
*
%*
#! '
)
*
を考えることで，一様乱数 ++,から !を通じて 'に従う乱(

ところで，!!")+*は常に解析的に表せるわけではなく，数+,(,を生成するという方法です．以下，具体的に見てい不連続な関数である場合も考えられます．そのような場合

きましょう．

は，生成したい乱数+の値域,%!
&.を十分細かく分割し，(,数値積分を使うなど し て +
（)#!
!
!
!
"），
"
#
/!
,
) #!)
)*

ここでは簡単のため，確率密度関数の定義域が %!
&.の

+
，+
，の離散データのテーブルを用意
!#!
" #"

連続な関数であるとします．もし下限，上限の両方あるい

ただしはいずれかが !'，"' となる場合は，*"%，*#&で十分します．そして，一様乱数ルーチンで生成した +
(の値が

に'
*
)
*
&!とみなせるような下限，上限を設定しておきま+) $+)""に入る時，例えば線形補間を用いて，("+しょう．さて，確率密度関数の性質として，'
*
)
*
%!である

,
)""!,
)
,
+
)
)"
)*
(#,
(!+
+
)""!+
)

こと，定義域全体で積分すると "になることの２つがあります．したがって，累積分布関数 !),*は，値域を !
!
".
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

Journal of Plasma and Fusion Research Vol.96, No.8 August 2020!"を生のようにして 4の近似値を得ることで+
4
3
(,
($! )
(*

/

#
0
!/

#
!/

1$!)
/
*
$% '/0/0& # $"!& #成することが可能です．

!

なお，+
4
(,が離散値しか取らない場合も，上と同様の考えか ら 逆 変 換 法 を 用 い て 生 成 で き ま す．確 率 変 数 4がの値域が !
!
".なので先ほどの逆変換法が使えそうです．

+ 通りの離散値 4!
!
+）を取り，それぞれの発
#
1!
*（*$"

しかしこのままでは逆関数が簡単に記述できないので，生確率が -* であるとします（ただし # +
）．この
* $"
*$"-

10$"!1と置き換え，10$& # とすると，その逆関数は/$&!#)
10*で あ る こ と が わ か り ま す．と こ ろ で
*)
/
!
'*
3
!
4*の逆変換は)/)#!/

時，区間 !
!
".一様乱数 3(に対して，
!"
,
#,
* "3
*
(%# *$!*$!-

3$/&
+,'，

を満たす,を求め，この,に対して 4,を採用すること
($4

4$/,
(
*'

で与えられるので，一様乱数 +1(!
2(,を用いて，で，+4(,を生成します．ただし，上式の総和を取る時に便

，
3
1(*&
&2(*
*)
+,)
#
!#)
($&
，
4
1(*,
&2(*
*)
(
*)
#
!#)
($&

宜的に -!$!とおきます．
もう１つの代表的な方法は棄却法と呼ばれるもので，一様乱数を使った試行錯誤によって得たい分布に従う乱数をと与えると，+
は２次元空間の正規分布 '
3
3
!
4*に従う
4
)
(!
(,

発生させるものです．まず，先ほどと同様に，生成したい分布を持つ乱数となります．

，その確率密度関数を '
とし，さら
乱数 +
#!
$.
.(,の範囲が-他にも様々な関数形の確率分布を一様乱数から生成するに'0"%であることが既知であるとします．ここで，%)*方法が知られていますが，ここでは紹介しきれないので，は厳密な最大値である必要はなく，確実に '
0
"%となる
)
*

興味のある方は文献［１，２］などを参照してください．

ような適当な定数を設定すれば大丈夫です（近い方が乱数（佐竹）の生成効率はよいですが）．次に，一様乱数発生ルーチンを使って，二組の一様乱数 ()(#!
$.と%!
%.を用意し!
)(-５．
２ 乱数の利用例ます．これに対して，'()*と %))の大小関係を判定します．

ここからは，磁場によって閉じ込められた核融合プラズ

もし '
$%
()*
)
($(
)を採用し，そうでなければ)であれば 4

マの理工学における乱数を用いたシミュレーション手法をに対して同様の判定を行います．そし
破棄し，次の)()!
%)*

いくつか紹介します．シミュレーションにおいてどのようて，必要な個数の +4(,が採用されるまでこの判定を繰り返に乱数が用いられているのかを説明する前に，なぜシミュせば，確率密度関数を 'に従う乱数 +4(,が生成されます．

レーションで乱数を用いるのかについて考えてみます．

!0

この他にも，'
0)*$& や正規分布の場合など，特定の関まずは，シミュレーションによって扱われるプラズマの数の場合に使える乱数の変換法が色々と知られています．

諸現象が，どのような方程式により記述されているのか，

ここでは代表的なものとして，２次元空間における正規分

その概略を見ていきましょう．議論を単純化するために，

布'
3!
4*を与える，Box-Muller 法を紹介します．
)

完全電離プラズマが，それぞれ "##個の水素イオンと電子

Box-Muller 法では区間 !
!
".一様乱数を二組用意しまで構成されているとします．ここで，"&"!#"程度です．

す．それらを +1(!
2(,とします．この手法では 1(と 2(は無このとき，荷電粒子の運動は，Newton 方程式（および電磁相関で独立な一様乱数であることを前提としますが，第４場を解くために Maxwell 方程式）で記述することができま章で説明したような「品質の良い」疑似乱数であれば，実す［３‐５］．電場と磁場をそれぞれ "および !とし，)番目の際に独立な疑似乱数列を２系統使う必要はなく，１系統の荷電粒子の質量を +)，電荷を .)として，Newton 方程式か

疑似乱数列からその奇数番目を 1(，偶数番目を 2(，というら粒子の位置 $))00*および速度 #))*は，具合に取れば十分です．さて，平均が )，分散共分散行
!
!
!*

/
$
0
0
*
$#))
*
))

列が !" !"となる２次元正規分布関数は
! "

#

および

#

" !3 "4
')
3
!
4*
$ & #
&
#

/
#
0
*
$
))

で 与 え ら れ ま す が，こ れ を 極 座 標 系（/$&3#"4#，）で表現すると，
4#
3*
'$%
*!")

.)
")
0
!$))
0
0
0
!$))
0
+
*
*
"#))
*
#!)
*
*
,
+)

で与えられます．ここで，)$"!
!
!
" です．この New#$1!
ton 方程式と Maxwell 方程式を解くことができれば，原理

!/#

/
')
/
!
'*
$ &#
&
#

的には，プラズマの温度，密度，閉じ込め時間など，磁場閉じ込め核融合炉のプラズマ性能を評価するために必要なであり，'に関しては !
&*に一様に分布していることに!
#なります（ここで因子 /は )
3
!
4*
/
!
'*の座標変換のヤコ/)情報を得ることができるはずです．しかし，現実には，この方程式系を解くことは不可能です．なぜなら，全粒子間'ビアンから出てきます）．また，$//&!/##$"を満たす! 'の相互作用を考慮して，" 個すべての荷電粒子の Newton

ことは容易に確認できます．そこで，/方向に関する分布

方程式を解く必要があるからです．

#

を乱数から生成する方法を考えると，累積分布関数

そこで，プラズマの状態を近似的に記述する基礎方程式
４０６

Lecture Note

5. Transformation of Random Numbers and its Applications

S. Satake et al.

としては，水素イオンと電子のそれぞれに対して，２体相に近いと仮定しています．一方，プラズマの力学的平衡状

互作用のみを考慮した６次元位相空間内の１体粒子分布関態における磁場や，温度・密度分布そのものは，計算コス数%%!
$)の時間発展方程式である Boltzmann 方程式：)!
#(

トを考えて，通常，Boltzmann 方程式から導かれる３次元空間の流体方程式とMaxwell方程式を自己無撞着に解くこ$$ (# " )
$
*)
)
!
*)
%
(!
"$#!(
)
$ "%
! "$$ " (# %"#(
#)
$
)
$% &#
$$

とで求めます［１２］．この流体方程式に基づいて温度・密度分布を決定する輸送方程式が非線形偏微分方程式であるたが用いられます．ここで，#は粒子種を表し，"および !

めに，解を数値的安定に求めることが困難になる場合があは外場とプラズマの集団運動が作り出す電磁場とします．

ります［１３］．そのような場合に対処するため，５．
２．
２節で$個の Newton 方程式から Boltzmann 方程式を導出する議紹介する遺伝的アルゴリズムを用いた手法が開発されてい論については，文 献［３‐５］を 参 照 く だ さ い．上 式 の右辺ます．この場合，乱数は，物理現象のモデル化に用いるとは，２体相互作用（Coulomb 衝突）によって起こる速度空

いうより，非線形偏微分方程式を解くアルゴリズムにおけ間での拡散を表現する Landau 衝突項で，Fokker-Planck

るツールとして使用されているわけです．

型の式（＊１）として記述されます［３］．さらに，左辺にお以上のように，乱数は，運動論方程式を解く際の膨大ないて，荷電粒子のジャイロ運動に対して近似的な取り扱い

計算コストを減じるためのモデリングにおいて利用されたを行うことで，５次元位相空間の１体粒子分布関数の運動

り，流体方程式により記述される非線形な輸送方程式を解

論方程式（近似の取り方によって，ジャイロ運動論方程式，

くために利用されたりしています．

またはドリフ ト 運動 論 方 程 式）を 得 る こ と が で きます（菅野・本多）５．
２．
１ 拡散過程とその応用［５，６］．この運動論方程式によって，現実的にプラズマの諸現象を扱うことができるようになります．

本節では，乱数を用いた基礎的なモンテカルロシミュ

プラズマの振る舞いは，運動論方程式によって記述されレーションの代表的な例として，拡散方程式の解法とるわけですが，次に，シミュレーションにおいてなぜ乱数

Dirichlet-Poisson 混合問題への適用の２つを取り上げ，乱が用いられるのかを Landau 衝突項に注目して考えていき

数がどのように利用されているのか紹介します．紹介するます．例えば，プラズマ性能を決める主な要因の１つであ

シミュレーション手法の数学的基盤は，確率過程論です．

る微視的乱流輸送を扱うシミュレーションにおいては，運「確率過程」とは，時間と共に変化する確率変数（確率的に動論方程式の主に左辺の寄与により生じる１体粒子分布関値が定まる変数）です．この後に示すように，シミュレー

数の微細揺動に対する数値計算の高精度化のために，分布ションにおける乱数の役割は，数値計算上の単なるツール

関数そのもの を５次 元 の 連 続 体 と し て 扱 う こ と が多くではなく，確率過程論に基づいて拡散現象を粒子的にモデ［７］，Landau 衝 突 項 は，速 度 空 間 の 微 積 分 方 程 式 と しリングしたときに導かれる数学的性質の具現化を担っていて，差分スキームなどを使って数値的に解きます．このよます．

うな微視的乱流シミュレーションは，理化学研究所のスー

以下では，確率過程論における数学的結果（定理など）

パーコンピュータ「京」などの大規模計算機を利用して１００を用いた手法の説明に重点を置きますので，数学用語の厳

時間程度の計算時間を必要とします［８，９］．一方，最近の研密な定義や数学的結果の導出については，文献［１４‐１７］な究［１０］で，プラズマ性能を劣化させる不純物の輸送においどを参照ください．なお，本節で扱う手法では，簡単のたて，Coulomb 衝突の寄与による輸送現象（新古典輸送）のめ，空間を１次元としますが，'次元空間（'""）へ拡張

重要性が指摘されていますが，そのような現象に対するシできます．

ミュレーションにおいては，微視的乱流シミュレーション!拡散方程式の解法

ほどの速度空間分解能は要求されていないので，計算コス

プラズマ・核融合分野におけるモンテカルロシミュレー

トの低減を優先した粒子的なシミュレーション手法を用いションの代表的なものの１つは，拡散方程式の解法ではなることが多いです［１１］．その場合，Landau 衝突項によるいでしょうか．ここでは，荷電粒子の Coulomb 衝突や，熱拡散現象は，分布を構成するテスト粒子の速度空間におけ

拡散などのランダムネスが内在する物理現象に対する粒子

るランダムウォークによってモデル化されます．そのよう的なモデリングに基づいたシミュレーション手法を紹介しな Coulomb 衝突の効果を考慮したシミュレーションにおます．以下のような拡散方程式（Fokker-Planck 方程式）をいて，粒子のランダムウォークを具現化するため，乱数が解くことを考えることにします．

用いられるのです．このような乱数の利用に対する簡単な$%)$" $# #(
*)
!(
*)
%
)
!
*)
%
)
!
*)
(!
%! *
(
+
"
* *)
(
+
$
)
$
*
#$
*#
（１）

例を５．
２．
１節で紹介します．また，このシミュレーション手法を発展させると，科学技術分野で扱われることが多い
Dirichlet-Poisson 混合問題に適用することができ，そのよ

ここに，時間パラメータを )&!，空間を *とし，)!
*'#

うな応用についても紹介します．
上記のような運動論方程式を解く際には，多くの場合，で，#は 実 数 と し ま す．ま た，!(
，#(
*)
*)
'#お よ びプラズマを閉じ込めている磁場や，プラズマの温度および#(*)&!です．式（１）を解くため，分布 %を構成する粒子密度の分布は，既知である，もしくは，初期に与えたものの位置 *()))という量を導入し，*()の時間変化 の式を（＊１） Fokker-Planck 型とは，どのようなものかについては，５．
２．
１節で紹介します．

４０７

Journal of Plasma and Fusion Research Vol.96, No.8 August 2020

＋ノイズ項，と表現しましょう．ここ
+
#
-*
+
%-*
%+#!*
+
+
+

具体的な例で，拡散方程式を解くモンテカルロシミュ

で，この「ノイズ項」から生じるランダムウォークの歩幅レーションがどのようなものかを見ていきましょう．プラが拡散係数 #*-+を与えるように設定します．多数の粒子

ズマ中のテスト粒子（例えば，電子とします）が，背景プに対して，この式を解くと，時刻 +におけるそれぞれの粒

ラズマ（水素イオン と し ま す）と の Coulomb 衝 突 に よ子の位置 -*++が求まりますので，粒子の位置分布も得らり，速度空間をピッチ角散乱する場合を考えます．電子のれ，それが式（１）の解です．数学的に，もう少し厳密に言速さを ,，磁力線方向の速度を ,/として，*#,/#,で速度空えば，式（１）の解 &+!
-+を求めるために，式（１）に対する*間を表現することにします．*は，電子のピッチ角 &（速度

粒子描像の微分方程式：*+&を表します．ここベクトルと磁力線のなす角）の余弦 $

+
-*
+
-*
+
+
%-*
%+
%%*
+
#!*
+
+
")*
+
+
+

で，ピッチ角散乱によって速さ ,は変化しないと仮定して（２）います．*の取り得る値は，!"%*%"です．電子の速度空に注目し，この式を解いて得られる「確率過程 -*」の分++

間 *における分布 &を与える拡散方程式は，以下のように

布が解 &
+
!
-+であること（＊２）を利用します．ここで，*表現されます［１９］．

-+)*-+です．この確率過程 -*+#*-++は，拡散過程と#)*++&+" +# ( " # &
*+
*
&
*!
#! 0
!(
0 , !* -1 （６）
1
"
&'
+
*
+
+
#+
*# &'

も呼ばれます．式（２）の右辺第２項における%*++は，数学における意味での「Brown 運動」で，Wiener 過程とも呼ばれます［１４］．Wiener 過程 %*++は，①連続で，②離散時刻

簡単のため，ここでは，電子‐イオン衝突周波数 (
&'を定数

!#+
!+
#!
!"+
""+
#"4"+
(!""+
( を 適 当 に 取 っ て，%*と 仮 定 し ま す．式（６）に 対 応 す る 拡 散 過 程 **
，つ ま
+
+

として，任意の時刻 +++!%*( に対して %*(+(!"+が平均 !

り，テスト粒子の速度空間における位置の時間発展は，時

で分散 *
+
!
+
!
(!+
(!"+の正規分布 $*(!+(!"+に従い，③ま間ステップを %+として，式（２）のように与えられる確率

と%*
が互い
た，()
#'として%*
+
+
+
+
!%*
!%*
(+
(!"+
!"+
'+
'

微分方程式から自然に導かれる計算スキーム：に独立となるという性質を持ちます［１４］．式（２）の「ノイ

+
**
+
%+
%% （７）
**
+
+
"!*#*
+
+
"#(
"%+
+
#**
+
!(
&'
&',

ズ項」において %%*+#+%+のように書かない理由は，%*++

により与えることができます．ここで，%**
+
+
+
+
'**
"%+
+
!**
+

が至るところ微分不可能であるためです．
式（２）のような方程式は，確率微分方程式と呼ばれ，日としました．また，テスト粒子は，時刻 +#!で位置 **!+本の数学者，伊藤 清氏によって，１９４２年に数学的基礎が与から出発するとします．式（７）のような計算法は，Euler‐

えられました［１８］．式（２）は，確率過程 -*++を求めるた丸山スキームと呼ばれています［２０，２１］．もちろん，常微分

めの「積分」を表現していて，方程式の解を求める数値計算スキームと同様に，Runge-

+

+

!

!

Kutta スキーム（後で紹介します）など高次の近似法もあ

-*
*
-*
+
-*
*
*
%*
!+
%%*
+
+
"" )*
+
#-*
"" !*
+
+
+ （３）

ります［２０‐２２］．ここでは，Euler‐丸山スキームを利用することにします．

のように与えられます．式（３）の右辺第３項は，伊藤積分

式（７）のような計算スキームは，確率微分方程式の時間と呼ばれ，)*-**++における時刻 *の取り方に気を付けて，離散近似です．時間離散近似には，大別して，確率微分方+

程式の解の経路 ,**+++&!をよく近似する「強い近似」と解

)!"

%*
+
+
-*
+
+
,
!%*
" )*
-*
*
*
%%*
'
) ! )*
+
+
+
#(
(+
(""+
(+
%+3 !(#!
!
（４）の平均や分散などの分布特性を近似する「弱い近似」があり［２０］，%% に対して，強い近似では正規乱数 'を用いて%% #'%
2 +，弱い近似では，例えば，%% #$2%+のように取ります．ここで，弱い近似における $ は，コイン投げなのように定義されます．後で数値計算することを意識して，%+#+
%+，%*
#
)とし，+
+
+
#%*
"%% とし
( #(
(""+
(+

どで発生させた，値が ""または !"の二値乱数を用いて，

ま す．こ こ で，%% は，平 均０で 分 散１の 正 規 分 布

それぞれ確率 1/2 で与えられるとします．
（弱い近似にお$*!
!
"+に従う正規乱数'を用いて，%% #'%2 +と表現できます．確率微分方程式（２）の解の存在とその一意性のた

いては，%% に対して他の取り方もあります．
）どちらの

め，!*
-+および )*-+の満たすべき条件は，適当な定数切に判断すればよいでしょう．

近似を採用するかは，この後に示すように，問題ごとに適"!および "" に対して
!*
-+
)*
-+
-.
"".
.
.
".
.
%"!,
-

拡散方程式
（６）を数値シミュレーションで解く場合には，Euler‐丸山スキーム（７）における強い近似の採用は，（５a）あまり有効ではないことがわかっています．以下では，このことについて説明します．式（７）において強い近似を採

および
!*
-+
.+
)*
-+
.+
-!..
.
!!*
.
".
!)*
.
%"".

用すると，数値的に .
*.
$"となる不具合が生じます．すな（５b）

わち，数値計算において %+は有限な値であるので，１スです．ここで，１５］など!
.(!です．詳しくは，文献［１４，テップの計算の際に発生させた正規乱数の絶対値が，あるで論じられていますので，そちらを参照ください．

値'*.
$"が発生します．ここで，'!より大きくなると .
!

（＊２）証明については，例えば，文献
［１
４］
の問題８．
３を参照ください．この後で紹介する，式（１０）
を用いる計算法との関係も説明されています．

４０８

Lecture Note

5. Transformation of Random Numbers and its Applications

は，ある時刻 )における .の値を .&)"!('として，
'
#$&
)
!#

(",
"!,'
&
(&
#!('
#,

S. Satake et al.

は+
!
"!'
!""'
,の範囲に存在することになります．ここ
で，'#*
%)
#
#です［２３］．また，有限個のテスト粒子を用(*いた手法なので，得られた分布に統計的な誤差は生じま（８）

す．統計誤差は，テスト粒子数 #.に対し，"##- .に比例しで す．た だ し，!",#*%)""，!"("#で，ま た，(*て小さくなります．

.#$&"!('の 正 符 号 は .$"と な る ケ ー ス，負 符 号 は以上が，テスト粒子の位置分布によって拡散方程式の解."!"となるケースで選ぶこととします．.#"!(からを与えるモンテカルロシミュレーションの紹介ですが，こ出 発 す る テ ス ト 粒 子 の１ス テ ッ プ の 計 算 で 不 具 合れとは別の計算法もあります．関数 *&)!
+'を )$!で，（.$"）が発生する確率 $(,,は，$(,,#&
#
)
"
#
(,
)
&&
'
- 'で
!##

*&
)
!
+'
&&
+&
)
#%++
'
'
,

す．出発点を .#!""(とした場合も同じです．ここに，

（１０）

+'
+'で，(,+'は，誤差関数です．発生確率
(,
)
&&
)
)
#"!(,
&
&

であると定義すると，このとき，*&)!
+'に対して以下の関$(,,から，不具合は，*.*#"の近傍で比較的起こりやすいこ

係式が成り立つことを利用します［１４］．

とがわかります（＊３）．時間ステップ %)を小さくするこ/*#!*および *&
!
+'
+'
!
#&&
/
)

とで，１ス テ ッ プ 毎 の 発 生 確 率 $(,,は 小 さ く で き ま すが，ゼロにはならないので，根本的に解決したとは言えま（つまり，せん．数学的に %). !の極限を取れば，$(,,. !

（１１）

ここで，演算子!は，式（２）に対応して，式（７）を確率微分方程式として見たときには測度ゼロのイ#+'/"""&
+'/#
! #!&
/
+ #
/
+

ベント）となりますが，シミュレーションにおいてはゼロでない確率で発生するレアイベントとなります．そのた（１２）ありますし，どのような処理でも，拡散方程式（６）には含であり，%++&&+&)'',は，式（２）に従って時刻 (#!に位置+から出発した確率過程 (
+&
(
+&
)
'
)
'
'
(%!の分布による &&め，*.*$"となった場合の処理をどうするか考える必要がまれていない操作なので，*.*#"近傍のテスト粒子の分布の期待値で，この場合には，式（２）が与える時刻 )におけへの影響が心配です．また，拡散方程式（６）の基となるる +か ら ,へ の 遷 移 確 率 密 度 を '&)!
+!
,'と し て，Landau 衝突項のモデリングでは，小角度散乱を前提とし

!&&
,'
'&
)
!
+
!
,'
',と表現することもできます．最初に紹介ていますので［３，１９］，有限な%)の場合に，%% #)%- )におした計算法と式（１０）を用いる計算法は，Kolmogorov の前ける )の値によって散乱効果が大きくなりすぎること自進方程式と後退方程式の関係になっています［１４］．

体，レアイベントとは言え，モデリングの前提に反するの

先ほどのピッチ角散乱の問題に適用してみましょう．解ではないかという疑問もあります．

くべき方程式は以上の議論から，この問題では，Euler‐丸山スキームに

/#&
/
&
/&)
.'
. ""*
"!.#) #
&!
#!*
(
(*
/
. # (*/)/.

おける強い近似の採用は，あまり有効ではないことがわかりました．そこで，文献［２３］に従い，弱い近似で扱うことにし，%% #$-%)として，（１３）

で，対応する拡散過程は，式（７）（数値計算する際には，
!
.'
.'
!
&
#&&
#
式（９））で す．例 え ば，初 期 分 布 を &

.&
)
)
)
*
)
*
%)
%) （９）
"!.#&
"%)
'
#.&
'
!.&
'
'
)
$"(
(*
(*

+.'
#
#!*
+&
$で あ る と す る と，各 時 刻 )に お け る 解()

例えば，１ステップごとに !から "の範囲の一様乱数を

&
)
!
.'
&&
.&
)
&
#%.+
'
'
,は，図１のようになります．ここで，速度空間の位置 .における時刻 )の分布の値 &
)
!
.'は，式
&

のように計算することにします．ここで，$については，##
Tausworthe（トーズワース）法で発生させ，その値が "（９）により与えられる，時刻 )#!に .から出発したテス未満なら負符号，1/2 以上なら正符号とします．式（９）のト粒子の時刻 )での位置 .&).&)'における &&''の平均値でように計算すれば，*.*$"となる不具合は起こりません．

す．当然のことながら，図１の分布は，時刻 )#!でテスト

時刻 )#!で初期分布 &!
.'となるようにテスト粒子の初!
&粒子を初期分布 &&（９）によっ.'となるように配置し，式期位置 .&（９）によって与えられる時!'を設定すれば，式て時間発展させたテスト粒子の位置の分布と同じものにな

刻 )におけるテスト粒子の位置の分布が，拡散方程式（６）ります．最初に紹介した計算法の方が，シンプルでわかりの数値解 &)!
.'になります．時間ステップ%)を十分小さく&やすく感じるかもしれませんが，式（１０）を用いる計算法の取れば，分布の計算に関して，式（９）は，確率微分方程式

考え方は，より一般的な初期値境界値問題への応用に発展として見たときの式（７）の良い近似となり，シミュレー

させることができ，次に紹介する手法とも関連していまションにおいて，分布 &)!
.'は，初期分布 &
!
.'に依ら
!
&
&

す．詳しくは，例えば，文献［１５］で論じられていますので，!"（#*ず，１衝突時間 '以降に理論通り一様な分布に緩(*(*そちらを参照ください．
!Dirichlet-Poisson 混合問題への適用

和します．ただし，どのような初期分布であっても，.&)'

この節の後半では，科学技術分野で扱われることが多いの数値計算における最初の１ステップ以降は，テスト粒子（＊３）出発点を *（５b）を満たしていないため，解の一意性に問題が#"とした場合は，*.*.*$"となる不具合は発生しませんが，式あるように見えます．しかしながら，後で示す式（９）に従って，有限な時間ステップ %)で数値的に解く際には，右辺第２項がある
ために，そのような問題は生じません．

４０９

Journal of Plasma and Fusion Research Vol.96, No.8 August 2020

り，Runge-Kutta スキームを採用することで，これらを与える式（１６）に対する計算精度を向上させることができます．ここで，式（１６）に対する Runge-Kutta スキームは，以下の通りです［２１，
２２］．
" !
"
-'
+
+
)"")#(
%
"#+
(
$-'
(
" #+
" -#+'' ""!#(
"
#
#
（１８a）
!"$!'
-'
+
(
(
（１８b）

図１

ピッチ角散乱による電子分布 f (t, +)の時間発展を式（１０）に基づき計算．このシミュレーションにおいて，'ei#t ＝ 10−3と設定し，−1!+!1 の範囲の格子点（分割数５０）のそれぞれに対して，105 個のテスト粒子を用いて期待値 E+[ $(+(t))]
を求めました．実線は，初期分布 $(+) ＝ { 2 − sin ((+)}/4，記号＋は時刻 t/*eit ＝ 0.5 に お け る 分 布，記 号!は 時 刻
ei ＝ '
t/*ei ＝ 5.0 における分布です．

（１８c）

!#$!'
-'
+
%
(
"!"#+
")"-#+
#(
-'
+
%
)#$)'
(
"!"#+
")"-#+
#(

（１８d）
（１８e）

ただし，%
!
!
"(に従う互いに独立な正規乱数
"と %
#は，%'
です．これらの正規乱数は，例えば，Tausworthe 法で発生させた !から "の範囲の一様乱数を使って，
５．
１節で説明した Box-Muller 法で生成します．
以下の例題に適用してみましょう．
*#'*'-#""-" #!''%"-( !#'""-(!-$"'*# *-'"!-($&!#

Dirichlet-Poisson 混合問題に対して，モンテカルロ法をいかに適用するか紹介します．次の方程式を考えます．
!"&'
-(
,'
-(
-(
)
*
$#'

)"$)'
-'
+
(
(

（１４）

（１９）

!
），境 界 条 件 はた だ し，-&+!
!
",と し て（!$!
"$"
$'
!(
"(
$"お よ び $'$!と し ま す．式（１９）の 解 は，

ただし，演算子!および境界条件は，以下の通りとします．

'
-(
-"
-#"
#(
#(です．シミュレーションの結果を
'
$"!'
!'

#
,
!,$ ""'
-(,,
-(,
"!'
,
#
,
-#

図２に示します．ここで，図２では，得られた数値解の平（１５a）

滑化のため，最小二乗法を用いた多項式近似を行っていま

境界条件：-$-(-( （１５b）'$$'!!
" において，,す．多項式の項数の決定には，AIC（Akaike Information

ここで，空間 -については，-&+!!
",を仮定しています．

Criterion，赤池情報量規準）［２４］を利用しました．AIC は，式（１５a）から，これに対応する確率微分方程式は，次のよ

統計モデルの良さを評価し，良いモデルを選択するための

うになります［１４，１５］．

指標です．数値的に得た解 'を微分したい場合には，その

+
-'
+
+
-'
+
*-'
*&'
*+
(
$)'
(
(
(
"!'
(
(

ままでは数値解の統計誤差がノイズになるので，AIC など（１６）を利用しながら数値解を平滑化することが有効です．また だ し，"'-(-()'-(で す．方 程 式（１４）の 解 は，式$)'た，得られた数値解の妥当性は，異なる乱数列を使って複（１６）に従う確率過程 -'+(により数回，数値解を求めて確認します．例えば，２回計算して得られた数値解をそれぞれ '"，'#として，以下のような相

*

-

-'
*
-'
*
,'
-(
)'
+-,%' &'
**
(
(
(
(
$( #
&$
!

*

*

!

!

' #'
-'
*
-'
)
+-,%' &'
*)&**
(
(
(
(
!(-#
（１７）
$
と与えられます［１４，１５］．ここで，(-+#,は，式（１６）に従って，時刻 +$!に位置 -か ら 出 発 し た 確 率 過 程 )
-'
+
(
*
+%!
の分布による期待値を意味します．また，*は，式（１６）に従って動くテスト粒子（モンテカルロ法におけるサンプルの外へ最初に出た時粒子）それぞれに対して，区間 +!!
",刻によって与えられます．式（１７）により数値解,'-(を求めるには，すべてのテスト粒子それぞれに対する時刻 *の値が必要となるので，式（１６）の !'-(の値などから，計算時間がどれくらいか推定できても，乱数を用いてシミュレー

図２

ションしているために，計算時間を正確に定めることはできないという欠点があることに注意ください．
式（１７）を 用 い る 計 算 法 で は，各 テ ス ト 粒 子 の 経 路
-'
+
)
*
#!
$!
. と境界に最初に到達した時刻 *が重要であ
((
($"!
４１０

方程式（１
９）
の解析解と数値解．解析解を実線，数値解を破線で示しています［２５］．数値解は，式（１７）に基づき計算．
区間[0, 1]における格子点（分割数１００）のそれぞれに対して
103 個のテスト粒子を用いて，数値解を求めました．また，#t ＝ 10−4 と設定．得られた数値解の相対誤差は，-＜1/200
でした．

Lecture Note

5. Transformation of Random Numbers and its Applications

S. Satake et al.

敗し，!%
%&
%&かつ "#!とすれば解析解と良く一
#!###
%

対誤差
%
"

%
"

%
!

%
!

#
#
$#! '
%%! '
%%
"!#
#'
"'

致する数値解が得られました．もちろん，例題（１９）においても，左辺第３項を右辺に持っていき，ソース項として（２０）扱って問題無く解くことができます［２５］．
以上のように，本節では，確率過程論に基礎づけられたが十分小さいなら，数値解として妥当と言えます．
例題（１９）では解がわかっていましたが，もちろん，事前

モンテカルロシミュレーションを紹介しました．プラズに解がわからない問題を解くのが通常です．第４章で議論

マ・核融合分野では，ここで紹介したもの以外に，いろい

した，乱数の統計的品質が良いと判定される疑似乱数を利ろなモンテカルロ手法が開発・利用されています．例え用していれば，疑似乱数の性質が計算結果に影響すること

ば，文献［２６‐４４］などです．文献［２６‐３１］は Coulomb 衝突のは一般的には考えづらいです．それでも乱数の検定法で見モデリング，［３２‐３８］は運動論方程式の解法，［３９‐４２］は流抜けないような，疑似乱数の隠れた非ランダム性が計算結体方程式の解法，［４３］はプラズマ‐壁相互作用研究におけ

果に影響を与える可能性を排除したいのであれば，物理乱る動的モンテカルロ法，［４４］は中性子・光子輸送に関して数を利用することで，数値解の妥当性を判断する際に，そ書かれています．興味を持たれた方は，これらの文献を参のような影響を考えなくてよくなります．

照ください．
（菅野）

ところで，例題（１９）では，式（１７）における指数関数の引数に含まれる "が負の値で，計算の実行に問題は無かった

５．
２．
２ 輸送方程式求解への応用のですが，"が正の値を取る場合には，式
（１７）
のままでは

プラズマ輸送研究や統合モデリングにおいて，プラズマ

不具合（テスト粒子によっては，指数関数部分の数値計算の巨視的な時間発展や定常状態における分布を求めるためが破綻するなど）が発生することがあります．そのようなに輸送コードが広く用いられています．とりわけ，運転シ

不具合への対処の１つとして，式（１４）において "%
%&
$%
%&

ナリオ開発や実験解析においては定常状態を対象とするこを右辺に持っていき，ソース項 !%%&%&$%%&として処!"%とが多く，非定常輸送コードでプラズマの分布が変化しな理すると，数値解が落ち着くまで再帰的に計算を行う必要

くなるまで時間発展させて定常状態の分布を求めていまはありますが（つまり，ソース項における $%%&として，最す．輸送コードは放物型偏微分方程式である１次元の輸送

初は適当な近似解を与え，それにより得られた近似解をま方程式に基づいており，通常は有限差分法や有限要素法にたソース項の $%%&に代入して計算を繰り返すことで），う

よって時空間方向に離散化して数値的に解きます．温度分

まく克服できる場合があります［２５］．例として，以下のよ布を求めたい場合，圧力（内部エネルギー）を従属変数とうな方程式を考えます．

した輸送方程式を解きます．拡散係数（輸送係数とも言い

%## #
#!
"##
%%#

ます）がプラズマの内部状態に対して不変である場合，またはその応答が線形である場合は求解は容易です．しか（２１）

し，現実のプラズマの輸送は新古典輸送や乱流輸送によった だ し，%$(!
!
")と し て，境 界 条 件 は "%!&#"お よ びて支配され，それらは温度勾配などの熱力学的力に非線形

"%
%&
#%&
"&
&'%
#!"とします．式（２１）の解析解は， #%#$に依存しています．そのため，輸送方程式は非線形の偏微です．計算結果を図３に示します．ここで，#%&として最%分方程式となり，解析的に一般解を求めることは困難です

初に与えた近似解は，#%&%""です．シミュレーショ
%#!#し，数値的に解けたとしても数値的安定に求めることも困ンにおけるその他の設定は，図２と同じです．図３に示し

難になります．従来型解法の枠内で数値的安定に定常解をたように，式（１７）で !%%&#!かつ "###とすると計算に失求めるアルゴリズムが提案されてきましたが
［４５，
４６］，期
待した結果が得られないこともしばしばです．
本節では，輸送方程式を解いてプラズマの定常分布を求めるという問題を，定常状態を満たす温度と温度勾配の適切な組み合わせを求める大域的最適化問題へと転換し，その解法としてメタヒューリスティックなアルゴリズムである遺伝的アルゴリズム［４７，４８］を適用した例を紹介します．
ここまで読まれてきて，この話が一体どう乱数と関係するのか疑問に思われたかもしれませんが，遺伝的アルゴリズ
ムの中で乱数は重要な働きをします．以下では，まず遺伝的アルゴリズムの簡単な紹介をしたのちに，遺伝的アルゴ
リズムを用いた輸送方程式の求解法へと進んでいきます．
!遺伝的アルゴリズムとは

図３

遺伝的アルゴリズムは進化的アルゴリズムの中でも代表方程式（２１）の解析解と数値解．解析解を実線，式（１
７）
で
７）
で
F (x) ＝ 0 か つ "＝ #2 と し た 場 合 は 点 線，式（１F (x) ＝ − #2f (x)かつ "＝ 0 とした場合の数値解（得られた数値解の相対誤差 $＝ 1/100）を破線で示しています［２５］．

的なものであり，自然界における自然淘汰，遺伝，交叉，突然変異などの生命の進化過程を数値的アルゴリズムで模したものになります．遺伝的アルゴリズムでは必然的に進化
４１１

Journal of Plasma and Fusion Research Vol.96, No.8 August 2020

遺伝学の用語が多く使われますが，あくまで数値的な最適言います）の集団を生成します．扱う問題にも依りますが，

化アルゴリズムであり，生物学の文脈で使われる概念とは一般的に多くのテスト粒子数を必要とするモンテカルロ法

必ずしも厳密に一致しませんのでその点はご留意くださとは異なり，個体数は典型的に１００程度です．アルゴリズム

い．一口に遺伝的アルゴリズムと言っても実装は様々かとの流れは以下の通りです．
（１）パラメータ空間の中で個体思いますが，以下では PIKAIA［４９，５０］という Fortran 言語

をランダムに生成し，各個体の適応度を評価します．
（２）

における実装を基に，遺伝的アルゴリズムがどのようにし次に，適応度に応じて現在の集団からつがいを選択し交配て大域的最適解を求めるのかを見ていきましょう．

させ，次世代の集団を作ります．
（３）次世代の集団を現世最適解，すなわち解を求めるためには，当然ですが最適代と置き換え，（４）新しい集団の個体に対して適応度を評さの基準が必要になります．この品質の尺度を適応度と呼価します．集団の中で最適な（最大の適応度を持つ）個体

び，適応度を求めるための関数を適応度関数 #!*と言いま)が，設定した基準を満たすか設定した世代数に到達するます．大 域 的 最 適 化 問 題 を 解 く と は，パ ラ メ ー タ 空 間

で，（２）‐（４）の過程を繰り返します．最終世代の集団で最!*を最大化するパラメータ !を見つけるこ
!(.
"
!
"/で #
)
とを言います．なお，!は適当な規格化によって.
"
!
"/に収

適な個体の持つパラメータ &
!*
,
)!+&
-#
)
"が最適解とな
!""/
!(.
ります．
以上が遺伝的アルゴリズムのまるようにしています．

アルゴリズムになりますが，これだけだと乱数との関わりがはっきりしないかも

まず，これから生存競争を繰り広げる個体（表現型とも"#########################

乱

数

茶

話

５．
２．
１節で紹介したモンテカルロシミュレーションでです．

"#########################"#########################
は，拡散方程式を解く際に，対応する確率微分方程式に

試しに，拡散方程式（１３）で考えてみましょう．*)%*

従うテスト粒子の運動を利用しました．例えば，拡散方を−1/10 から＋1/10 の範囲の一様乱数として，次のよう

程式（１３）に対応した確率微分方程式（の計算スキーム）に式（９）に加えてみます．

として式（７）（数値計算する際には，式（９））を用いましたが，これ以外にも式（１３）に対応する確率微分方程式は

%
+
.)
%
%
.)
%
%
+
&%
&%
"!.#)
*
,
"&%
*
#.)
*
!+
*
"*)
*
,
$$+
(*
(*

あるのでしょうか？このことに関して，文献［１４］に興味

この式を用いて，%..
')
.)
%
*
*
/を計算すると，例えば，深い問いが載っています．それは，「いつ伊藤過程は拡%$!
$における分布は，図４のように図１と同じ結果(*#!

散過程となるか？」です．伊藤過程とは，次のように与になります．.%"とならない範囲で，*)%*-の最大値をえられる確率過程 ')%*を 指 し ま す．よ り 正 確な定義変えてみても，結果は変わりません．つまり，このようは，文献［１４］を参照ください．

な一様なノイズは，分布 #
%
"
.*
')
.)
%
)
#%..
*
*
/に影響し
ません．

%
%
%
%
'')
'%
'")
*
#()
*
"))
*
*

ここで紹介した文献［１４］の数学的結果は，微視的運動

ただし，()%%*お よ び ))*も 確 率 過 程 で，任 意の時刻

から巨視的スケールの分布が従う方程式へのモデリング%&!に対して確率１で以下を満たすとします．

にとって，重要な意味を持つように思えるのですが，い

%

%

!

!

かがでしょうか．
##()$$'$#' および # ))
*
*'$#'

（菅野）

文献［１４］による「問い」への答えですが，以下の拡散過
程 &)
%
*
%
&)
%
&)
%
%
'&)
'%
'")
*
#!)
*
*
",)
*
*
*

と ')
%
%
*との間に，()*に対する次の条件付き期待値につ
いて
!%/
()
%
'&)
%
%&.
*
#!)
*
*

かつ，))
%
*について
%
,)
'&)
%
))
%
))
%
'&)
*
*
*
*
*
*
#,)
図４

が，ほとんどすべての時刻 %，およびほとんどすべてのテスト粒子に対して成立する場合，&)
%
%
*と ')
*は，法則
の意味で一致する（つまり同じ分布を持つ）というのです．ここで，'&)%%*は %#!で位置 &から出発した ')*を意味し，!%は +; $%%‐加法族［１４］')$*,の生成する ,ピッチ角散乱によって時間発展する電子分布に対する一様なノイズの影響．実線は図１の時刻 t/ei ＝ 0.5 における分布，記号〇は一様 な ノ イ ズ を 加 え た 場 合 の 時 刻t/ei ＝ 0.5 における分布．一様乱数 *(t)は，Tausworthe法を利用して，−1/10 から＋1/10 の範囲で与えました．初期分布，各パラメータの設定値，およびテスト粒子の総数は，図１と同じです．

"#########################
４１２

Lecture Note

5. Transformation of Random Numbers and its Applications

S. Satake et al.

しれません．そのため，乱数との関連という観点から各段としましょう．座標値は有効数字４桁で取ったので，合わ階をもう少し詳細に見ていきたいと思います．

せ て８つ の 遺 伝 子 を 持 っ た 個 体 か ら 各 々 の 染 色 体

まず，個体とは何かについて説明します．個体とは遺伝

SP1: 12345678，SP2: 86427531 を作ります．このエンコード子の集合からなるものとして捉えられており，最適化した過程は自明かと思います．次に，一点交叉と呼ばれる演算

いパラメータの組の数だけ遺伝子の組があると思ってくだを行います．両親の染色体のつがいに対して乱数によってさい．例えば，２次元空間 &（ここで，%
）に
%
"
&'
"
&%(
"
!
")

決まった位置で染色体を切断し，付け替えます．この場合

おいて&
!
"
!
%
"
&'を求めるとします．念!
&!
&'に最も近い !$&ですと，一様整数乱数 "%(""))を生成し，例えば４が選ばのためですが，理想的な答えは&
!
"
!
!
&
!
&'なので自明です！

れたとすると，３番目と４番目の遺伝子座の間で組み替え

あくまで遺伝的アルゴリズムの概念 を わ か り や すく伝が 起 こ る た め，染 色 体 は そ れ ぞ れ SO1: 12327531，

え る た め に 設 け た 問 題 で す．適 応 度 関 数 の 設 定 方

SO2: 86445678 となります．PIKAIA では世代ごとに集団の法は色々あり ま す が，ユ ー ク リ ッ ド 距 離 の 二 乗 の逆数

数は増減しないため，２人の親から２人の子が生まれるこ

#

# !"

!'
$
%!!
!
&!!
!
&'"&
&') とすれば良いことがわかり&#(&とになります．つまり，これらは子の染色体となります．

ます．%， &は何桁かの数字（遺伝子）からなっており，今更に，ある変異率（通常は 0.005 程度です）で突然変異を起回の例では座標値になる訳ですが，その初期の座標値をこします．染色体を構成する遺伝子ごとに生成した一様乱"!
")の一様乱数から生成し，個体数分だけ繰り返します．
(

数 #%(
"
!
")が変異率を下回った場合に，その遺伝子を一

これが（１）で行っている作業です．オリジナルの PIKAIA

様整数乱数 "%("!
*)で置き換えます．これを一様な一点では線形合同法による乱数を用いていますが，筆者はそれ突然変異と呼びます．突然変異のさせ方は他にも色々あり

を Mersenne-Twister（メルセンヌ‐ツイスター）に換装しますが，詳細な説明は割愛します．ここでは SO1 に対してて使用しています．

は６番目の遺伝子が９に，SO2 に対しては２番目の遺伝子（２）の過程は遺伝的アルゴリズムの中でも最も重要なもが３になったとしましょう．これをデコードすると，子はのです．個体の集団の中から，１組のつがい，すなわち父それぞれ O1:&
，O2:&
!
"
!
!
"
!
!
"#$#
!
(*$"'
!
)$%%
!
&'()'となりまと母を選択します．父と母をランダムに選択するのではな

す．こうして生まれた子世代は一旦プールされ，親世代とく，適応度に応じて選ばれる確率を変える必要がありそう同数になるまで上記の過程が繰り返されます．なお，ランです．そのために，まずルーレット選択と呼ばれるサンプ

キング選択方式からわかるとおり，適応度の高い個体は１リング手法を用いると，どうなるか見てみましょう．そのつの世代で何度でも親に選ばれることがあります．

方法ではまず，全個体の適応度の和 !と，個体番号が若い

子世代に個体が出揃うと，（３）の世代交代を行います．

方から順に定義される適応度の累積分布を作成し準備して最も簡単で広く使われている方法は，集団を完全に子世代おきます．そして，一様乱数 #%("!)を生成すると，その!

の集団で置き換えてしまう，というものです．但し，エ乱数#は必ず離散的な累積分布のどこかの区間に対応する

リート選択という手法を使う場合は，親世代で最も適応度値を取り，その区間を示すインデックス値から特定の個体が高い個体は無条件で子世代へと引き渡されます．交叉や番号と一対一対応が付けられるために，乱数によって個体突然変異によってその世代での最適個体を壊してしまうのを選択できます．ルーレット盤は個体数で分割されておを避けることで，収束を悪くしないための工夫です．子世り，各個体の適応度の大きさによって区間の面積が決定さ

代からは１個体がランダムに排除されてしまいますが，仮れているイメージです（イメージしにく い 方 は［４８］のに選ばれた子世代の個体が親世代の最適個体よりも適応度

図２．
５をご覧ください）．これで問題は全くないように思いが高い場合は，親世代の最適個体は次世代に引き継がれまますが，実は適応度の大きさに応じて決めてしまうと，多せん．

くの問題が生じることが知られています．わかりやすい例

最後に，こうして作られた子世代の適応度を評価するのを挙げると，１つか２つの個体が極端に高い適応度を持っが（４）です．P1，P2，O1，O2 の適応度を評価すると，O2ておりその他がそうでない場合，このサンプリング方式でが最も大きな値を示していることがわかります．これははほぼ常に適応度の高い同じ個体が選ばれてしまうことにO2 が &!
"!
!
&!
&'に最も近い個体であることを示しており，なり，最適化過程にとって重要な集団の多様性が失われてこの例では，１世代を経てよりしまいます．それを避けるために，ランキング選択と呼ばたことになります．参考までに，このシミュレーションのれる方法が使われます．適応度の高い順に個体を並べるの結果を図５に示します．図５では各世代において最適個体

ですが，個体の選ばれやすさはランク（順位）の高さで決の持つ座標値が世代ごとに示されています．既に第１世代

まり，適応度の値そのものでは決まらない，という方式でから正解にかなり近く，１０世代程度で正解と見なして申しす．つまり個体間での選択確率の差が同じである，という

分ない精度に到達し，第２５世代で完全に正解に到達しま

わけです．ランクを横軸に取り，相対適応度を縦軸に取る

す．なお，適応度関数 $の定義から，&
!
"
!
&
!
&'における適!

と，右肩下がりの直線を引くことができますが，この勾配応度は無限大になりますが，数値的に工夫することで無限を制御することによって，淘汰圧の強さを調整することが大を避けるようにしています．ただし，遺伝的アルゴリズ

できます．

ムはメタヒューリスティックアルゴリズムであるため，正優秀な子どもが生まれこ の よ う に し て 選 ば れ た つ が い の 父 を解に必ず到達する保証はないことに注意してください．こ

P1: &，母を ,#+
%
"
&'
!
"
!
%
"
&'
!
"
!
!
"#$%
!
&'()'
!
)'%#
!
(&$"'
#&
&
#&

の例では &!
"!
!
&!
&'にたどり着きましたが，選択する乱数
４１３

Journal of Plasma and Fusion Research Vol.96, No.8 August 2020が算出する輸送係数に大きな影響を与えてしまいます．そこで微分に依らない解法である遺伝的アルゴリズムに着目し，大域的最適化手法を用いた定常輸送方程式の解法の開発に取り組みました［１８］．
以下では，簡単のため１粒子種のみを考えます．定常の輸送方程式は'# '$ &""%!#! ". '.!!))
+
&%+
*
' '
'
%
%

（２２）

と書けます．%は規格化小半径，' は体積，.は %微分，)は密度，$は拡散係数，&は温度，%は加熱パワーなどに対#応するソース項になります．)+&%+*は平衡計算で求められ

図５

各世代における最適個体の持つ遺伝子（x，y座標）．最初の１０世代でほぼ正解の(0.5, 0.5)に収斂し，２
５世代で完全に正解と一致します．

る メ ト リ ッ ク で あ り，こ こ で は 既 知 で あ る と し ま す．

シードによっては 0.4999 などで止まってしまうこともあり

分方程式となります．'")"%が既知である場合に，この方

)"
)."
$#$'
&"
&."
/(は一般に密度，温度やそれらの勾配などに非線形に依存しているため，式（２２）は非線形の偏微ます．しかし，それでも十分正解と見なせるでしょう．こ

程式を満たす &を求める，というのがここで考える問題でこまで簡単に遺伝的アルゴリズムの実装を見てきました

す．ここで，式（２２）の括弧内は熱流束 $ であり，が，数多の過程において一様乱数が用いられていることが

# '
$'
%(
$&
#!))
+
&%+
*
'
%

わかるかと思います．
!大域的最適化手法を用いた輸送方程式の求解

（２３）

を使って

遺伝的アルゴリズムは大域的最適化手法の一つとして位置付けられており，多峰性を持ったパラメータ空間におい

'
'.$(
!#! '
"%
'
'

て最適解（最大値とします）を見つけることを得意としています．多峰性の例として，*%,"!
"-における正規分布の重# #ね合わせからなる関数 +'
*(
*!*
&(+*,
###
!'
((#
(#"!()（２４）と書き換えられます．式（２４）を空間積分すると，を 考 え，そ の 最 大 値 を 見 つ け て み ま し ょ う．こ こ で，

%

'.$'
%(
%(
#$ '.%(%$#'

!
"
!
"
!
"
!
"
!
"
!
*
&""
!#"
*
&#!""
!
'
!
#
!
!#%
!
$
!
&
!
"%-と し ま
#,
,
""
#"

!

!
#を中心とする急峻な正規分布の頂す．最大値 0.9 は *#!

（２５）点に位置する一方，幅広い裾野を持った正規分布に極大値となります．ある磁気面を横切って外側に流れ出る熱流束あります．Newton 法などに代表される局所最適化手法では，その磁気面より内側に吸収されたパワーと等しいこと

は，極大値ではなく最大値を見つけられるかどうかはひとを表しています．式（２３）の温度 & を変化させることで式

えに初期値の選択に掛かってきますが，問題を解く前に最（２５）を満たすように$を変えることから，##'.を目標熱流

大値にたどり着く適切な初期値を知っていることは稀で束と呼びます．

す．局所勾配に依存する最適化手法では，*!0.282 に初期

では，どのように式（２５）を満たす &を見つければ良いの値を取ると極大値 0.3 へと向かってしまい真の最大値 0.9でしょうか．通常の手法ですと，空間方向にメッシュを作を見つけることはできません．言い換えると，ランダムに成し，有限差分法や有限要素法で空間離散化し，非線形項

初期値を選んだ場合，７割強の確率で最大値を求められなである $を線形化して & に対する連立一次方程式を作いことになります．大域的最適化手法は局所最適化手法とり，それを直接法ないし間接法で解き，得られた解（温度）比べ一般に計算コストが掛かる一方で，アルゴリズムにもを再び $に代入して連立一次方程式を解き直し，その作業依りますが局所最適化手法よりも広くパラメータ空間を俯を温度が収束するまで繰り返します［５１］．ここではそのよ瞰する特徴があります．極大値近傍に適応度の高い個体がうな一般的な方法を採らずに解&を求める方法を考えてみ!
#近傍に個集まっている状況下でも，突然変異などで *#!

ましょう．温度勾配を規格化した "
&#
'
%(
#
& で式
#
"& $'
!'

体が生じるとそちらの方が高い適応度を持つため，局所最（２３）を表記すると，

適解から抜け出して最大値を見つけ出せます．

#
$'
%(
$'
&"
#
"& (
& " # #.
"
#!))
+
&%+
*
"& '

遺伝的アルゴリズムは最大値を見つける際に微分値を必要としません．そのため，微分を高精度に評価することが（２６）となります．最後の等式は式（２５）から来ました．式（２６）で困難であったり，僅かな微分値の違いが結果を大きく左右してしまったりするような問題であっても，安定に解を見は $の&"#"& の依存性を陽に書きました．温度と温度勾配"つけることができます．そのような問題の例としては，分には物理的に現実的な定義域がありますから，そのパラ布の勾配に強い依存性を持つ輸送モデルを用いる場合が挙メータ領域内で各%において式（２６）を満たす'&"#"& (の組"

げられます．分布の勾配という微分の精度が，輸送モデル

み合わせを見つけることができれば，それが式（２２）の解と４１４

Lecture Note

5. Transformation of Random Numbers and its Applications#

S. Satake et al.
#

$'!$)
"!%,#
'
"， '
"
"#!
##!
$'
"

なります．この手法だと，直接温度の微分値である ""!$を探しているため，微分方程式の解を求めるのに微分が不（２８）を評価し，最終的に遺伝的アルゴリズムによって要であることがわかると思います．
解の組み合わせを求める手法ですが，極端なことを言え

!"
'
'
'
($
*%
#*
+
"!
#&

ば，そのパラメータ領域内を極度に細かく分割して全ての（２９）組み合わせに対して当てはめを行えば，どこかで解は見つを 最 大 化 す る 個 体 を 見 つ け ま す．そ の 個 体 の 持 つかるでしょう．現実的にはそこまで細かな分割はできませ$!
"!$ &が最終的な式（２２）の解となる訳です．'
"
%
#の意味は

んので，計算精度と掛かる計算時間を秤にかけて妥当な刻わかりにくいかもしれませんが， '#が小さいとはすなわ

み幅を決めざるを得ません．非線形性の強い輸送モデルをち，ある温度 $)を入力として用いて式（２７）を計算した時に用いる場合，僅かな勾配の差が大きな熱流束の差を生んで得られる$'が入力の$)とがほぼ同じである，ということをしまうため，先程まで極端に細かくないにしてもやはり十意味しています．そのような $ は式（２７）を 満 た し て い

分細かな分割が必要になります．規格化温度勾配の定義域

る，ということがわかるかと思います．

を*!
!
"!+として，刻み幅を 0.01 とかなり粗く取ったとして!計算例

も，１０００メッシュ必要になります．温度に対しても同じく

上記のアルゴリズムで温度分布を解いた計算例を一つ示

6１０００メッシュ用意すると，組み合わせは 10 にも上ります．

します．計算は，上記のアルゴリズムに基づく定常輸送

これを，空間分割数分繰り返さなくてはいけないため，膨

コード GOTRESS で行いました．加熱分布や平衡は JT-60U

大な試行が必要になります．さらに，この組み合わせの中のとある放電実験のデータから持ってきました．輸送モデには解として全くあり得そうに無い組み合わせも存在してルには CDBM［５２］を用い，径方向５０点で１粒子種（電子にいるため，無駄の多い手法であることがわかります．その相当します）の温度分布の計算を PC クラスタで実行した

ような無駄を省くために，この解の組み合わせを見つける

ところ，約 2.77 秒で計算を終えました．平均すると，各径アルゴリズムとして大域的最適化手法である遺伝的アルゴ

方向点当たり２１２世代で収束解が得られていることになりリズムを用いるのです．遺伝的アルゴリズムでは，最初こます．なお，世代ごとの個体数は１００であり，CDBM の評価

そモンテカルロ的にランダムに値を生成しますが，そこかを行った回数の総数は１
０６万６千回になります．図６が示らはより解に近い（適応度の高い）組み合わせを持った個

すとおり，各径方向位置において目標熱流束に対応する加体がより生き残るため，急速に最適解へと収斂していきま

熱パワー累積分布 "（実線）と輸送熱流束に対応するパす．

ワー累積分布 %,#（丸印）が一致していることがわかりま,

熱流束を求める時は， 正解 である目標熱流束 ""%がす．これは，式（２６）を満たしていることを示しています．

事前にわかっていたため，# を目標熱流束に合わせることその時の規格化温度勾配分布は図７に，温度分布は図８にができたのですが，正解がわかっていない温度 $はどのよ

示されています．温度分布や規格化温度勾配分布が共に滑うに求めたら良いでしょうか．一つの考え方は，""!$ の値らかに得られていることがわかるかと思います．この時のが定まったのだったら，それは温度に対する常微分方程式

規格化温度勾配分布は温度勾配を何らかの微分スキームでになるため，境界の温度さえわかれば後は全て積分で求め微分したものではなく，温度分布と独立に得られたものでられる，というものです．しかし，この手法は上手くいかあることを強調しておきます．また，図８には拡散係数分ないことがあることがわかっており［１８］，別の手法を考え布も示されています．GOTRESS 内部では新古典熱拡散係,る必要があります．式（２３）で % ##"であることに留意し数も計算しているため，CDBM で評価された乱流熱拡散係て外側境界から積分してみましょう．すると，数との和が示されています．なお，磁気軸においては規格

$&

$%
$&
$%&
#$%
"#

$

&$
"
,
#
('
#%
$!
"
!$ &%
"
)
$$)
(

化温度勾配，つまり輸送熱流束は０であることが自明であ（２７）という関係が得られます．$%$%&は計算領域の外側境界$#$%での固定された温度であり，Dirichlet 条件に相当するものです．式（２７）は非線形の積分方程式であり，被積分
関数中の #に陰に含まれる$（$)とします）と左辺の$（$'）
が同一となるような $を見つけなければなりません．本来
このような解を見つけることは難しいのですが，大域的最適化手法によって比較的容易に見つけ出すことができます．ところで，式（２７）は ""!$ にも陰に依存しているため，"!$ も適切に定まらなければ収束解にはたどり着きませ"ん．つまり，式（２６）も同時に満たさなくてはいけません．
（２６）
および（２７）を計算
各個体の持つ %$!
"!$ &に対して式"し，得られた # と $ に対して

図６

４１５

目標熱流束に対応する加熱パワー累積分布 P と輸送熱流束に対応するパワー累積分布 V´Q ．

Journal of Plasma and Fusion Research Vol.96, No.8 August 2020

参考文献

図７

図８

［１］津田孝夫：モンテカルロ法とシミュレーション（三訂版）（培風館，１９９５）．
［２］伏見正則，逆瀬川浩孝（監訳）：モンテカルロ法ハンド
ブック（朝倉書店，２
０１
４）．
［３］D.R. Nicholson, Introduction to Plasma Theory (Wiley, 1983);
小笠原正忠，加藤鞆一
［訳］：プラズマ物理の基礎（丸
善，１
９
８
６）
．
［４］D.G. Swanson, Plasma Kinetic Theory (Chapman and Hall/
CRC, 2008).
［５］プラズマ・核融合学会［編］：プラズマシミュレーショ
ン（京都大学学術出版会，２０
１
８）．
［６］渡邉智彦，洲鎌英雄：プラズマ・核融合学会誌 81, 534(2005).
［７］渡邉智彦：プラズマ・核融合学会誌 81, 686 (2005).
［８］前山伸也：プラズマ・核融合学会誌 91, 589 (2015).
［９］S. Maeyama et al., Phys. Rev. Lett. 114, 255002 (2015).
［１０］M. Nunami et al., Phys. Plasmas 27, 052501 (2020).
［１１］K. Fujita et al., Plasma Fusion Res. 14, 3403102 (2019).
［１２］J.P. Freidberg, Ideal Magnetohydrodynamics (Plenum, 1987).
［１３］M. Honda, Comput. Phys. Commun. 231, 94 (2018).
［１４］B. !ksendal, Stochastic Differential Equations: An Introduction with Applications (Springer, 2003); 谷口説男
［訳］
：確
率微分方程式―入門から応用まで（丸善，１
９
９
９）
．
［１
５］A. Friedman, Stochastic Differential Equations and Applications (Dover, 2004).
［１６］保江邦夫：確率論（数理物理学方法序説４）
（日本評論
社，２
０
０
１）
．
［１
７］保江邦夫：物理数学における微分方程式（数理物理学方法序説別巻）（日本評論社，２
０
０
２）
．
［１
８］高橋陽一郎［編］：伊藤清の数学（日本評論社，２
０
１
１）
．
［１
９］P. Helander and D.J. Sigmar, Collisional Transport in Magnetized Plasmas, Chapter 3 (Cambridge University Press,
2002).
［２０］P.E. Kloeden and E. Platen, Numerical Solution of Stochastic Differential Equations (Springer, 1999).
［２１］T. Misawa and H. Itakura, Phys. Rev. E 51, 254 (1995).
［２２］J.R. Klauder and W.P. Petersen, SIAM J. Numer. Anal. 22,1153 (1985).
［２３］A.H. Boozer and G. Kuo-Petravic, Phys. Fluids 24, 851(1981).
９８３）．
［２４］坂元慶行 他：情報量統計学，第４章（共立出版，１［２５］R. Kanno et al., Plasma Fusion Res. 6, 2403066 (2011).
［２６］T. Takizuka and H. Abe, J. Comput. Phys. 25, 205 (1977).
［２７］K. Nanbu and S. Yonemura, J. Comput. Phys. 145, 639(1998).
［２８］C. Wang et al., J. Comput. Phys. 227, 4308 (2008).
［２９］Y. Masaoka and S. Murakami, Plasma Fusion Res. 8,2403106 (2013).
［３０］G. Zhang and D. del-Castillo-Negrete, Phys. Plasmas 24,
092511 (2017).
［３１］S. Satake et al., Comp. Phys. Comm. 255, 107249 (2020).
［３２］X.Q. Xu and M.N. Rosenbluth, Phys. Fluids B 3, 627 (1991).
［３３］W.X. Wang et al., Plasma Phys. Control. Fusion 41, 1091
(1999).
［３４］S. Brunner et al., Phys. Plasmas 6, 4504 (1999).
［３５］S. Murakami et al., Nucl. Fusion 40, 693 (2000).
［３６］S. Murakami et al., Nucl. Fusion 46, S425 (2006).
［３７］S. Matsuoka et al., Phys. Plasmas 22, 072511 (2015).

規格化温度勾配 1/LT の定常解．

温度 T の定常解とその時の拡散係数 !．拡散係数は CDBM
で評価された乱流熱拡散係数に新古典熱拡散係数が重畳されたものです．磁気軸での拡散係数は０としています．

るため，磁気軸では計算しておらず，磁気軸での温度は温度勾配零の条件から算出しています．そのため磁気軸での拡散係数を算出する必要はなく，０と置いています．
本節では，一見およそ繋がりがなさそうな遺伝的アルゴ
リズムと偏微分方程式の求解が，解法の工夫によって結びつけられ，輸送方程式を解くという実際の問題へと適用できることをご 紹 介い た し ま し た．輸 送 コ ー ド と してのGOTRESS は複数粒子種を扱えたり並列計算が行えたりす
るなど多くの特徴を備えています．詳細は［５３］を参考にしてください．
（本多）

謝

辞

５．
２．
１節 に お け るシ ミ ュ レ ー シ ョ ン で 利 用 し た Tausworthe 法に基づく疑似乱数発生プログラムを提供いただ
いた，高丸尚教氏（中部大学 工学部 教授）に感謝いたします．また，５．
２．
１節に対して，松山顕之氏（量子科学技術
研究開発機構 六ヶ所核融合研究所 主幹研究員）より，貴重なコメントをいただきました．ここに感謝申し上げます．

４１６

Lecture Note

5. Transformation of Random Numbers and its Applications

S. Satake et al.

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

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

しんすけ

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

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

たけ!!!!!!!!!!!!!!!!!!!!!!!!!!

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

ほん

だ

みつる

充

量子科学技術研究開発機構 那珂核融合研究所 先進プラズマ研究部 上席研究員．主にトカマク中の輸送現象や輸送シミュレー
ションを研究対象としていますが，機械学
習やベイズ推定，最適化問題にも興味を持っていて研究への応用を進めています．家族でカラオケボックスに行くことが多かったのですが，昨今の情勢でなかなかままなりま せん．９
５点以上を獲れる曲を増やすのが目標です．

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

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

本 多

かん

の

りゅう

た

ろう

菅 野 龍 太 郎
自然科学研究機構 核融合科学研究所 ヘリ
カル研究部 核融合理論シミュレーション
研究系 准教授．乱数を用いた確率論的な計算手法であるモンテカルロ法全般に興味があります．最近の研究では，プラズマの衝突輸送現象に対するモンテカルロ法に基づいたドリフト運動論シミュレー
ションを行っています．

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

さ［５２］M. Honda and A. Fukuyama, Nucl. Fusion 46, 580 (2006).
［５３］M. Honda and E. Narita, Phys. Plasmas 26, 102307 (2019).

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

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

An Introductory Analysis with Applications to Biology, Control, and Artificial Intelligence, second ed. (MIT Press, Cambridge, 1992).
［４８］C.R. Reeves, J.E. Rowe, Genetic Algorithms: Principles andPerspectives (Kluwer Academic Publishers, New York,
2002).
［４９］P. Charbonneau, Astrophys. J. Suppl. Ser. 101, 309 (1995).
［５０］http://www.hao.ucar.edu/modeling/pikaia/pikaia.php［５１］J. Candy et al., Phys. Plasmas 16, 060704 (2009).

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

［３８］H. Yamaguchi and S. Murakami, Nucl. Fusion 56, 026003(2016).
［３９］Y. Feng et al., J. Nucl. Mater. 266-269, 812 (1999).
［４０］A.M. Runov et al., Phys. Plasmas 8, 916 (2001).
［４１］M. Kobayashi et al., Contrib. Plasma Phys. 44, 25 (2004).
［４２］R. Tatsumi et al., Plasma Fusion Res. 15, 1403003 (2020).
［４３］プラズマ・核融合学会［編］：プラズマシミュレーショ
ン，第７．
４．
１節（京都大学学術出版会，２
０
１
８）
．
［４
４］桜井 淳：日本原子力学会和文論文誌 2, 555 (2003).
［４５］S.C. Jardin et al., J. Comput. Phys. 227, 8769 (2008).
［４６］G. V. Pereverzev and G. Corrigan, Comput. Phys. Commun. 179, 579 (2008).
［４７］J.H. Holland, Adaptation in Natural and Artificial Systems:!!!!!!!!!!!!!!!!!!!!!!!!!!

４１７

J. Plasma Fusion Res. Vol.96, No.8 (2020)４
１８‐４１９講座 シミュレーションのための乱数入門

６．おわりに
6. Summary
宇佐見俊介
USAMI Shunsuke
自然科学研究機構 核融合科学研究所（原稿受付：２
０
２
０年５月２
２日）

本講座では，コンピュータシミュレーションにおける乱ました．その一方で，「これを通過したら，その乱数の品数の実践的な入門書として，様々な乱数発生方法の原理か質は保証される」という万能検定は存在せず，多くの検定ら始めて，その高速化・並列化手法，乱数の品質・検定法，をクリアすることで，「良い乱数らしい」という確度を上

そして，乱数の変換・利用例まで紹介しました．以下では，

げるしかない，ということが強調されました．また，乱数これらの内容を振り返って，本講座のまとめとさせていたを並列化する場合，各並列プロセスに配分された乱数列同だきます．

士の独立性が重要な品質ファクターとなることを論じまし

第２章では，コンピュータにおいて乱数を生成する原理

た．その独立性を議論する文献がないことから，著者が独を，決定論的な演算によって作られる疑似乱数とランダム自に行った独立性に関するテストについて紹介しました．

な自然現象を利用する物理乱数の２つに分けて説明しまし

第５章では，乱数の変換例と利用例を示しました．まずた．２．
１節では，疑似乱数の発生法として，線形合同法，

５．
１節では，一様乱数から任意の確率密度分布に従う乱数

Tausworthe（トーズワース）法，Mersenne-Twister（メルを作る方法として，逆変換法（直接法）と棄却法を紹介しセンヌ・ツイスター）法（MT 法）を取り上げました．特ました．逆変換法では，作りたい確率密度分布の累積分布に，MT 法については，サンプルプログラムを例示して実関数を利用しました．棄却法では，２組の一様乱数を用意際の使い方を解説しました．一方２．
２節では，物理乱数発生

し，一方は，他方の乱数要素の採用・棄却を判定するため器として，定電圧ダイオードの出力信号に含まれるゆらぎに用いました．その後，例として Box-Muller 法を取り上を利用したものと，レーザー光の偏光における量子ゆらぎ

げ，２次元空間で正規分布を生成する具体的な方法を解説を利用したもの（量子乱数）をピックアップし，それぞれしました．次に５．
２節では，乱数をシミュレーションにおいの原理を紹介しました．また，それらの物理乱数が利用でて，どのように活用しているのかについて紹介しましきるサービスについても紹介しました．

た．５．
２．
１節では，拡散方程式の解法と Dirichlet-Poisson

第３章では，並列プログラムで乱数を使用することを想混合問題を例として，乱数を用いた基礎的なシミュレー

定して，乱数発生の高速化チューニング方法を解説しまし

ション手法がどのようなものかを見ました．ここでは，乱た．MT 法による乱数生成の並列化方法として Dynamic

数は，Coulomb 衝突や熱拡散といったランダムな物理過程

Creator，および Jump Method を挙げ，特に後者についてを，確率過程論に基づいて具現化するために用いられていは，数値計算ライブラリ KMATH̲RANDOM をベースにました．また，５．
２．
２節では，プラズマ・核融合分野のみなして，その使い方を説明しました．また，乱数利用の効率

らず，様々な科学技術分野で注目を集めている「遺伝的ア化のため，「水瓶方式」というアイデアを示しました．使

ルゴリズム」を取り上げました．その適用例として，非線える物理乱数発生器が１個の場合においても乱数を並列プ形な輸送方程式を解く問題を大域的最適化問題へと転換

ログラムで効率よく利用する方法として，「源泉かけ流しし，遺伝的アルゴリズムを用いて解を求める方法を紹介し方式」というアイデアも説明しました．これら「水瓶方ました．こちらでは，乱数は，偏微分方程式を解くための式」，および「源泉かけ流し方式」は，本講座の著者である

アルゴリズムにおけるツールという役割を果たしていまし佐竹真介氏が考案したオリジナル手法です．

た．

第４章では，乱数らしさについて，周期性，ビット，多この第６章は，本講座のすべての原稿が出揃った段階で次元均等分布性（あるいは結晶構造）といった観点から詳書いていますので，いわば本講座の紹介一覧，あるいは目しく論ずるとともに，乱数が持つべき統計的性質を利用し次のような役割を果たしているとも言えます．読者の皆様ての乱数の検定方法を解説しました．具体的には，上昇連，が，本講座を振り返って読み返したい箇所，または自身の下降連，近接値の出現率，平均値，二乗平均値が挙げられ研究に必要な箇所にたどり着きやすいことを目標に，まとNational Institute for Fusion Science, Toki, GIFU 509-5292, Japanauthor’s e-mail: usami.shunsuke@nifs.ac.jp

４１８!2020 The Japan Society of Plasma
Science and Nuclear Fusion ResearchLecture Note

6. Summary

S. Usami

めました．最後に，第１章で述べたことの繰り返しになりの役割を果たし，知的好奇心を刺激するものとなっている

ますが，乱数を用いようとする，あるいは，乱数自体に興ことを願っております．

味がある研究者，学生などにとって，本講座が入門として!!!!!!!!!!!!!!!!!!!!

う

さ

み

しゅん すけ宇 佐 見 俊 介
自然科学研究機構 核融合科学研究所 ヘリ
カル研究部基礎物理シミュレーション研究
系 准教授．
粒子シミュレーションを用いて，プラズマの様々な複雑な現象を調べていて，現在の主な研究テーマは
磁気リコネクションです．SF が好きなためか，現実世界はコンピュータ上のシミュレーションでは？と密かに妄想しています．今回の講座をとりまとめている際に，世界がシミュ
レーションなら，量子ゆらぎに基づく乱数も実は疑似乱数
で，何らかの規則性があったりしないか？とさらに妄想を膨らませています．

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

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

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

４１９