数学アルゴリズム研究ノート『衝突エントロピー検定による擬似乱数の評価』

数学アルゴリズム研究ノート（３）

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

2003年８月　田沼英樹

３－１．衝突検定について

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

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

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

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

　点を $r$ 回プロットした時点で $n$ 個のセルが埋まっている確率を $p_{m}(n, r)$ とすると

$\displaystyle p_{m}(n, r+1) = \frac{n}{m}p_{m}(n, r)+\frac{m-(n-1)}{m}p_{m}(n-1, r)$

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

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

　まず， $r$ 回の試行の後埋まっているセルの個数（確率変数）を $F_{r}$ とします．このとき，衝突回数は $r-F_{r}$ となります．ここで， $F_{r}$ の特性関数を

$\displaystyle f_{r}(x)=\sum_{n}p_{m}(n, r)\cdot x^{n}$

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

$\displaystyle f_{r+1}(x)=xf_{r}(x)+\frac{x-x^{2}}{m}f_{r}^{\prime}(x)$

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

$\displaystyle f_{r+1}^{\prime}(1)=1+\frac{m-1}{m}f_{r}^{\prime}(1), f_{r}^{\prime}(1)=1$

$\Rightarrow E(F_{r})=f_{r}^{\prime}(1)=m-m\left(\frac{m-1}{m}\right)^{r}$

が導かれ， $E(r-F_{r})=r-E(F_{r})$ より衝突回数の期待値が求められます．

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

$\displaystyle f_{r+1}^{\prime\prime}(1)=\frac{2m-2}{m}f_{r}^{\prime}(1)+\frac{m-2}{m}f_{r}^{\prime\prime}(1), f_{1}^{\prime\prime}(1)=0$

より $V(r-F_{r})=V(F_{r})=f_{r}^{\prime\prime}(1)-f_{r}^{\prime}(1)-{f_{r}^{\prime}(1)}^{2}$ として具体的に求めることができます．

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

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

３－２．衝突エントロピー検定

　前記と同様にセルの総数を $m$ とし， $r$ 回の試行の後埋まっているセルの個数（確率変数）を $F_{r}$ とします． $r$ 回目の試行で衝突が起きると $1$ ，それ以外は $0$ となる確率変数を $C_{r}$ とすると，

$\displaystyle F_{r}=F_{r-1}+1-C_{r}, P(C_{r}=1)=\frac{F_{r-1}}{m}, P(C_{r}=0)=\frac{m-F_{r-1}}{m}$

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

$\displaystyle I_{r} = -C_{r}\log P(C_{r}=1)-(1-C_{r})\log P(C_{r}=0) = -C_{r}\log\frac{F_{r-1}}{m}-(1-C_{r})\log\frac{m-F_{r-1}}{m}$

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

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

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

$\displaystyle E(I_{r}) = -\frac{F_{r-1}}{m}\log\frac{F_{r-1}}{m}-\frac{m-F_{r-1}}{m}\log\frac{m-F_{r-1}}{m}$

は平均情報量または情報エントロピーとよばれるものになります．自己情報量と平均情報量の差 $I_{r}-E(I_{r})$ は期待値が $0$ の確率変数になりますが，式をまとめると

$\displaystyle I_{r}-E(I_{r}) = \left(C_{r}-\frac{F_{r-1}}{m}\right)\log\frac{m-F_{r-1}}{F_{r-1}}$

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

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

$\displaystyle H_{n} = \sum_{r=2}^{n}\left(I_{r}-E(I_{r})\right) = \sum_{r=2}^{n}\left(C_{r}-\frac{F_{r-1}}{m}\right)\log\frac{m-F_{r-1}}{F_{r-1}}$

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

　 $H_{n}$ の和を構成する確率変数 $\left\{I_{r}-E(I_{r})\right\}$ が互いに独立でないことから，このままでは $H_{n}$ の確率分布について詳しく解析することができません．そこで， $F_{r-1}$ を期待値に置き換えて $H_{n}$ を近似したものを

$\displaystyle\hat{H}_{n} = \sum_{r=2}^{n}\left(\hat{C}_{r}-\frac{E(F_{r-1})}{m}\right)\log\frac{m-E(F_{r-1})}{E(F_{r-1})}, P(\hat{C}_{r}=1)=\frac{E(F_{r-1})}{m}, P(\hat{C}_{r}=0)=\frac{m-E(F_{r-1})}{m}$

とします．ここで， ${\hat{C}_{r}}$ が互いに独立であるとすると，級数の各項も独立となるので，その分散は

$\displaystyle V(\hat{H}_{n}) = \sum_{r=2}^{n}V\left(\left(\hat{C}_{r}-\frac{E(F_{r-1})}{m}\right)\log\frac{m-E(F_{r-1})}{E(F_{r-1})}\right) = \sum_{r=2}^{n}\frac{E(F_{r-1})\left(m-E(F_{r-1})\right)}{m^{2}}\left(\log\frac{m-E(F_{r-1})}{E(F_{r-1})}\right)^{2}$

となります．そして，期待値の一般式

$E(F_{r-1}) = m-m\left(\frac{m-1}{m}\right)^{r-1}$

を代入すると，

$\displaystyle V(\hat{H}_{n}) = \sum_{r=2}^{n}\left(\frac{m-1}{m}\right)^{r-1}\left(1-\left(\frac{m-1}{m}\right)^{r-1}\right)\left(\log\frac{\left(\frac{m-1}{m}\right)^{r-1}}{1-\left(\frac{m-1}{m}\right)^{r-1}}\right)^{2} = \sum_{r=2}^{n}\left(\frac{m-1}{m}\right)^{r-1}\left(1-\left(\frac{m-1}{m}\right)^{r-1}\right)\left(\log\left(\left(\frac{m}{m-1}\right)^{r-1}-1\right)\right)^{2}$

が導かれます．

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

　そして， $\hat{H}_{n}$ による近似が十分に良いのであれば，おそらく $H_{n}$ も正規分布に近い分布をもつであろうと予想されるので，その標本値を分散の近似値 $V(\hat{H}_{n})$ を用いて評価することが可能になります．これを，衝突エントロピー検定とよぶことにします．

３－３．乱数検定の実例

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

　まず，多次元分布の偏りを調べるために，点をプロットする空間を $6$ 次元の超立方体とします．各次元の座標には $4$ ビットずつを割り当て，セルの数を $m=2^{24}=16777216$ とします．各試行ごとに，乱数列 ${x_{r}}$ から値を $6$ 個ずつ取り出した組 $(x_{6r}, x_{6r+1}, x_{6r+2}, x_{6r+3}, x_{6r+4}, x_{6r+5})$ の，それぞれ通常は上位 $4$ ビットずつを用いることによって，点が入るセルを決定します．

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

　さらに， $H_{n}$ を $\sigma$ 値で評価するために，標準偏差を $\sigma=\sqrt{V(\hat{H}_{n})}$ で近似して比を求めます．そして，一応参考のために衝突回数の期待値からの偏差 $E(F_{n})-F_{n}$ を求め，同様に標準偏差 $\sigma=\sqrt{V(F_{n})}$ で割って $\sigma$ 値を計算します．

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

● Numerical Recipes ran2

　Numerical Recipes [1] で推奨されている ran2 とよばれる擬似乱数です．２つの周期の異なる線形合同法と切り混ぜ機構を組み合わせた擬似乱数で，約 $2.3\times 10^{18}$ の周期を持ちます．

$H_{n}=1091.23 \downarrow$

$F_{n}=0 \rightarrow$ $\leftarrow F_{n}=m/2=8388608$

$H_{n}=-2482.73 \uparrow$

　　 $\begin{array}{l}最大箪mathrm{l}―mathrm{F}　H_{n}=　1091.23=　1.66173\cdot\sigma, \\最十mathrm{I}箪mathrm{l}―mathrm{F}　H_{n}=-2028.11=-0.954386\cdot\sigma, \\最小箪mathrm{l}―mathrm{F}　H_{n}=-2482.73=-1.19837\cdot\sigma, \end{array}\begin{array}{l}n= 737874, F_{n}= 721659, \\n=11627974, F_{n}=8388608, \\n= 7393042, F_{n}=5980909, \end{array}\begin{array}{l}E(F_{n})-F_{n}=　224.193=　1.82566\cdot\sigma\\ E(F_{n})-F_{n}=-552.829=-0.487316\cdot\sigma\\ E(F_{n})-F_{n}=-1705.83=-1.9243\cdot\sigm;%<$

　グラフの横軸は $F_{n}$ ，縦軸は $H_{n}$ を表します．グラフの右端では $F_{n}=m/2$ となり，セルが半数埋まった状態に相当します．

　緑色の線は $H_{n}$ の変化を表していて，縦軸の範囲は $H_{n}$ の最大値と最小値で正規化されています．

　赤い線は $H_{n}=0$ を表していて，緑色の線がこの線より上にある場合は理想的なランダムより衝突が多く，下にある場合は衝突が少ないことを示しています．

　衝突エントロピー $H_{n}$ はビット単位で表され，青い線は $H_{n}=\pm 1024$ （bit）および $H_{n}=\pm\sigma=\pm\sqrt{V(\hat{H}_{n})}$ の範囲を示しています．

　この ran2 の場合，序盤に若干衝突が多く， $H_{n}$ が一旦 $1.66\cdot\sigma$ まで上昇しますが，その後は大体 $\pm\sigma$ の範囲に収まっていることから，この検定結果では特に異常はありません．

● 線形合同法　 $x_{r+1}\equiv 1664525\cdot x_{r}+1013904223$

　　 $\begin{array}{l}最大箪mathrm{l}―mathrm{F}　H_{n}=　143.656=　1.30224\cdot\sigma, \\最十mathrm{I}箪mathrm{l}―mathrm{F}　H_{n}=-20074.2=-9.44651\cdot\sigma, \\最小箪mathrm{l}―mathrm{F}　H_{n}=-20074.8=-9.44679\cdot\sigma, \end{array}\begin{array}{l}n= 74820, F_{n}= 74636, \\n=11602379, F_{n}=8388608, \\n=11506242, F_{n}=8340390, \end{array}\begin{array}{l}E(F_{n})-F_{n}=　17.4156=　1.35336\cdot\sigma\\ E(F_{n})-F_{n}=-13360.9=-11.789\cdot\sigma\\ E(F_{n})-F_{n}=-13426.2=-11.89\cdot\sigma\e;%<$

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

● 二次合同法　 $x_{r+1}\equiv 4\cdot x_{r}^{2}+5\cdot x_{r}+1$

　　 $\begin{array}{l}最大箪mathrm{l}―mathrm{F}　H_{n}=　1767.15=　0.989679\cdot\sigma, \\最十mathrm{I}箪mathrm{l}―mathrm{F}　H_{n}=　 696.95=　0.327971\cdot\sigma, \\最小箪mathrm{l}―mathrm{F}　H_{n}=-1041.51=-0.875398\cdot\sigma, \end{array}\begin{array}{l}n= 4208776, F_{n}=3721702, \\n=11630639, F_{n}=8388608, \\n= 1837778, F_{n}=1741072, \end{array}\begin{array}{l}E(F_{n})-F_{n}=　671.084=　1.13737\cdot\sigma\\ E(F_{n})-F_{n}=　779.653=　0.687191\cdot\sigma\\ E(F_{n})-F_{n}=-372.37 =-1.28566\cd;%<$

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

　　 $\begin{array}{l}最大箪mathrm{l}―mathrm{F}　H_{n}=　62.6422=　1.13363\cdot\sigma, \\最十mathrm{I}箪mathrm{l}―mathrm{F}　H_{n}=-4958.33=-2.33329\cdot\sigma, \\最小箪mathrm{l}―mathrm{F}　H_{n}=-5486.16=-3.33824\cdot\sigma, \end{array}\begin{array}{l}n= 32888, F_{n}= 32849, \\n=11624227, F_{n}=8388608, \\n= 3432244, F_{n}=3105640, \end{array}\begin{array}{l}E(F_{n})-F_{n}=　6.78723=　1.19742\cdot\sigma\\ E(F_{n})-F_{n}=-2426.66=-2.13939\cdot\sigma\\ E(F_{n})-F_{n}=-1711.13=-3.42266\cdot\sigma;%<$

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

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

　Ｃ言語式 x = (0x80000057 &−(x&1)) ^ (x>>1) で生成されるシフトレジスタ列を ${x_{r}}$ とし，二次合同式 $y_{r+1}\equiv(4\cdot y_{r}+1)(y_{r}+x_{r+1}+1)$ で延長した擬似乱数列です． $2^{64}-2^{32}$ の周期を持ちます．

　　 $\begin{array}{l}最大箪mathrm{l}―mathrm{F}　H_{n}=　1212.19=　1.08955\cdot\sigma, \\最十mathrm{I}箪mathrm{l}―mathrm{F}　H_{n}=-1399.15=-0.658413\cdot\sigma, \\最小箪mathrm{l}―mathrm{F}　H_{n}=-1931.3 =-0.930614\cdot\sigma, \end{array}\begin{array}{l}n= 1641549, F_{n}=1563489, \\n=11628531, F_{n}=8388608, \\n= 7483398, F_{n}=6038740, \end{array}\begin{array}{l}E(F_{n})-F_{n}=　308.606=　1.18135\cdot\sigma\\ E(F_{n})-F_{n}=-274.315=-0.241802\cdot\sigma\\ E(F_{n})-F_{n}=-1538.97=-1.72271\cdot\si;%<$

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

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

　なお，延長の式が前回扱った $y_{r+1}\equiv f(x_{r}, y_{r})$ ではなく $y_{r+1}\equiv f(x_{r+1}, y_{r})$ という形になっていますが，実はこの方がプログラムが簡単になり，また擬似乱数の性質も ${x_{r}}$ が一つずれるだけで，本質的には変わりません．

● Mersenne Twister MT19937

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

　　 $\begin{array}{l}最大箪mathrm{l}―mathrm{F}　H_{n}=　3955.12=　1.90517\cdot\sigma, \\最十mathrm{I}箪mathrm{l}―mathrm{F}　H_{n}=　3841.43=　1.8077\cdot\sigma, \\最小箪mathrm{l}―mathrm{F}　H_{n}=-566.347=-1.33149\cdot\sigma, \end{array}\begin{array}{l}n= 7501457, F_{n}=6046175, \\n=11635257, F_{n}=8388608, \\n= 407271, F_{n}= 402458, \end{array}\begin{array}{l}E(F_{n})-F_{n}=　2580.37=　2.88404\cdot\sigma\\ E(F_{n})-F_{n}=　3088.12=　2.72141\cdot\sigma\\ E(F_{n})-F_{n}=-90.5311=-1.31393\cdot\si;%<$

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

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

　　 $\begin{array}{l}最大箪mathrm{l}―mathrm{F}　H_{n}=　1133.72=　0.715551\cdot\sigma, \\最十mathrm{I}箪mathrm{l}―mathrm{F}　H_{n}=　353.715=　0.166451\cdot\sigma, \\最小箪mathrm{l}―mathrm{F}　H_{n}=-1134.88=-1.23223\cdot\sigma, \end{array}\begin{array}{l}n= 3164414, F_{n}=2883400, \\n=11629035, F_{n}=8388608, \\n= 1214503, F_{n}=1171860, \end{array}\begin{array}{l}E(F_{n})-F_{n}=　498.107=　1.06647\cdot\sigma\\ E(F_{n})-F_{n}=-22.3108=-0.019666\cdot\sigma\\ E(F_{n})-F_{n}=-274.107=-1.38856\cdot\;%<$

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

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

　Mersenne Twister を ${x_{r}}$ とし，二次合同式 $y_{r+1}\equiv(4\cdot y_{r}+1)(y_{r}+x_{r+1}+1)$ で延長した擬似乱数列です． $2^{19969}-2^{32}$ の周期を持ち， $624$ 次元にわたってほぼ均等に分布することが理論的に証明されます．

　なお，Mersenne Twister ${x_{r}}$ はデフォルトの初期値で生成し，延長乱数列は $y_{-1}=0$ という初期条件で検定を行いました．

　　 $\begin{array}{l}最大箪mathrm{l}―mathrm{F}　H_{n}=　2413.78=　1.21356\cdot\sigma, \\最十mathrm{I}箪mathrm{l}―mathrm{F}　H_{n}=　1810.79=　0.852121\cdot\sigma, \\最小箪mathrm{l}―mathrm{F}　H_{n}=-597.184=-0.535592\cdot\sigma, \end{array}\begin{array}{l}n= 5973245, F_{n}=5024530, \\n=11632842, F_{n}=8388608, \\n= 1647508, F_{n}=1569444, \end{array}\begin{array}{l}E(F_{n})-F_{n}=　1095.2 =　1.42638\cdot\sigma\\ E(F_{n})-F_{n}=　1880.98=　1.65777\cdot\sigma\\ E(F_{n})-F_{n}=-243.79 =-0.930129\cdo;%<$

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

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

　　 $\begin{array}{l}最大箪mathrm{l}―mathrm{F}　H_{n}=　2066.34=　1.24866\cdot\sigma, \\最十mathrm{I}箪mathrm{l}―mathrm{F}　H_{n}=-511.946=-0.240911\cdot\sigma, \\最小箪mathrm{l}―mathrm{F}　H_{n}=-1093.27=-0.982056\cdot\sigma, \end{array}\begin{array}{l}n= 3487167, F_{n}=3147726, \\n=11627519, F_{n}=8388608, \\n= 1643231, F_{n}=1565571, \end{array}\begin{array}{l}E(F_{n})-F_{n}=　891.456=　1.7598\cdot\sigma\\ E(F_{n})-F_{n}=-780.347=-0.687883\cdot\sigma\\ E(F_{n})-F_{n}=-248.249=-0.949408\cdot\si;%<$

　途中若干衝突が多くなりますが，最終的には元の 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 ［日本語版］Ｃ言語による数値計算のレシピ』（技術評論社，1993）．

[2] Mersenne Twister: A random number generator

次に進む

戻る