総研大・KEK夏期実習:数値シミュレーション入門 --Ising Model のMonte Calro シミュレーション --

柴田章博、松古栄夫

2008年6月2日-4日

はじめに

「数値シミュレーション入門」は、確率的計算技法としてとしてモンテカルロ法に基づくシミュレーションの入門の実習を行う。モンテカルロ法は物理学においてはメトロポーリスらによって流体シミュレーションの方法として提唱され、統計物理学など様々な分野応用されるようになった。統計力学に基づくモンテカルロシミュレーションの手法の基礎を紹介する。磁性体の統計力学のモデルは、本来の磁性体の性質を解明のみならず、さまざまな広範な問題に応用されている。物性物理学やのや物理化学のシミュレーションや素粒子物理学の格子ゲージ理論シミュレーションののみならず、統計的手法に基づく情報処理の分野やニューラルネットワークの学習理論などさまざまである。

本実習は、モンテカルロ積分の基本的なアイデアから出発して、統計力学の基本的なモデルあるイジングモデル（Ising model）を題材にその統計物理学的な解析的法（解析的び数値的方法（数値シミュレーションの方法））を解説する。また、実際にシミュレーションをおこないさまざまな量を測ることによって数値シミュレーションの手法を体得することを目的とする。

確率分布

分布関数・確率密度

ある変数のとる値が事前に予知されず確率$p$にしたがって出現するとみなすとき、これを確率変数といい $X$,$Y$,$Z$などの大文字で表現する。確率変数には、サイコロのように有限個の整数などの値をとる離散型と身長、体重などのように取り得る値の範囲を持った連続型とがある。

確率分布は、ある変数がどの値をとる可能性が大きいか、小さいかを実現値と確率の間の関数として現したものである（通常 $x,y,z$などの小文字で表現する）。ある変数が $x$以下である確率を知りたいとき、MATH と現す。この関数はある$x$を与えるとこめられることであり、累積分布関数（単に分布関数）と呼ばれる。この関数は単調非減少関数で、有限の値を仮定して次のようのにとる。MATH 分布関数は、微分可能な場合が多く、その導関数を確率密度関数と呼び、$f(x)$と書くことにする。$f(x)$はすべての$x$に対して、$f(x)\geq0$であり、MATH と現される。特に、MATHである。また、$a\leq X\leq b$である確率は、MATH とあらわせる。連続変数の場合、一点の積分はゼロであるので、$\Pr(X=x)=0$である。

一方、離散変数である場合には、$p_{i}=\Pr(X=x_{i})$であって、分布関数は積分記号でなく、和記号で現される。MATH

確率変数の分布の様子を示す指標として、分布の中心やバラつきの度合いを示す指標があると便利である。分布の中心を現すものとして、代表的なのが期待値と分散である。MATHたたし、MATH一方、連続変数に対して次で与えられる。MATH

連続分布関数の例

連続型一様分布

連続型確率変数$Y$が区間$[a,b]$間で確率密度一定である分布を連続型一様分布（Unifoem distribution）といい、$Y\sim U(a,b)$と略記する。これは、確率密度MATH で与えられる。$p(y)\geq0$、かつMATHである。

一様分布の平均、分散は次で与えられる。MATH

正規分布

正規分布(nomal distribuytion)は次のように与えられ、$N(\mu,\sigma^{2})$で略記する。MATH 特に、$N(0,1)$のものを標準正規分布(図normDist)とよび、この分布に従う確率変数は、 しばしば$Z$と関数$\Phi(z)$を用いて現される。MATHと $Z\sim N(0,1)$は、MATHの関係でいつでも変換することができる。正規分布の平均、分散は次で与えられる。MATH

多変数の確率分布（2次元）

連続型確率変数の組MATHについて分布を考える。それぞれの区間$(a_{1},b_{1}]$、$(a_{2},b_{2}]$に変数が同時に入る同時確率は、同時確率密度関数 $\ p(y_{1},y_{2})$をもちいてMATH で与えられる。MATH を満たさなければならない。$Y_{1,}Y_{2}$の周辺分布は次のように得られる。MATH このように、積分によって変数を消去することを周辺化という。

$Y_{1,}Y_{2}$の同時確率密度（同時分布）$p(y_{1},y_{2})$が与えられたときMATHMATHとするときMATHは,$Y_{1},Y_{2}$の共分散とよぶ。また、MATHは相関関数と呼ばれる(MATH)。この量は、しばし$Y_{1},Y_{2}$との間の線形相関関係の強さを示す量と解釈される。相関係数が最大（最小）の$\pm 1$をとる場合は、MATHときである。言い換えると、すべてのデータが直線(eq:line)上にあるときである。

一方、$Y_{1},Y_{2}$とが独立ならば、必ず$\rho =0$となる。実際、 MATHの時、MATHである。しかし、逆に$\rho =0$であっても必ずしも、$Y_{1},Y_{2}$は独立とは限らない。たとえば、半径１の円周上に点が分布する場合MATH（ここでMATHはディラックのデルタ関数）には、MATH,MATHであるが、$f(y_{1},y_{2})$は独立分布ではない。

独立な確率変数の和

統計処理においては、ランダムな誤差を含んだ測定値や観測量の集計を行う。その基本的な、操作である確率の変数の和$X+Y$についてその性質を調べる。期待値については、加法性が成立する。MATHさらに$X,Y$が独立であるときは、分散の加法性が成立する。MATH$X,Y$が独立でない時は、MATHである。分散の加法性が成立するためには、無相関であれば十分となる。

ｎ個の変数MATHにたいして期待値の加法性が成立する。MATH また、MATHが独立な場合には、分散の加法性が成立する。MATH

同一分布

いま、ｎ個の変数MATHが同一分布に従うとき、これらの期待値、分散を$\mu ,\sigma ^{2}$であるとすれば、$E[X_{k}]=\mu $,MATH,$\quad (k=1,..,n)$であるのでMATHである。したがって標準偏差は、$\sqrt{n}\sigma $と$\sqrt{n}$に比例することが重要である。MATHの相加平均MATHとおくと相加平均の期待値、分散は次で与えられる。MATH

和の分布関数

確率変数 $X,Y$が独立で、それぞれ$f(x),g(y)$に従うとすると、和 $Z=X+Y$の従う確率分布$h(z)$を考える。$X+Y=z$となるのは、$X=x,$$Y=z-x$の同時分布のすべての和であるから、離散分布のときMATH であり、連続分布のときは、和を積分に置き換えればよくMATH となる。関数$f(x),g(y)$から、$h(z)$をつくる操作を畳み込みという。

たとえば、$X,Y$が正規分布MATH、MATHに従うとき、その畳み込み（記号 $\ast$を使って表現する）は、MATH となる。すなわち、MATHの正規分布になる。

多変数$X_{k}$がそれぞれMATH（$k=1,..,n$）であるとき、MATHはMATHに従う。とくに$X_{k}(k=1,...,n)$が同一分布MATHに従うとき、MATHにそれぞれ従う。

大数の法則と中心極限定理

大数の法則

いま、$X_{k},k=1,...,n$が互いに無相関で期待値$\mu $、分散$\sigma ^{2}$であるとすると、$\bar{X}=$ MATHにたいして、MATHが成り立つ。このときチェビシェフの不等式MATHを適用すると、MATHに対してMATHが成立する。ここで、$\delta $がどんなに小さな正定数であっても、十分大きな$n$をとることで、式（eq:LargeNth）の右辺はいくらでも小さくすることができる。すなわち、$n$を大きくとることで、$\bar{X}$が$\mu $の近辺に集中していることが２次の積率（moment）までの確率でいえる。この事実を大数の弱法則（大数の法則）という。標語的にいえば、「大標本では観察された標本平均を母集団の真の平均とみなしてよい。」となる。

中心極限定理

$X_{k}$,$\quad k=1,...,n$が互いに独立かつ同一の確率分布に従う確率変数列で、MATH とする。このとき、MATHにたいして、つぎが成立する。MATHすなわち、$S_{n}$は常に正規分布に近づいていくことが知られている。（証明略）。これを中心極限定理と呼ぶ。標語的にいえば「和MATHは$n$が十分大きければ、大体正規分布であると考えてよい。」となる。

モンテカルロ積分(Monte Calro integration)

モンテカルロ法 Note_1 というと乱数を使った円周率$\pi$の計算の例を聞いたことがあるかもしれない。統計物理学に基づくモンテカルロシミュレーションは、必ずしも単純な積分を意味しないが、モンテカルロ技法の一端を見ることができるるので、モンテカルロ積分について概観する。

実際の円周率の計算には、計算誤差が大きくモンテカルロ法が用いられることはないが。しかし高次元の数値積分においては、計算手続きが指数的に増大し、積分区間を$M$点に分割する$d$-次元の積分公式を使うと計算量が$O(M^{d})$のオーダになる。非常に自由度の大きな積分を実行するには、一般に計算量が積分の次元に依存しないモンテカルロ法が有利となる。

領域サンプリング

最も簡単な例を考える。積分領域$\Omega $と、関数値域$f(x),(x\in I)$を囲む領域$A$を考え、領域$A$にランダムな点を振る。MATHが分かっているとすれば、曲線$f(x)$の下側の領域の体積は、生成したランダムな点が曲線$f(x)$の下側の領域に含まれる割合として計算できる。すなわち、$\int_{A}f(x)dv$が次のように計算される。

MATH この方法は、曲面（曲線）$f(x)$が複雑な形状をしていても、曲線$f(x)$の下側の領域の点の生成と積分を容易に行うことが出来る利点がある。しかし、一般に効率のよい方法ではなく、見積もりの誤差も議論されなければならない。

一様サンプリング（Uniform Sampling）法

$d-$次元空間の領域領域$\Omega$(体積$V)$の変数MATHの関数 $f(\QTR{bf}{x})$の多重積分MATH を求めることを考える。モンテカルロ積分（一様サンプリング法）では、積分領域$\Omega$から一様ランダムに抽出した$N$個の点MATHを使って積分値を推定する。（さいころを振って積分する。）代表点$X$の大きな$N$の極限で近似値MATH は、大数の法則によって、積分MATHに収束する。ここで、角括弧は、$N$個の点における算術平均を意味する。一様ランダムサンプリング（同一分布）に対する、相加平均の期待値、分散は式（eq:EV-Ensamble）で与えられるから、$\tilde{I}$の見積もりの誤差はMATH となる。台形公式における誤差オーダーは $O(M^{-2})$であるのに対し、モンテカルロ法の誤差の見積もりは $O(M^{-1/2})$であり、とても悪く思われる。しかしながら、積分区間を$M$点に分割する$d$-次元の積分公式を使うと計算量が$O(M^{d})$のオーダであるため、積分公式に基づく積分値の推定には高次元ではべき級数的に増大する。一方で、モンテカルロ法では$O(M)$の計算量となる。

実際の計算を行う際には、与えられた領域$\Omega$の一様サンプリングされた点と体積$V$がわかっていることが前提となるが、実際の計算できるのは$n$次元直方体の領域など限られた場合である。$\Omega$が複雑な形状を持った領域の場合は、体積の見積もり及び一様サンプリングか容易に行える$\Omega$を囲む領域$\Omega^{\prime}$を導入することで、$\Omega$内の点を生成することが可能である。また、領域$\Omega^{\prime}$(体積$W$)で定義された関数$g$を$\Omega$内の点$x$に対して、$f(x)$、外の点においてゼロであるように定義すればMATH として積分を推定することができる。

重点サンプリング（important sampling）

一様サンプリングによるモンテカルロ積分は次のような問題点がある。高次元空間の中の一様のサンプリングを生成することは難しい。一様にサンプリングポイントが生成できたとしても、被積分関数の評価値がほとんどの点で小さく結果的に無駄な計算となり、評価誤差が大きくなってしまう。

このような状況を改善するために積分を次のようにある確率分布関数 $\rho(x)$を用いて書き換えることを考える Note_2 ：MATH ここで、$h(x)=f(x)/\rho(x)$とおいた。 いま一様な重み $dv$でサンプリングする代わりに、重み$\rho(x)dv$でサンプリングができるとしたらどうなるであろうか。このようサンプルされた点は、確率分布$\rho(x)$で点が生成することに対応する。

相加平均の期待値、分散の式（eq:EV-Ensamble）によれば、MATH であるから、積分MATHの値とその誤差は、関数$\rho(x)$に従うサンプリング点$Y$で関数の値$h(y)$の期待値と分散で与えられる。MATH 一様サンプリングの場合は、$\rho=1/V=($一定$)$ の場合である。

この式は$\rho(x)$の選び方によらないことに注意。したがって、一様サンプリング法を改善するためには、$\rho(x)$をうまく選んで積分に寄与するところを重点的にサンプルすればよいように思われる。

では、もっとも効率のよい$\rho(x)$の選び方は何であろうか。それは、積分の推定誤差$\rho$が与えられた$f$にたいして最小となるように選べばよい。すなわち、拘束条件（eq:ProbCond）をみたし、MATH が最小になるような$\rho $を求る変分問題MATH をとけばよい。これを解くと$\rho(x)$は次の最適解MATH を満たせば良い。

実際には、重点サンプリングでモンテカルロ積分を行うには、$\rho(x)$の分布を与えるアルゴリズムを構築しなければならないが、一般の関数MATHに対して構成することはできない。この節では、一般の$f(x)$に対する重点サンプリングのアルゴリズムについては触れることはしない。重点サンプリング法を拡張した、適応型モンテカルロ積分の方法として、VEGASアルゴリズムが知られている。（たとえば、Numerical Recipies in C/C++ NRinC/C++ を参照）。

本実習で扱う統計力学に基づくのモンテカルロシミュレーションが、多次元の重点サンプリングの分布関数を生成するアルゴリズムを与えるあたること、さまざまな物理量を測ることが重点サンプリングによる多次元積分を実行することに対応することをみる。

  1. Donald E. Knuth, "The art of Computer Programing" Vo2. Seminumerical Alogorithm. Addidon-Wesly, 2nd edition 1981. (サイエンス社からの邦訳あり)

  2. 奥村晴彦 「C言語によるアルゴリズム事典」 技術評論社ISBM4-87408-414-1

  3. 伏見正則「乱数」東京大学出版,1989

  4. G.E.P. Box amd M.E. Muller, An. math. Statist. 26: 610 (1958); G.Mardagla and T.A.Bray, Rev. Soc. Ind. Appl. Math. 6 260(1964)

  5. Numerical Recipies in C/C++ the art of Scientific Computing 2nd edition, W.H. Press 他著、Cambridge Press

  6. 統計的情報処理と統計力学, 田中和之編著, 臨時別冊・数理科学 SCGライブラリ50