数値シミュレーション入門 --計算機実習--

柴田章博、松古栄夫

計算機での分布関数の生成

一様乱数の発生：線形合同法

線形合同法は、整数の一様乱数を生成するポピュラーな方法で、適当な初期値 $I_{o}$を出発点とする次の漸化式で与えられるbook:Okumura MATH この漸化式で、$0\leq I_{k}<M$の整数を次々と生成することができる。 ここで演算$\quad p\func{mod}q$は 整数$p$を整数$q$で割ったときの剰余をあらわす。$M$が、2の累乗であるならば、 $a\func{mod}8$を5ないし3（5が安全）として、定数項は奇数とする。このとき周期がちょうど$M$である、その周期に $0$から$M-1$の整数が1つづつ現れることが知られている。ref:ICM

線形合同法の 下位$k$ビットを見ている乱数の周期は高々$2^{k\text{\quad }}$ビットであり、線形合同法ではMATH ことに注意する必要がある。$M=2^{32}$のとき、$a$のとりかたとして、MATH などがある。Cの unsigned long の算術で $M=2^{32}$の線形合同法の演算は正しい下位ビットを与えることが保障されている。高速にするには、$c=0$と選ぶ方法がある。

しかしながら、「高級言語」での実装は一般には簡単ではない。なぜなら、積をとると$M$よりも大きな数となり桁あふれを起こし、正しい下位ビットを与えないからである。この問題を回避する方法として、Schrage'sのアルゴリズムが知られている。MATH と $M$を分解して計算をする。ここで、記号$[p]$は、$p$を超えない最大の整数、$r$は$M$を$a$で割った剰余である。もし$r$ ($r<q)$ が小さく、$z$が$0<z<M-1$満たす整数であるとすれば、$a(z\func{mod}q)$と$r[a/q]$はともに、$0$と$M-1$の間の数とすることができる。まとめて書くと次のように与えられる。MATH このようなポータブルな生成ルーチンのパラメータは次のものが知られている。NRinC/C++

$a$ $q$ $r$ $M$
$16807$ $127773$ $2836$ $2^{31}-1$
$48271$ $44488$ $3399$ $2^{31}-1$
$69621$ $30845$ $23902$ $2^{31}-1$

さらに、乱数の切り混ぜを行って周期の長い乱数を生成する方法がある。

正規分布の生成

簡便な方法

最も簡単な正規分布の生成法は、0と１の間の一様乱数 Uを12個をたして6を引く方法がある。これは、一様乱数に限らず、ランダムな数を多数加えると正規分布に近づくという中心極限定理を利用したものである。一様乱数$U$は平均$1/2$ 分散$1/12$ となるので、12個をたして6を引くとちょうど平均０分散１の乱数になる。

ボックス・ミュラー法

よりよい正規分布の生成法として、ボックス・ミュラー法が知られている。BoxMuller このアルゴリズムの概要は次のとおり。

それぞれ独立な平均ゼロ、分散１の正規分布を $x,y$とすると、同時分布はMATH であたえられる。これを局座標で表示すると、MATH であり、同時分布は$\exp (-s/2)ds$に比例する。密度関数が、$\exp (-s/2)$ （$s\geq0$）に比例する分布関数は、$s$で積分してつぎであたえられる。MATH MATH as $s\rightarrow\infty$ であるような、乱数$s$ を作るには、$0\leq U<1$の一様乱数を分布関数の逆関数で変換して、MATH として得ることができる。したがって、これから $x,y$をえるには、$0\leq\theta<2\pi$ の一様乱数をもとめ、MATH とする。 または、MATHの一様乱数$\xi,\eta$を生成し、MATHをみたす$\xi,\eta$を選ぶ。この$u$から$s=F^{-1}(u)$を計算し、$x=s\xi/\sqrt{u},$ $y=s\eta /\sqrt{u}$をとると、正規分布が2つ同時に得られる。

演習：確率・統計/モンテカルロ積分

  1. サンプルプログラムを参考に一様乱数$U(0,1)$を生成するプログラムを作成せよ。一様乱数を発生しヒストグラムを作成し、一様乱数が発生をしていることを確かめよ。

  2. 正規分布 $N(0,1)$を発生させるプログラムを作成せよ。

  3. 2種類の方法を用いて、円周率$\pi $を計算せよ。それぞれの場合において、サンプリング数を変えて真の値との差をプロットせよ。

  4. [発展 ]区間$[0,1]$を$M(=20)$分割して、各区間の中点$x$で$\sqrt{1-x^{2}}$となるようなヒストグラムで与えられる分布関数を$\rho (x)$とする。

演習：統計力学に基づくシミュレーション

一次元イジングモデルの数値シミュレーション（J>0）

サンプルプログラムを参考に、次の手順によって、２次元$(N\times N)$ の格子上にあるイジングモデルのシミュレーションを実行せよ。

  1. 熱浴法(heat bath)によって配位を生成するプログラムを作成せよ。

  2. 配位を、cold stat で生成し、配位の変化する様子を可視化して調べよ。

  3. Metropolis 法によって配位を更新するプログラムを作成し、配位の更新の様子を可視化して調べよ。

  4. 各モンテカルロスウィープ（sweeep）において、スピンの時空間平均MATH を計算し、その変化をプロットせよ。

  5. すべてのスピンが上向きでそろった状態（cold start）から,さまざまなMATH) でシミュレーションを実行し磁化$\bar {S}$を計算せよ。また、$\beta $と磁化をの関係をグラフしなさい。

  6. それぞれの$T$において、磁化を計算せよ

  7. それぞれの$T$において、比熱を計算せよ。

  8. それぞれの$T$において、磁化率を計算せよ。