数理解析研究所講究録
1011 巻 1997 年 77-88

77

無理数回転を利用した擬似乱数生成法–

確率論による数値解析的アプローチ

杉田九州大学大学院数理学研究科
Abstract. Lebesgue 確率空間 $(\Omega,P),$ $\Omega=[0,1),$$P=\mathrm{L}\mathrm{e}\mathrm{b}\mathrm{a}\mathrm{e}1^{\mathrm{e}}$洋 (Hiroshi SUGITA)

測度, の上で定義された $\{0,1\}$ -値定常過程
$m$$x_{n}^{(m)}..( \omega):=\sum_{i=1}d_{t}(\{\omega+n\alpha\})$($\alpha$は無理数, {X} は実数$x$(mod 2),$\{0$,$n=0,1,$$\ldots$の 2 進小数展開における第 桁を表す) は $marrow\infty$,1}-値独立同分布確率変数列) に分布収束する。我々は $m$ をの小数部分, また $d_{i}(x)$ はのとき硬貨投げの確率過程 (平均 1/2 の$\omega\in\Omega$

$i$

$x$

十分大きく取ったとき, この定常過程のサンプルパスを擬似乱数として用いることを提唱する。実用には$\alpha=(\sqrt{5}-1)/2$

1のとき $m=90$ 程度で十分である。

序「擬似乱数生成法」 というと, 数学的に純粋でない「まがい物」, という印象を持つ人が少な,「計算機プログラムによる乱数生成
くない。「乱数」をどのように定義しようともは実現不可能である」 ことには間違いなく, 従って実際に計算機プログラムによって生成される
擬似乱数は確かに 「乱数の近似物」 にすぎない。
$(\mathrm{C}\mathrm{f}.[14,5,13])$方で「近似物』が立派に通用している例もある。 たとえば, 無理数や無限級数なども, 有限の記憶領域しか持たない計算機では実現不可能であって, 実際は無理数を有理数で近似して, あ
るいは無限級数を有限和で近似して計算する。 こうした近似を 「まがい物」 と思う人はほとんどいない。数値解析では「近似」 と言うとき, 次の二つの条件を要請する:(A.1) 誤差について定量的に述べることができること。
(A 2) 誤差をいくらでも小さくするアルゴリズムが存在すること。
これらの要請が満たされれば, 実際的な計算において 「真の値」 と「近似値」 の間に機能的な差はなくなる。無理数の近似有理数あるいは無限級数の近似有限和はまさにこれらの条件を満たしている。 -方, 現在主流の各種の擬似乱数は残念ながら数値解析における上の二つの要請を満たしているとは言い難い。
.
..
本稿では, 擬似乱数を数値解析における近似のレベルまで引き上げる–つの実用的な手法について述べる。すなわち, 我々は擬似乱数を力学系によって定まる定常過程として位置付けし, その誤差を確率論の言葉でもって表し, さらにその誤差が に収束するような擬似乱数生成アルゴ$0$

リズムの「列」を構成する。

本稿の構成を明らかにしておく。第 2 節では力学系による擬似乱数生成の–般論を述べ, 我々の方法の位置付けを行う。 第 3 節では我々の擬似乱数の定義と主定理を述べる。 第 3 節では我々の擬似乱数の多次元分布を求めるためのアルゴリズムを紹介する。第 4 節ではその前の節で与えたアルゴリズムによって多次元分布を求め, その–様性に関して評価する。第 5 節は我々の擬似乱数を生成するための によるプログラム例を掲載した。
$\mathrm{C}$2$\text{力学系によ条擬似乱数生成}$
,

21

.

$,.-$

$\cdot$

硬貨投げの確率過程

78

.

.

$\mathrm{r}:..\cdot$$i$.

$\cdot$$,\mathrm{e}_{\wedge}.$

乱数のモデルとしては原理的にも実際的にも 2 点集合 $\{0,1\}$ に値をとる平均 1/2 の独立同分布, すなわち硬貨投げの確率過程, を考えれば十分である。硬貨投げの確率過確率変数列$\gamma$$\{X_{n}\}_{n=0}^{\infty}$

程はいろいろな性質を持つが, ここでは次の三つの性質に注目する。

-値強電過程である。すなわち, 任意の $k,$と $(X\iota, \ldots,x\iota+k-1)$ が同分布である o.

$(\mathrm{B}^{-}.1)\{0,1\}$$l\in \mathrm{N}$(B.2) エルゴード的である。 とくに, 確率 1 で, 任意のに対して Rk-値確率変数 ($k\in \mathrm{N}$と任意の$X_{0},$$\epsilon\in\{0,1\}^{k}$$\ldots$,Xk-l)に対して 1,

Probab. $((x_{0}, \ldots,xk-1)=\epsilon)=$ hm $\underline{1}\#\{0\leq n\leq N-1|(X_{n}X_{n}k_{-}1)\mathrm{t}^{k-1})’\cdots,=\epsilon\}$$Narrow\infty N$(B.3) 確率 1 で見本は \nu 進正規列である。すなわち, 確率 1 で, 任意のに対して$k\in \mathrm{N}\text{と任意^{の}}$ $\epsilon\in\{0,1\}^{k}$$\lim\underline{1}\#\cdot\{0\leq n\leq N-1|(x_{n\mathrm{t}^{k-}}1)’\ldots,xnk-1)=\epsilon\}=2-k$$Narrow\infty N$

となる。
.

.

.

.

また, 以上三つの性質を持つ確率過程は硬貨投げの確率過程に限られる。以下, 本稿で扱う擬似.

乱数は硬貨投げの確率過程をモデルとするものとする。

22

定常過程としての擬似乱数

擬似乱数を計算機で生成しようとするとき, 数学的には力学系と呼ばれる枠組み (集合$\Omega$

とそ-回擬似乱数のの組のこと) を利用するのが一般的である。
の上の変換$\Omegaarrow\{0,1\}$$f$
および初期値 (擬似乱数の「種」 と言うことも多い), 写像 :
場合は, 力学系
を適当に設定して,$T:\Omegaarrow\Omega$$\{0,1\}$$(\Omega,T)$$(\Omega, T)$$\omega\in\Omega$$X_{n}(\omega):=f(T^{n}\omega)$と定義される数列
そこで問題は「力学系義される数列,$n=0,1,$(1)$\ldots$

として得られる。
をどのように設定すれば, (1) で定, 写像 $f$ および初期値が良い擬似乱数になるか」 ということである。

$\{X_{n}(\omega)\}^{\infty}n=0$$\omega\in\Omega$$(\Omega, T)$$\{X_{n}(\omega)\}_{n=}^{\infty}0$

はユーザが自由に指定できるようにしておく。 これは–つの擬似乱数生成
アルゴリズムによって様々な擬似乱数系列を生成できるようにするためである。確率論的には初期値 をある上の確率測度 $P$ に従って選択するという状況を想定するのが自然だろう。そう症の確率過程と見なされ, 乱数のモデル「硬貨投げ
は確率空間
すれば数列
.
,:
の確率過程」 と同じ土俵に乗ることができる。
確率測度 $P$ は原理的にはどのようなものを仮定してもよい。 しかし, 様々な計算をうまく行うために $T$-不変であると仮定するのがよい。 このとき,は強定常過程となって条件
さて, 初期値$\omega$$\omega\in\Omega$$\Omega$$.\{X_{n}(\omega)\}^{\infty}n=0$$(\Omega,.P)$

$-$

$r$

$\{X_{n}(\omega-.)\}_{n=}^{\infty}0$(B. 1) を満たす。
.
.
.
また, 擬似乱数のサンプルからその分布を推定できるという実際的な仮定は是非とも必要だか
ら, エルゴード性 (B.2) も仮定しよう。
$-$$1\#\{\cdots\}$

は集合 ...
$\{$$\}$の要素の個数を表す。

79

2.3

-意エルゴード性以上の仮定の下で (B.3) を満たすようにできる, すなわち初期値の選択だけがランダムに行われるような力学系による確率過程で硬貨投げの確率過程を実現できる。
le (Borel の例) $\Omega=[0,1),$とする 2。このとき, 関数 $f$ を$1.\mathrm{E}\mathrm{x}\mathrm{a}\mathrm{m}\dot{\mathrm{p}}$$P=\mathrm{L}\mathrm{e}..\mathrm{b}\vee$aele 測度, さらに変換 $T$ を,
1,$\omega<1/2$$0$$f(\omega):=d_{1}(\omega)=\{$とすれば, 確率空間 $(\Omega, P)$ 上の確率過程$T\omega:=\{2\omega.\},.\omega\in\Omega$$\omega\geq 1/2$$\{X_{n}(\omega)\}_{n=0}\infty$$X_{n}(\omega)=f(T^{n}\omega)$,$\omega\in\Omega$,$n=0,1,2,$$\ldots$

は硬貨投げの確率過程の–つの実現である (たとえば [2] の最初の十数ページを見よ)。
を \nu 進有限小数に設定せざ
もっとも, Borel の例は乱数生成には使えない。計算機では初期値に収束してしまう。 この難点を乗り越えるを得ないが, それだとすべての軌道がたちどころに$\omega$

$0$

るために,

香田氏らは 2 次以上の多項式を用いたカオス的力学系による擬似乱数生成法を考案している $([9, 11]^{3})$ 。しかしながら,カオス的力学系では長い軌道の追跡が数値計算上不可能であるの場合を除いて絶望的である。
ため, 性質 (B.2) が数値的に成り立つことを保証するのは小さい
たとえ長い軌道の追跡が可能になったとしても, カオス的力学系では性質$k$$(\mathrm{B}.2)^{\text{を満たさ}}.\text{ない初}$期値$\omega$が必ず存在することにも注意しなければならない。

そのため, ここではカオス的力学系を擬似乱数生成には採用しない。そのかわり, 任意の初期値$\omega$に対して性質 (B.2) を要請する。 そのためには–意エルゴーード性:(B.4) 任意の有界連続関数$g:\Omegaarrow \mathrm{R}$

.

に対して

$N arrow\infty \mathrm{h}\mathrm{m}\frac{1}{N}\sum^{N1}g(\tau^{n}\omega n=0-)=\int_{\Omega}g(\omega)\prime P(\omega’)$,$\forall\omega\in\Omega$(2)を仮定しよう。 このとき, 初期値 の選び方によって擬似乱数の統計的性質は左右されないので大変都合がよい。 しかし, この場合は理論的にも決して硬貨投げの確率過程は得られないことに注意しよう。
$\omega$24

擬似乱数の誤差と統計的検定

我々の定式化では擬似乱数は定常過程であるから, 硬貨投げの確率過程からの隔たり, すなわ
ち擬似乱数の誤差, を確率論の楓念によって述べることができる。それは「擬似乱数の分布と硬貨投げの確率過程の分布の距離」 と定義すればよい。一般に位相空商 $S$ 上の確率測度全体には標準的な位相として「弱収束の位相」 と呼ばれる位相が知られている。それは, すべての連$\mathcal{M}_{1}(S)$続関数 $f:Sarrow \mathrm{R}$ に対して, 写像$\mathcal{M}_{1}(S)\ni\mu\mapsto\int_{S}f(s)\mu(S)\in \mathrm{R}$はの小数部分を表す。
2 実数 に対して
3 本講究録の香旧氏の記事も参照せよ。
$x$$\{x\}$$x$80

が連続になるような最弱の位相である。 とくに, 可分完備距離空間上の「弱収束の位相」は距離付
け可能であり, その位相を与える距離を–般に Prohorov の距離 (詳しくは [1] を見よ) という。 そ上の Prohorov の距離で計る。
こで擬似乱数の誤差はProhorov の距離をもって擬似乱数の誤差と定義する理由は, 擬似乱数の統計的検定と関係がある。実際, 擬似乱数の善し悪しを見極める各種の統計的検定は, 擬似乱数の分布が硬貨投げの確率過程とどれくらい離れているかを確かめる作業であると言ってよい。
.
に硬貨投げの
一般に統計的検定とはどういうものかを述べよう。まず, 関数 4については, その分布が良く分かっを代入して得られる確率変数
確率過程ているとする。 そこで $\{0,1\}$ -値擬似乱数の–つのサンプル{xn};。を同様に代入して得らは擬似
を X のサンプルと仮定したときにおよそ起こり得ない値ならば,が,れた値を擬似乱数として採択す乱数としてふさわしくないと判断し, また, 十分起こり得る値ならばる。 このような作業が擬似乱数の検定である。
我々の定式化の下では, $\{0,1\}$ -斗強定常過程である擬似乱数 X’{Xn’}n\infty =。を同様に代入し擬似乱数のサンプルが上記の意味で採択されの分布に近ければ近いほど,の分布がる確率が高くなる。すなわち, Prohorov の距離で測った誤差が小さい擬似乱数は各種統計的検定に合格する確率が高くなる。
の分布が各サンプルの値は何らかのサンプル平均であることが多いので,実際にはの漸近的相対度数に遺伝していること, すなわち, –意エルゴード性 (B.4) が成り立つと都合が$\mathcal{M}_{1}(\{0,1\}^{\infty})$$F:\{0,1\}^{\infty}arrow \mathrm{R}$$F(\mathrm{X})$$\mathrm{X}=\{X_{n}\}_{n}\infty=0$$\mathrm{x}=$$F(\mathrm{x})$$\mathrm{x}$$\mathrm{x}$$\mathrm{x}$$=$$F(\mathrm{X}’)$$F(\mathrm{X})$$\mathrm{X}’$$F(\mathrm{X}’)$

よい。

25

問題設定

–

確率論による数値解析的アプローチ

以上の概念を準備すれば擬似乱数生成法の問題設定を述べることができる。序で述べたようにの条件を満たすようにしなけれ擬似乱数を数値解析的レベルまで引き上げるためには (A. )ばならないが, 前節までで (A.1) はクリアした。 (A 2) のために, 次のように問題を設定しよう。
$1$

を確率空間, $T$ をは次の性質を持つ写像列$(\Omega,P)$$\Omega$$f^{\langle m)}$上の P-不変な変換で–意エルゴード的であるとする。そこで問題: $\Omegaarrow\{0,1\}$ を構成することである:$x_{n}^{(m)}(\omega):=f^{(m)}(T^{\iota}\omega)$

で定まる強定常過程の列収束する。

$(\mathrm{A}.2)$,$\{X_{n}^{\mathrm{t}}m)\}^{\infty}n=0’ m\in \mathrm{N}$$\omega\in\Omega$,$n=0,1,$$\ldots$(3), が $marrow\infty$ のとき, 硬貨投げの確率過程に分布

このような写像列 $f^{(m)}$ を構成できれば理論的には (A 2) をクリアできる。擬似乱数としては十分大きな $m$ に対して (3) で定まるものを採用すればよい。 しかし, 実用的な観点からはそのよう.
な写像 $f^{(m)}$ が計算機によって高速に計算されることが望ましい 5。
の積位相に関して連続であるとしてよい。
4 現実的な検定では有限個のサンプルしか扱わないから $F$ は
5 実際, [7] では我々と同様の枠組みの下で Gauss 型独立同分布確率変数列に分布収束するようなものを実現しているが収束の早さが遅い。そのため, 擬似乱数生成に利用しようとすると膨大な計算量を必要とするので実用的でない。
$\{0,1\}^{\infty}$

81

3

無理数回転を利用した擬似乱数生成

3.1

定義と極限定理

それでは, 標題の擬似乱数の生成アルゴリズムを述べよう。 $(\Omega, P)$ を Lebesgue 確率空間, すはの 2 進小数展開におなわちに対して測度とする。
$\Omega=[0, \mathrm{L}),$

ける小数点以下第定常確率過程$\omega\in\Omega$$P=\mathrm{L}\mathrm{e}\mathrm{b}\mathrm{e}\mathrm{S}\mathrm{l}\mathrm{e}$$i$

桁を表すこととし, 各を次のように定める。

$m\in \mathrm{N}$$d_{i}(\omega)$に対して$(\Omega, P)$$\omega$上で定義された $\{0,1\}$-値強$\{X_{n}^{(m)}\}^{\infty}n=0$2Definition(mod 2),$x_{n}^{(m)}( \omega):=\sum_{i=1}^{m}d_{i}(\{\omega+n\alpha\})$ここに$\alpha$は無理数であり, また記号 $\{x\}$ は実数$x$$\omega\in\Omega$,$n=0,1,$$\ldots$(4)の小数部分を表す。

このとき次の主定理が成り立つ。

は $marrow\infty$ のに対して 6, 確率過程
ほとんどすべての無理数
とき硬貨投げの確率過程に分布収束する。すなわち, 任意のと任意の’0h . . . , $\epsilon_{k-1}\in\{0,1\}$

3.

$\mathrm{T}\mathrm{h}\mathrm{e}\mathrm{o}\mathrm{r}\mathrm{e}\mathrm{m}$$([18])$$\{X_{n}^{(m}\}^{\infty})n=0$$\alpha$$k\in \mathrm{N}$に対して次式が成り立つ。
$\ldots,X_{k-}^{(}(m_{1}))\omega=\epsilon_{k-1})=2^{-k}$$\lim_{marrow\infty}P(X_{0}^{(m)}(\omega)=\epsilon_{0},$無理数回転の–意エルゴード性 (変換プルパスについて次の定理が成り立つ。

$T:\omega\mapsto\{\omega..\cdot+..,\alpha.\}$4.
$\epsilon_{0},$$\ldots,$$=$

ほとんどすべての無理数 , 任恵の$\epsilon_{k-1}\in\{0,1\}$ に対して次式が成り立つ。

$\mathrm{T}\mathrm{h}\mathrm{e}\mathrm{o}\mathrm{r}\mathrm{e}\mathrm{m}$him$\alpha$hhm

が (B.4 ) を満たすこと) によって, 各サン$\omega\in\Omega,$.

任意の$k\in \mathrm{N}$および任意の$\underline{1}\#\{0\leq n\leq N-1|X_{nk}^{(m)}(\omega)=\epsilon_{0}, \ldots , X_{nk+k}^{(m})-1(\omega)=\epsilon_{k-1}\}=2^{-k}$$marrow\infty Narrow\infty N$

は, $m$ が大きいとき擬似乱数
4.
.
均等に分布するということを意味する。
$\mathrm{T}\mathrm{h}\mathrm{e}\mathrm{o}\mathrm{r}\mathrm{e}\mathrm{m}$$\{X_{n}^{(m.)}*.\cdot.

\cdot(.\cdot\omega\backslash .)\}$の各サンプルパスは多次元にわたってほぼ.

測度に関して。事実としてはすべての無理数で成り立つと思われるが, 厳密な証明のためには現在のとに必要。なお, 3 次元以下の周辺分布はすべての無理数に対して硬貨投げの確率過程のそれにころ技術的な条件が収束する。
$,\mathrm{L}\mathrm{e}\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{g}\mathrm{u}\mathrm{e}$$\alpha$82

3.2

多次元周辺分布の計算アルゴリズム

2Definition で与えた強定常確率過程$P(X^{(m)}=\epsilon 00’\ldots$$\{X_{n}^{\mathrm{t}^{m})}(\omega)\}_{n=}^{\infty}0$, $x^{(m)}=\epsilon-1)$$k-1k$ ,の多次元周辺分布$k\in \mathrm{N}$,$\epsilon_{0},$$\ldots,$$\epsilon_{k-1}\in\{0,1\}$(5)について考えよう。我々の擬似乱数では (5) をすべて算出するためのアルゴリズムが存在する。

5.

$\mathrm{L}\mathrm{e}\mathrm{m}\mathrm{m}\mathrm{a}$$([18])$(i) 多次元周辺分布 (5) は次の量から算出することができる。

$E_{0,k,\ldots,k\iota_{-}1;0}^{(m_{1}}:=)\mathrm{d}\mathrm{d}P(x^{\mathrm{t}m)}+0x+k_{1}+x_{kl-1}^{\mathrm{t}m)}=\text{奇}\mathrm{t}m)\ldots \text{数})$$0<k_{1}<\ldots<k_{l-1}$(皿)$l\in \mathrm{N}$が奇数ならば

$E_{0}^{\langle)},m_{k_{1},\ldots,k\iota_{-1;}\mathrm{o}\mathrm{d}\mathrm{d}^{=1}}/2$,,$l\in \mathrm{N}$である。

を求めるアルゴリズムがあればよい。そのため
5Lemma から が偶数のときにに対して,に, いくつかの記号を導入しなければならない。無理数回転で用いられる無理数$l$$E_{0}^{(m)},k_{1},\ldots,k1-1;\mathrm{o}\mathrm{d}\mathrm{d}$$\alpha$$\alpha_{\mathrm{j}}:=\{k_{j}\alpha\}$,$j=1,$$\ldots,$( は偶数)$l-1$ ,

$l$

とおき,
$($

$\alpha_{\dot{\tau}}^{\mathrm{t}m)}L$$:=$$[2^{m}\alpha_{j}]/2^{m}$$:=$$\{[2^{m}\alpha_{j}+1]/2^{m}\}$$:=$$2^{m}(\alpha_{j}-\alpha_{j}^{\mathrm{t}})m)L$I は整数部分を表す。次にとする。 ここに, 記号

上の置換$\{1, \ldots, l-1\}$に定める。
$1>\beta_{\sigma}^{\mathrm{t}^{m)\mathrm{t}m)}}\mathrm{t}m,1)\geq\beta\sigma \mathrm{t}m,2)\beta^{(}\sigma\geq\ldots\geq\iota-1)\geq m)\mathrm{t}m,0$

ただし, 便宜上$\beta_{\sigma(}^{(m)}m,0$)$:=1,$)$\beta_{\sigma(}^{(m}m,\iota$

)

$:=0$

と約束しておく。 そして$\alpha_{j}^{\langle m)U}$$\alpha_{j}^{(m)_{S}}’:=\{$$\alpha_{j}^{(m)L}$, if $\sigma(m,j)\leq s$, if $\sigma(m,j)>s$とした上で$\alpha^{(m)_{S}}’:=(\alpha_{1}^{(m),s}, \ldots, \alpha_{\iota_{-}1})\mathrm{t}^{m}),S$,$s=0,1,$$\ldots,$$l-1$

とおく。最後に, 2 進有限小数の集合を以下のように定義する。
$D:= \bigcup_{m\in \mathrm{N}}\{\frac{n}{2^{m}}\in[0,1)|n=0,$$\ldots,$$2^{m}-1\}$

6Theorem ([18])レ$E_{0}^{\mathrm{t}^{m_{k_{1},\ldots,k}}},)l-1$

; odd

)

),$)s+1$)$-\beta\sigma \mathrm{t}^{m}$$= \sum_{s=0}(\beta_{\sigma(m,S}^{\langle m)\langle m}$$(\alpha^{\mathrm{t}m),s}B)$$\sigma(m, \cdot)$

を次のよう

83

レ1

ここに $B(\cdot)$ は$D^{l-1}=\overline{D\cross\cdots \mathrm{x}D}$の上で定義されたある実数値関数で, $B(\alpha^{(m),S})$ の値は,$B(\alpha^{\mathrm{t}^{0}),s})=0$$\cdot$$s=0,1,$$\iota-1$$\ldots,$

および次の漸化式で計算される。
$\frac{1}{2}B(\alpha^{(m-1)}’)S_{2}+\frac{1}{2}B(\alpha^{\mathrm{t}^{m-}})1,S1+s_{2})$if,$s_{1}$

is even,$B(\alpha^{(m),S})=\{$$\frac{1}{2}$ただし,$s_{1},$$s_{2}$(1-$B( \alpha^{(1),S}-2)m)+\frac{1}{2}(1-B(\alpha^{()_{S}}-1,1+s2)m)$$s_{1}$

is odd.

は次で与えられる。
$s_{1}:=j1 \sum_{=}^{l-1}dm(\alpha)\mathrm{t}jm),S$

4

, if

,

$s_{2}:= \sum_{j=1}dSm(\alpha)\sigma \mathrm{t}^{m}\dot{s})$多次元周辺分布の事前評価を用いると統計的検定を待つまでもなく, 多次元周辺分布に関する統計的性質を調べることができる。以下に挙げる例では, 無理数回転に用いる無理数と我々の擬似乱数の場合は$6.\mathrm{T}\mathrm{h}\infty \mathrm{r}\mathrm{e}\mathrm{m}$して黄金分割の比として知られる次の数を採用した 7。
$\alpha=\frac{\sqrt{5}-1}{2}$412 項間相関はじめに, 我々の乱数の 2 項間の相関について.
$a^{(m)}(K):=1 \leq k\leq K\mathrm{m}\mathrm{a}\mathrm{x}|E_{0,k;\mathrm{o}\mathrm{d}}^{(m}-)\frac{1}{2}\mathrm{d}|$と定義する。

$N_{\mathrm{c}}^{\langle m)}(K)$$N_{c}^{(m)}(K):= \frac{1}{16(a^{\langle m})(K))^{2}}$(6)を臨界サンプル数 8 と呼ぶ。次の各々の帰無仮説
1$E^{(m)}$0, ;odd$k$$=\overline{2}$’$k=1,.2,$ . $.,$ $,K$ ,以下ならば, 上の各々の仮説はそに関して検定 (危険率 5%) を行うとき, サンプル数がれぞれ 93% 以上の確率で採択されることが期待できる ([181)。
のときに計算し, 表にしたのが 7Table9 である。 $a^{(m)}(K)$ の値のすぐ右側(6) を$N_{\mathrm{c}}^{\langle m)}(K)$$K=1\mathrm{O}\mathrm{O}\mathrm{O}\mathrm{O}$の $()$ の中は最大値がどのような $k$ によって達成されたかを表わす。

の主張を成り立たせる無理数かどうか, 筆者は厳密な回答ができない。
および
7 この数が
しかし, 実用上は 2 進小数展開で 100 桁程度あればよく, 厳密な議論は大して問題にならない。
の 4 倍である。
8[18] で述べられている critical sample number はここでののときほぼ指数的に減少することが読み取れ$b^{(m)}(16)$$marrow\infty$が) とと 8Table において97.
についてそのことが示せる
る。実際, 理論的にもほとんどすべての無理数$3.\mathrm{T}\mathrm{h}\infty \mathrm{r}\mathrm{e}\mathrm{m}$$4.\mathrm{T}\mathrm{h}\infty \mathrm{r}\mathrm{e}\mathrm{m}$$N_{\mathrm{c}}^{(m)}(K)$$\mathrm{T}\mathrm{a}\mathrm{b}\mathrm{l}\mathrm{e}$$a^{(m)}\mathrm{t}10\mathrm{o}\mathrm{o}0$$a$$\langle[181)\circ$84

7Table0.4860680 (5473)
0.1084934 (1449)(305)
0.0435756
0.0029834 (305)0.0001943 (610)0.0000136 (8484)

10
$20$
$30$

$40$ .
$50$
$60$

4.2

$70$

$1.2\cross 10-\epsilon$

$80$

2 Ox

$90$

$8.5\cross 10-9$$100$$2.9\cross 10^{-}9$$2.6\cross 10-1$$5.3\cross 10^{0}$$3.3\cross 10^{1}$$7.0\cross 1.\mathrm{o}^{3}$$1.7\cross 10^{6}$$3.4\cross 10^{8}$(7264)(7697)(165)$10^{-7}$$4.1\cross 1010$$1.6\cross 10^{1}2$$8.7\cross 10^{1}4$(5201)$7.7\cross 10^{1}5$

多項間相関

次に–般の多次元周辺分布 ( $K$-次元以下) の評価を考えよう。 このとき各偶数$E^{(m)}$$0,k_{1},\ldots,k_{l-1;}$

odd’,

$l$

について

$1\leq k_{1}<\cdots<k\iota-1\leq K$を評価すればよい。 これらを全部調べることは比較的小さな $K$ についてさえ計算量が莫大にな
り大きな $K$ では絶望的に思えるが, それでも少しは望みがある。 8Table は $K=16$ の場合を計算したものである。左の欄は,$b^{(m)}(16):=_{1\leq k_{1}<\ldots k} \max<\iota-1\leq 16|E_{0,k1,k,\ldots,k_{l-1;}}^{\langle m})2$を表わし, 右の欄は最大値がどのような$k_{1},$$\ldots$odd$- \frac{1}{2}|$,によって達成されたかを表わす。 8Table からは次の仮説が成り立つように見受けられる。

$9.\mathrm{H}\mathrm{y}\mathrm{P}^{\mathrm{o}\mathrm{t}}\mathrm{h}\mathrm{e}\mathrm{s}\mathrm{i}\mathrm{s}$各$K\in \mathrm{N}$に対して $m$ が十分大きいとき,$1 \leq k_{1}<\cdots<kl-1\leq K\mathrm{m}\mathrm{a}\mathrm{x}|E_{0,k1k\mathrm{I}}^{\mathrm{t}m}),\ldots,--1;\mathrm{o}\mathrm{d}\mathrm{d}\frac{1}{2}|=1\leq k\leq K\mathrm{m}\mathrm{a}\mathrm{x}|E_{0,k;\mathrm{o}\mathrm{d}}^{(m)}-\mathrm{d}\frac{1}{2}|$

実は 9.

は

の証明を詳しく見ると成り立つことが十分期待できるのであるが現在のところ厳密な証明はない。 もし 9H\psi othes おが正しければ, 我々は 2 項間の相関の最大値さえ評価すればよいことになる 10。

5$\mathrm{H}\mathrm{y}\mathrm{p}\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{s}\mathrm{i}\mathrm{S}$$3.\mathrm{T}\mathrm{h}\infty \mathrm{r}\mathrm{e}\mathrm{m}$プログラミング

5.1$\mathrm{C}$によるプログラム例

数値解析の常として, 精度を高めるためには計算時間が長くなり, 計算時間を短くするためには精度を落とさなければならない。実用に供するためには精度が高ければそれでよい, というわ.

けには行かない。
$1\mathrm{o}_{9.\mathrm{H}\mathrm{y}}\mathrm{P}^{\mathrm{o}\mathrm{t}}\mathrm{h}\mathrm{e}\mathrm{s}\mathrm{i}\mathrm{s}$が正しければ$3.\mathrm{T}\mathrm{h}\infty \mathrm{r}\mathrm{e}\mathrm{m}$および$4.\mathrm{T}\mathrm{h}\infty \mathrm{r}\mathrm{e}\mathrm{m}$がすべての無理数に対して成り立つこともわかる。

85

ここでは $\alpha=(\sqrt{5}-1)\backslash /2$ および $m=90$ の場合に, 我々の擬似乱数を生成するための$\mathrm{C}$によるプログラムの例を挙げる。前節の 7Table で見るように擬似乱数の精度としては実用上十分であると期待される。 これは後述のように生成速度の観点からも十分実用になる。
$/\subset===========.==---===---=.=.================^{-}--===================*/$$/*$

Inplementation of Pseudo-random number generator by$/*$$\mathrm{m}90$-method with the irrational number(sqrt (5) $-1$ ) $/2$

$*/$

.

$*/$

$/\subset==================_{---}--------------------------------------*/$#include#define LIMIT#define CARRY$<\mathrm{s}\mathrm{t}\mathrm{d}\mathrm{i}\mathrm{o}.\mathrm{h}>$$\mathrm{O}\mathrm{x}3\mathrm{f}\mathrm{f}\mathrm{f}\mathrm{f}\mathrm{f}\mathrm{f}\mathrm{f}$Ox40000000
omega[5]:unsigned longvoid$/*\mathrm{C}\mathrm{u}\mathrm{r}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{t}$$\mathrm{m}90\mathrm{S}\mathrm{e}\mathrm{t}\mathrm{S}\mathrm{e}\mathrm{e}\mathrm{d}\mathrm{S}(\mathrm{S}\mathrm{o}.

unsigned long$\mathrm{s}0.\mathrm{s}\mathrm{l}.\mathrm{S}2,\mathrm{S}3.\mathrm{S}4$$\mathrm{o}\mathrm{m}\mathrm{e}\mathrm{g}\mathrm{a}[0]=\mathrm{s}0$.omega\mathrm{s}3.

$/*$\mathrm{s}4)$

Initialization $*/$

;

.

$\{$

$[..3]$\mathrm{s}1_{*}\mathrm{s}2.

seeds $*/$

& LINIT;$\mathrm{o}\mathrm{m}\mathrm{e}\mathrm{g}\mathrm{a}[1]=$s3 & LIMI.T;$=$$\mathrm{o}\mathrm{m}\mathrm{e}\mathrm{g}\mathrm{a}[4]=$sl &LIMIT;$\mathrm{o}\mathrm{m}\mathrm{e}\mathrm{g}\mathrm{a}[2]=$s4 &LINIT;

s2 &LIMIT:

$\}$

char

$\mathrm{m}90\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{o}\mathrm{m}\mathrm{b}\mathrm{i}\mathrm{t}$

$/*$

$()$

Returns$0$

or 1 at random $*/$$\{$

static unsigned long;char data-byte;
union bitarray.

{$\mathrm{a}\mathrm{l}\mathrm{p}\mathrm{h}\mathrm{a}[5]=$$0\mathrm{x}\mathrm{l}7\mathrm{f}4\mathrm{a}7\mathrm{c}\mathrm{l}$$0\mathrm{x}278\mathrm{d}\mathrm{d}\mathrm{e}6\mathrm{e},$$/*$Irrational number (sqrt (5) $-1$ ) $/2*/$$0\mathrm{x}\mathrm{l}7\mathrm{c}\mathrm{e}7301.0\mathrm{x}205\mathrm{C}\mathrm{e}\mathrm{d}\mathrm{c}8,0\mathrm{x}0\mathrm{d}042089$

$\}$

$\{$

unsigned longchar$\mathrm{o}\mathrm{f}_{-}32\mathrm{b}\mathrm{i}\mathrm{t}\mathrm{s}$$\mathrm{o}\mathrm{f}_{-}8\mathrm{b}\mathrm{i}\mathrm{t}\mathrm{S}[4]$;:data-bitarray:
int :

$\}$

$\mathrm{j}$for$(\mathrm{j}=4;\mathrm{j}>=1:)\{$$\mathrm{o}\mathrm{m}\mathrm{e}\mathrm{g}\mathrm{a}[\mathrm{j}]+=$

if (omega
;
else

;
alpha
& CARRY ){ omega$[\mathrm{j}]$$[\mathrm{j}]$$[\mathrm{j}]\iota=$

LIMIT;$\mathrm{o}\mathrm{m}\mathrm{e}\mathrm{g}\mathrm{a}[-\mathrm{j}^{]}++$; }$-\mathrm{j}$

.

$\}$

omega $[0]+=$ alpha $[0]$ omega $[0]t=$ LIMIT;;
omega $[0]\wedge$ omega [11data-bitarray.
data-bitarray.
data-bitarray.
data-byte
data-bitarray.
data-byte;(data-byte $>>4$ )data-byte
data-byte;(data-byte $>>2$ )data-byte)):
return(1 a (( data-byte $>>1$ )$\wedge \mathrm{o}\mathrm{m}\mathrm{e}\mathrm{g}\mathrm{a}[2]$$\mathrm{o}\mathrm{f}_{-^{32\mathrm{b}\mathrm{i}\mathrm{t}}}\mathrm{S}=$$\mathrm{o}\mathrm{f}_{-}8\mathrm{b}\mathrm{i}\mathrm{t}\mathrm{S}[1]$$\mathrm{o}\mathrm{f}_{-}8\mathrm{b}\mathrm{i}\mathrm{t}\mathrm{S}[0]\wedge$$=$$\mathrm{o}\mathrm{f}_{-}8\mathrm{b}\mathrm{i}\mathrm{t}\mathrm{S}[3]$$\wedge \mathrm{d}\mathrm{a}\mathrm{t}\mathrm{a}-\mathrm{b}\mathrm{i}\mathrm{t}m\mathrm{a}\mathrm{y}.\mathrm{o}\mathrm{f}_{-8\mathrm{b}\mathrm{i}\mathrm{t}\mathrm{S}}[2]\wedge$$\wedge$$=$$\wedge$$=$$\wedge \mathrm{d}\mathrm{a}\mathrm{t}\mathrm{a}-\mathrm{b}\mathrm{y}\mathrm{t}\mathrm{e}$

$\}$

void

main $()$

$\{$

int

$\mathrm{j}$;$\mathrm{m}90_{\mathrm{S}\mathrm{e}}\mathrm{t}\mathrm{S}\mathrm{e}\mathrm{e}\mathrm{d}\mathrm{s}$for$(\mathrm{j}=1.

printf$\}$(0.0.0*0*0):\mathrm{j}<=50.

$(^{\mathrm{l}\mathrm{I}}\backslash \mathrm{n}^{\mathrm{l}\mathrm{I}})$:\mathrm{j}++)$

printf$(^{\dagger \mathrm{l}}/.\mathrm{d}\mathrm{l}\iota \mathrm{m}*90\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{o}\mathrm{m}\mathrm{b}\mathrm{i}\mathrm{t} ())$

;

:

86

無理数回転を正確に実行するために, このプログラムでは多倍長加算を行う。 すなわち, 我々が必要としているのは匍 bit だが $(m=90)$ , ここでの加算は 150 bit で行っている。 このため,
少なくとも$2^{50}$個以上の擬似乱数を丸め誤差の影響を受けずに正確に生成することができるであろう 11。
このプログラムでは実行速度を上げるために, 次のようなトリツクを利用している: 関数$f^{(m)}( \omega):=\sum_{i=1}^{m}d_{i(\omega})$(mod 2) の値を計算する部分で, 排他的論理和の演算を用いている。た..

とえば, $\omega\in[0,1)$ の最初の 16 bit が0100111(禾)1011011) $(\omega)=1$ である。次に, この bit の並びであったとしよう。 このとき, 1 が 9 個あるから,を真二つに分けて, それらの排他的論理和 (XOR) をとってみると,$f^{\mathrm{t}^{16}}$

01001110 XOR 01011011$=$

00010101.

となる。演算結果 . は 8bit になるが, これは 1 を 3 個持っているから,$\omega’$$f^{\langle 8)}(\omega’)=1$である。

一般にこの手続きによって 1 の個数の偶奇は変わらないことに注意せよ。上のプログラムではこうしたトリックを何回も用いて計算速度を上げている (演算 XOR ぼきわめて早く処理される)。
もし, アセンブリ言語を使用できる場合はパリティフラグが有用であろう。

52
521

補足
擬似乱数の生成速度について$\mathrm{m}90_{\mathrm{r}\mathrm{a}\mathrm{n}}\mathrm{d}\circ \mathrm{m}\mathrm{b}\mathrm{i}\mathrm{t}$は 1 秒間にパソコン PC-9821Xa10 で約 1,000,000 個のランダムビットを生成する。 この生成速度は従来の擬似乱数生成法より遅いと感じるユーザも少なくないだろう。 しかし計算量の理論によって明らかにされたように (邦語の解説記事としてはたとえば [10]), より精度の高い乱数を生成するためにはより複雑なプログラムがどうしても必要になる。 そのため, 精度の高い乱数生成に時間がかかるのは理論上, 止むを得ない。 最近, 乱数の研究者のあいだでは簡便なプログラムで乱数を生成させることの限界が訴えられるようになってきた (cf. [15])。そのた
め多少時間がかかっても高精度の乱数が望まれている。実際, 乱数生成に時間が多少かかるという欠点はたとえば後述の並列計算によってハードウエア的に克服されるからである。

5.2.2

周期についてよく尋ねられる質問として 「その擬似乱数の周期はどのくらいか」 というのがある。無理数回転を利用した擬似乱数の場合, 理論的には周期は明らかに存在しない (あるいは無限大)。 もっとも, 無理数回転も実際には近似した 「有理数回転」で代用するから, もちろん, 周期はあるわけ
である。
では周期は
で, たとえば関数
は周期の存在しない擬似乱数を実用限界
しかしながら, 前節で述べたように,$\mathrm{m}90\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{o}\mathrm{m}\mathrm{b}\mathrm{i}\mathrm{t}$$2^{150}$$\mathrm{m}90_{\mathrm{r}\mathrm{a}}\mathrm{n}\mathrm{d}\mathrm{o}\mathrm{m}\mathrm{b}\mathrm{i}\mathrm{t}$を超えるほど大量に生成するので機能的には「周期を持たない乱数」を生成していると言える。
を超える周期について議論することはまったく不毛である。なぜなら, それほ
そもそも,ど多くの乱数を使うことは現実にはありえないから。
$2^{100}$

11 前節の 7Table によれば差が大きくなる。 もっとも, 87) $(1\mathrm{o}000)=8.7\mathrm{x}10^{14}<2^{50}$$N_{\mathrm{c}}^{(\mathrm{K}}$

使い切るには 3 カ月以上かかる。

$\mathrm{x}10^{14}$個も生成すると統計的には誤なので擬似乱数をと言う数は実用上十分大きく, 1 秒間にbit を使っても $8.7\cross 10^{14}$ bit を$2^{50}$$10^{8}$

87

5.23

有理数回転による擬似乱数生成

は厳密には「有理数回転による擬似乱数生成」である。 ただ, 余分な桁数関数で計算することにより機能的には無理数回転と同等というに過ぎない。 -方で, 有理数回転による擬似乱数生成自身も十分良い擬似乱数生成法なのである。具体的には
andombit においての関数 $B$ を計算するこで行う。そのとき, 擬似乱数は周期的になるが加算をとによって多次元分布を知ることができて, それは 6Table とほとんど変わらない。加算の負担が減る分だけ擬似乱数の生成速度が少し早くなる。
$\mathrm{m}90_{\mathrm{r}}\mathrm{U}\mathrm{d}\mathrm{o}\mathrm{m}\mathrm{b}\mathrm{i}\mathrm{t}$$\mathrm{m}90r$$4.\mathrm{T}\mathrm{h}\infty \mathrm{r}\mathrm{e}\mathrm{m}$$90\mathrm{b}\mathrm{i}\mathrm{t}$5.2.4

並列計算について無理数回転を利用した擬似乱数生成法は並列計算に大変適している。 いま, プロセヅサが $K-$
個あるとする。最初に擬似乱数の初期値 $\omega\in[0,1)$ を選び, 次にどれか–つのプロセッサで$\omega_{j}:=\{\omega+j\alpha\}$,$j=0,$$\ldots,$$K-1$の無理数回転を行い擬似で無理数を計算する。そして第 番目のプロセヅサには初期値の無理数回転による擬似乱数をで無理数乱数生成を生成させる。 これで, 全体では初期値生成していることになる。注目すべきことに, 初期値の設定後, 計算は各プロセッサごとに完全に独立して行われるので大変能率がよい。
$K\alpha$$j$$\omega_{\mathrm{j}}$$\omega$$\alpha$参考文献[1] P.Billingsley, Convergence of Probabih $ty$ measures, John Willey&Sons, (1968)[2] P.Billingsley, Probability and measure, 2nd edition, John Willey&Sons, (1986)[3] N.Bouleau and D.L\’epingle, Numerical methods for stochastic processes, John Wiley&Sons,(1994)[4] R.Burton and M.Denker, On the Central Limit Theorem for Dynamical Systems,$AMS$ . $108$ -2(1987), 715726$\mathrm{Z}$}$\mathrm{u}ns$.

[5]$\mathrm{G}.\mathrm{J}$.Chaitin, Algorithmic Information Theory, $IBM$ J. Res. Develop. 21 (1977), 350359[6]$\mathrm{J}.\mathrm{N}$.Franklin, Deterministic simulations of random processes, Math. Comp. 17 (1965), 28-

59[7] K.Fukuyama, The central limit theorem for Rademacher system, Proc. Japan Acad. 70,
Ser. , No.7 (1994), 243-246$\mathrm{A}$[8] 伏見正則, 乱数, (東京大学出版会), (1989)[9] 香田徹, “カオスの間接的時系列解析法とその応用”,(1993), 661-668

システム制御情報学会誌, 37, No 11[10] 小林孝次郎, コルモゴロフの複雑さとうンダムネス, bit vo1.27 No 5, 共立出版, May (1995),

88-95.

88[11] T.Kohda and A.Tsuneda, Pseudonoise Sequences by Chaotic Nonlinear Maps and Their
Correlation Properties, IEICE TRans., E76-B, No.8 (1993), 855-862[12] L.Kuipers and H.Niederreiter, Uniform distribution of sequences, Interscience, (1974)$[1.\cdot.3]$

Knuth D.E., The Art of Computer Programming, 2nd ed., Addison-Wesley, (19.81), (邦訳)
準数値算法/乱数 (渋谷政昭訳), サイエンス社, (1983)[14] Martin-L\"of, The definition of random sequences, Inform. Control 7 (1966), 602-619[15] H.Ni\’eerreiter, New developments in umiform pseudorandom number and vector generations, Lecture Notes in Statistics 106, Springer (1995), 87-120.
[16] S.Ogawa, Pseudorandom Functions Whose Asymptotic Distributions Are AsymptoticallyGaussian, Math. Anal. and Appl., 158, No.1, (1991)[17] S.K.Park and.Miler (訳:西村恕彦), “乱数生成系で良質なものはほとんどない”, bit(共立出版) 4 月号, 5 月号, (1993)$\mathrm{K}.\mathrm{W}$[18] H.Sugita, Pseudo-random number generator by means of irrational rotation, Monte Carlo
Methods and Applications, VSP, 1-1, 3557(1995).
[19] 数理解析研究所講究録 498, 乱数プログラムパヅケージ, (1983)[201 数理解析研究所講究録 850, 確率数値解析における諸問題, (1993)[21] 杉田洋, 無理数回転による擬似乱数生成, 数理解析研究所講究録 915, 数値計算アルゴリズ
ムの現状と展望 II, (1995)[22] H.Sugita and S.Takanobu, Limit theorem for symmetric statistics with respect to Weyltransformation: Disappearance of dependency, (1996), preprint[23].Tausworthe, Random numbers generat\’e by linear recurrence modulo two, Math.
Comp. 19, (1965), 201-209$\mathrm{R}.\mathrm{C}$[24] 津田孝夫, モンテカルロ法とシミュレーション, 培下館, 改定版 (1977)[25] P.Walter, An Introduction to Ergodic Theory, Springer, (1981)杉田四九州大学大学院数理学研究科 (工学部分室)
812-81 福岡市東区箱崎 6-10-1
E–mail: sugita@math kyushu-u .jp$.\mathrm{a}\mathrm{c}$