\documentclass[12pt]{jarticle}
\usepackage[dvips]{graphics}

\oddsidemargin 0mm
\evensidemargin 0mm
\textwidth 165mm
\topmargin -15mm
\textheight 230mm

\author{}
\title{計算物理学 No.7(モンテカルロ法)}
\date{}

\begin{document}
\maketitle

統計力学など、さまざまな物理現象(他に集団遺伝学、原子炉中の中性子拡散な
ど)を表すのに広く使われるのがモンテカルロ法
である。高次元の積分など他に有効な数値計算法が無い場合にも有効である。

\section{疑似乱数}

まず、モンテカルロ法を使う前に乱数を発生させる数値的方法を考察する。
厳密にはプログラムで生成される数列は決定論的なものであるから疑似乱数と呼
ぶべきである。

シミュレーションに用いる疑似乱数列は以下の性質を持つことが要求される。
\begin{enumerate}
 \item 統計的性質がよい。
 \item 計算効率が高い。
 \item 周期が長い。
 \item 再現性がある。
\end{enumerate}

\subsection{線形合同法}

漸化式
\begin{equation}
 I_{j+1} = {\rm mod} ( a I_{j} +c , M)
\end{equation}
により発生される数列 $I_{1},I_{2},I_{3},\ldots, (0 \le I_{j} < M )$ を疑似乱数列とするものを
線形合同法(1948年頃提案された)と呼ぶ。
ここで、$M$ は正の整数で {\bf 法} と呼び、$a,c$ は $0 \le a,c < M$ とい
う整数でそれぞれ {\bf 乗数}、{\bf 増分}と呼ぶ。
$ I_{j}$ を $M$ で割る
$x=I_{j}/M$
と、[0,1) 区間の一様乱数
$0\le x <1$ を得ることが出来る。

線形合同法により生成される疑似乱数列はいつかは元に戻り、その周期は $M$
以下である。パラメータ $a,c$ の選び方により $I_{j}$ の性質は大きく異なる。
例えば、$a=0$ では $I_{j}=c$、
$a=1,c=0$ では $I_{j}=I_{1}$、
$a=1,c=1$ では周期 $M$ ではあるが、とうてい乱数列とは呼べない。

[課題1]
線形合同法による疑似乱数列を生成するプログラムを作成してみよ。

\subsubsection{パラメータの推奨値}

パラメータ $a,c$ の選び方には以下のようなもの
が推奨されている(文献\cite{Okumura} より引用 )。

「
 M が 2 の累乗なら ${\rm mod} (a, 8)$ を 5 または 3 とし（5 の方が安全）、 定数項 c
 を奇数とする。このとき、周期はちょうど M となり、 その 1 周期分には 0
 から M - 1 までの整数が 1 個ずつ現れる。

高速にするには c = 0 とし、初期値 $I_{1}$ を奇数にする。 ただし、周期は M / 4
になる。 
」

\subsubsection{
線形合同法の欠点
}
線形合同法は簡便かつ有効な疑似乱数発生法ではあるが、いくつか欠点がある。

\begin{enumerate}
 \item 合同法乱数は、上位の桁はランダムだが、下位の桁はランダムでない。 
 \item 合同法乱数は、1 個ずつ使えばランダムだが、いくつか組にして使えば
 ランダムでない。 
\end{enumerate}
性質 1. から説明すると、線形合同法による疑似乱数列は
下位 k bit だけを見ると周期は高々 $2^k$ である。

\medskip

[課題2] 
[課題1] で作ったプログラムに対し、上記の性質を確認してみよ。
複数のパラメータに対し試してみること。

\bigskip

性質 2. に関しては
たとえば 0 から 15 までの乱数で、(x, y) の座標を生成することを考える。
ここで周期 16 の線形合同法で乱数を生成してみる。すると、 生成され
る (x, y) の座標は 16 / 2 = 8 通りにしかならない。 (x, y) の座標は
$16^2=256$ 個もあるのに、8 個の座標しか生成できないのです。 (x, y, z) の
場合は $16^3=4096$ 個の可能な座標があるのに、16 個の座標しか生成できません。 

\medskip

[課題3] 
[課題1] で作ったプログラムに対し、上記の性質を確認してみよ。
複数のパラメータに対し試してみること。
gnulpot の機能を使うと、2次元プロット、3次元プロット（splot を使う)で可
視化し易くなる。

\subsection{M 系列}

線形合同法は簡便かつ有効な疑似乱数発生法ではあるが、「1 個ずつ使えばラン
ダムだが、いくつか組にして使えばランダムでない」 などの欠点がある。

線形合同法のほかには、M 系列乱数 (M-sequence random numbers) というアル
ゴリズムが有名である。詳しくは述べないが、参考文献として
\cite{M-sequence-random-numbers} をあげておく。

\subsection{Mersenne Twister (MT)}

この他、1996-7 年に提案された、Mersenne Twister (MT) と呼ばれる疑似乱数
発生法もある。参考文献 \cite{MT} 参照 。

\section{モンテカルロ法}

乱数を使った無作為抽出を利用して数学や物理の問題を解こうという手法。

例： 円周率 $\pi$ の計算

正方領域
$0\le x <1, 0 \le y <1 $
内の点を乱数で多数発生させ、 1/4 円内
\[
 y \le \sqrt{1-x^2}
\]

に入る点の割合を調べる。
この割合は $1:\pi/4$ になる。

\begin{thebibliography}{99}
\bibitem{Okumura}
	「C 言語による最新アルゴリズム事典」 奥村 晴彦 著、技術評論社 (1991))
\bibitem{M-sequence-random-numbers}
http://www.ysr.net.it-chiba.ac.jp/data/rand/node6.html

http://www.edu.cs.kobe-u.ac.jp/Algorithm/pdf/Note03.pdf 
\bibitem{MT}
http://www.math.keio.ac.jp/~matumoto/mt.html

http://www.emit.jp/mt/mt.html
\end{thebibliography}

\end{document}
