数学アルゴリズム研究ノート (3)

3.衝突エントロピー検定による擬似乱数の評価

2003年8月 田沼英樹


3-1.衝突検定について


 前回まで考察した擬似乱数列について,周期に関する性質はすでに理論的に得られているわけですが,それが実際に擬似乱数として使い物になるかどうかを評価するためには,理論的または統計的に検定を行う必要があります.

 理論的検定としては,線形合同法のスペクトル検定やM系列乱数のランクによる検定などが有名ですが,一般に理論的検定は特定の乱数発生アルゴリズムの性質に依存していて,また大抵は全周期にわたる性質のみが得られるので,別種の擬似乱数を比較する場合や周期の一部のみの性質を調べる場合には適用できません.

 統計的検定は,一般の乱数列について特定の統計的性質に注目して検定を行うものですが,その中の一つに衝突検定というものがあります.

 正方形や立方体,あるいはさらに高次元の超立方体の領域に擬似乱数を使ってランダムに点をプロットしていくことを考えます.このとき,もし点の座標を定める擬似乱数に偏りがあれば,点は領域の全体に均等に分布せずに,不自然に同じ箇所に集まることになります.そこで,領域を同じ大きさの個のセルに分割し,毎回点をプロットする際に,すでに同じセルに別の点があった場合,衝突が一回発生したとしてその回数をカウントしていきます.

 点を回プロットした時点で個のセルが埋まっている確率をとすると


という漸化式が成り立ち,衝突回数が回である確率はとして求めることができます.このように衝突回数の確率分布が具体的に求まることから,ある決まった回数だけ試行した結果の衝突回数を評価することにより,擬似乱数の検定を行うことができます.これを衝突検定といいます.


 通常,試行回数はセルの数に対して十分に小さくしますが,試行回数を増やしていくと衝突回数の期待値が大きくなることから,正確な確率分布を求めるために必要な項数が多くなり,数値計算の誤差や計算時間あるいはメモリ使用量についての配慮が必要となってきます.そこで,試行回数が多いときに衝突回数を評価するために,衝突回数の期待値と分散を求めます.

 まず,回の試行の後埋まっているセルの個数(確率変数)をとします.このとき,衝突回数はとなります.ここで,の特性関数を


とおくと,上の漸化式より


という関係式が成り立ちます.これから



が導かれ,より衝突回数の期待値が求められます.

 衝突回数の分散については,漸化式


よりとして具体的に求めることができます.


 ここで,衝突回数は点がプロットされる順番によらないことから,試行回数を多くしたときに擬似乱数の短期間の偏りに関する情報が失われてしまう可能性があります.例えば,セルがまだ数個しか埋まっていない状態で衝突が十回続くのはかなり異常であるといえますが,セルが半分近く埋まっている状態で衝突が十回続くのは通常ありえることです.このどちらの場合でも,衝突判定では同じ十回としてカウントされてしまい,異常さが埋もれて検出できなくなってしまいます.このように,プロットの回数を増やしていくと確率的なノイズが増えてしまい,検定の精度がかえって落ちてしまうという問題があります.

 そこで,このような問題点に配慮した,さらに高精度な検定法についてこれから考察していきます.



3-2.衝突エントロピー検定


 前記と同様にセルの総数をとし,回の試行の後埋まっているセルの個数(確率変数)をとします.回目の試行で衝突が起きると,それ以外はとなる確率変数をとすると,


が成り立ちます.ここで,回目の試行で衝突が起きた/起きなかったという事象の自己情報量は


となります.この自己情報量はどれだけ珍しい事象が起きたかを表す尺度で,事象の起きる確率が低いほど,大きな正の値をとります.自己情報量の単位は,対数の底をとした場合はビット(bit),底をとした場合はナット(nat)とよばれるものになります.

 セルがまだ半分まで埋まっていないとき,衝突が起きる確率が起きない確率より低いことから,衝突が起きた場合の方が起きなかった場合に比べてより大きな自己情報量を持ちます.また,1回目の試行では明らかに衝突は起きないので,新たな情報がないことからとします.

 自己情報量を確率変数とみなし,その期待値をとったもの


は平均情報量または情報エントロピーとよばれるものになります.自己情報量と平均情報量の差は期待値がの確率変数になりますが,式をまとめると


となります.この値はいわば,理想的なランダムからの確率的なずれをエントロピーで表した量となりますが,実際,埋まっているセルの数が小さいときに衝突が起きると大きな値になり,セルが埋まっていくにつれて振れ幅が小さくなっていき,になると値はとなります.そして,となる範囲では衝突が起きたときに正,起きなかったときに負の値をとることから,衝突回数の大小を情報の重要さに応じて重み付きで表したものと見なすこともできます.

 そこで,回の試行について総和をとった値


を衝突エントロピーとよび、これを用いて擬似乱数の評価を行います.


 の和を構成する確率変数が互いに独立でないことから,このままではの確率分布について詳しく解析することができません.そこで,を期待値に置き換えてを近似したものを


とします.ここで,が互いに独立であるとすると,級数の各項も独立となるので,その分散は


となります.そして,期待値の一般式


を代入すると,


が導かれます.

 が大きくなると,の確率分布を実際に求めるのは困難になりますが,級数の各項が互いに独立で確率分布の変化があまり大きくないことから,部分的には中心極限定理の条件を満たしていると見なすことができます.したがって,が十分大きいとき,は正規分布に近い確率分布をもつであろうと予想されます.

 そして,による近似が十分に良いのであれば,おそらくも正規分布に近い分布をもつであろうと予想されるので,その標本値を分散の近似値を用いて評価することが可能になります.これを,衝突エントロピー検定とよぶことにします.



3-3.乱数検定の実例


 上記の理論に基づき,これから実際の擬似乱数列について衝突エントロピー検定を実施してみたいと思います.

 まず,多次元分布の偏りを調べるために,点をプロットする空間を次元の超立方体とします.各次元の座標にはビットずつを割り当て,セルの数をとします.各試行ごとに,乱数列から値を個ずつ取り出した組の,それぞれ通常は上位ビットずつを用いることによって,点が入るセルを決定します.

 回目の試行の結果,埋まっているセルの総数をとします.衝突検定では衝突回数をそのまま直接評価しますが,衝突エントロピー検定では,衝突エントロピーからどれだけずれているかによって,擬似乱数の評価を行います.試行回数となる範囲で任意なので,乱数の性質が視覚的にわかりやすいように範囲内でのの変化をグラフで表示するとともに,の最大値と最小値、および最終値を記録して評価します.

 さらに,値で評価するために,標準偏差をで近似して比を求めます.そして,一応参考のために衝突回数の期待値からの偏差を求め,同様に標準偏差で割って値を計算します.


 では,これからいくつかの擬似乱数の実例について,実際に検定結果を挙げていくことにします.



● Numerical Recipes ran2


 Numerical Recipes [1] で推奨されている ran2 とよばれる擬似乱数です.2つの周期の異なる線形合同法と切り混ぜ機構を組み合わせた擬似乱数で,約の周期を持ちます.



image001.png



  


 グラフの横軸は,縦軸はを表します.グラフの右端ではとなり,セルが半数埋まった状態に相当します.

 緑色の線はの変化を表していて,縦軸の範囲はの最大値と最小値で正規化されています.

 赤い線はを表していて,緑色の線がこの線より上にある場合は理想的なランダムより衝突が多く,下にある場合は衝突が少ないことを示しています.

 衝突エントロピーはビット単位で表され,青い線は(bit) およびの範囲を示しています.


 この ran2 の場合,序盤に若干衝突が多く,が一旦まで上昇しますが,その後は大体の範囲に収まっていることから,この検定結果では特に異常はありません.



● 線形合同法 


    image002.png


  


 これは典型的な線形合同法による擬似乱数ですが,明らかに異常な結果が出ています.おそらく,疎な結晶構造が出現しているために,衝突が非常に少なくなっているものと思われます.



● 二次合同法 


    image003.png


  


 これは初回に扱った二次合同法ですが,この結果を見る限り,序盤の変動が多少激しいことを除けば一見良い乱数に見えます.ところが,別の初期値を使ってグラフを描いてみると,今度は次のように異常な挙動を示します.


    image004.png


  


 このように,擬似乱数の短期的性質は初期値による変動が大きいので,乱数発生アルゴリズムの評価をする場合には注意が必要です.とはいえ,初期値によらず常に異常を示す線形合同法に比べれば,幾分ましであるとも言えます.



● シフトレジスタ列の二次合同式による延長


 C言語式 x = (0x80000057 &−(x&1)) ^ (x>>1) で生成されるシフトレジスタ列をとし,二次合同式で延長した擬似乱数列です.の周期を持ちます.


    image005.png


  


 種となるシフトレジスタ列はそのままでは分布に偏りが大きく,同様の検定を実施しようとしてもセルが半数まで埋まらないので結果を出すことができません.ところが,これを簡単な二次式で延長することにより,比較的良好な検定結果が得られます.

 この延長した列の擬似乱数としての性質はまだ未知数ですが,おそらくを定数とみなした単なる二次合同法より悪くなることはないと予想されます.また,シフトレジスタ列の生成多項式をさらに項数が多いものに変更すれば,より乱数として優れたものになるのではないかと思われます.

 なお,延長の式が前回扱ったではなくという形になっていますが,実はこの方がプログラムが簡単になり,また擬似乱数の性質もが一つずれるだけで,本質的には変わりません.



● Mersenne Twister MT19937


 いよいよ大御所の Mersenne Twister [2] です.の周期を持ち,次元にわたって均等に分布することが理論的に証明され,しかも高速に生成できるという,おそらく現時点で最強ではないかと思われる擬似乱数です.


    image006.png


  


 非常に良い結果が出るかと期待されたのですが,検定結果は意外にも偏りを示しています.特にセルが3分の1ほど埋まったあたりで急に衝突が多くなり,衝突エントロピーでは近くまで上昇しています.たまたま検定の条件がこの擬似乱数と相性が悪かった可能性もありますが,理論的な性質の良さはあくまでも周期全体にわたるものであって,それが必ずしも短期的な性質には反映されないのかもしれません.

 この検定結果はデフォルトの初期値,あるいは init_genrand(5489) とした場合のものですが,別の初期値として試しに init_genrand(0) としたものを検定してみると,次のような結果になります.


    image007.png


  


 今度は前に比べて格段に良い結果が出ています.この結果だけで結論を出すのは危険ですが,現在のデフォルトの初期値には多少問題があるかもしれません.



● Mersenne Twister MT19937 の二次合同式による延長


 Mersenne Twister をとし,二次合同式で延長した擬似乱数列です.の周期を持ち,次元にわたってほぼ均等に分布することが理論的に証明されます.

 なお,Mersenne Twisterはデフォルトの初期値で生成し,延長乱数列はという初期条件で検定を行いました.


    image008.png


  


 まだ若干衝突が多めですが,元の Mersenne Twister にあったような偏りは解消されています.

 次に,Mersenne Twister を init_genrand(0) で初期化した場合の延長乱数列について検定を行います.


    image009.png


  


 途中若干衝突が多くなりますが,最終的には元の Mersenne Twister と同様,良好な結果に落ち着きます.




参考文献&リンク


[1] Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P., Numerical Recipes in C: The Art of Scientific Computing (Cambridge University Press, 1988).

  丹慶勝市・奥村晴彦・佐藤俊郎・小林誠 訳 『Numerical Recipes in C [日本語版] C言語による数値計算のレシピ』 (技術評論社,1993).

[2] Mersenne Twister: A random number generator



次に進む


戻る