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

柴田章博、松古栄夫

2006年6月13日-15日

はじめに

「数値シミュレーション入門」は統計力学に基づくモンテカルロシミュレーションの手法の基礎を紹介する。磁性体の統計力学のモデルは、本来の磁性体の性質を解明のみならず、さまざまな広範な問題に応用されている。--物性物理学やのや物理化学のシミュレーションや素粒子物理学の格子ゲージ理論シミュレーションののみならず、統計的手法に基づく情報処理の分野やニューラルネットワークの学習理論などさまざまである。--

本実習では、もっとも基本的なモデルあるイジングモデル（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

ベイズの定理

２つの事象$A,B$があるとき、同時に起きる確率$\Pr (A,B)$と$B$がおきた条件のもとで$A$がおきる条件付確率$\Pr (A|B)$および、そ逆$A$がおきたもとで$B$がおきる条件付確率$\Pr (B|A)$のあいだには、次の関係がある。MATH したがって次がえられる。MATH

これは、ベイズの定理（ベイズの公式）と呼ばれる。

連続分布関数の例

連続型一様分布

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

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

正規分布

正規分布(nomal distribuytion)は次のように与えられ、$N(\mu ,\sigma^{2})$で略記する。MATH 特に、$N(0,1)$のものを標準正規分布とよび、この分布に従う確率変数は、 しばしば$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 は相関関数と呼ばれる。この量は、しばし$Y_{1},Y_{2}$との間の線形相関関係の強さを示す量と解釈され、$Y_{1},Y_{2}$とが独立ならば、必ず$\rho =0$となる。しかし、逆に$\rho =0$であっても必ずしも、$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はMATHに、MATHはMATHそれぞれ従う。

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

前小節で行った、確率変数の和に対して次の定理がえられる。

大数の弱法則

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

中心極限定理

$X_{k},k=1,...,n$が互いに独立かつ同一の確率分布に従う確率変数列で、MATH とする。このとき、MATHにたいして、つぎが成立する。MATHすなわち、$S_{n}$は常に、正規分布に近づいていくことが知られている。これを中心極限定理と呼ぶ。

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

統計物理学に基づくモンテカルロシミュレーションに入る前に、モンテカルロ法 Note_1 に基づく数値積分について概観する。モンテカルロ法は、端的にいえば「サイコロを振って積分する方法」ということができる。

一様サンプリング(uniform sampling)

はじめに、もっとも簡単な例を考える。積分領域$\Omega $と、関数値域$f(x),(x\in I)$を囲む領域$A$を考え、領域$A$にランダムな点を振る。曲線$f(x)$の下側の領域の割合と、領域$A$の体積（面積）との積で積分が推定することが出来る。MATH この方法は、曲面（曲線）$f(x)$が複雑な形状をしていても、曲線$f(x)$の下側の領域の点の生成と積分を容易に行うことが出来る利点がある。しかし、一般に効率のよい方法ではなく、見積もりの誤差も議論されなければならない。

一様サンプリングによる、モンテカルロ積分の推定法を天下り的にではあるが導入する。多次元空間の領域$\Omega $(体積$V)$の一様なランダムな点$N$個の点MATHをサンプリングしたとする。モンテカルロの基本定理によれば次の式で推定される。MATH ここで、角括弧は、$N$個の点における算術平均を意味する。すなわち、MATH 式（eq:MonteI）の$\pm$の項は、標準偏差に相当する誤差の推定値であり、厳密な限界を示す誤差ではない。 Note_2 また、誤差がガウス分布である保障もない。大雑把に見積もった値である。

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

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

一様サンプリングの方法は、積分の誤差は 一様サンプリングされた関数の分散を用いて、$\sqrt{V[f(x)]/N}$で与えられた。これは、非積分関数の値が大きく変化をする場合には、推定誤差が大きくなってしまうことを意味している。このような状況を改善するために積分を次のように書き換えることを考える。MATHとする。いま、一様な重み $dv$でサンプリングする代わりに、重み$\rho (x)dv$でサンプリングができるとしたらどうなるであろうか。

このようサンプルされた点は、確率分布$\rho (x)$で点が生成することに対応する。確率分布であるので$\rho (x)>0$であって、MATH を満たす関数である。いま、式（eq:EV-Ensamble）によれば、MATH であるから、積分MATHの値とその誤差は、$\rho (x)$分布に従うサンプリング点$Y$で関数の値$h(y)$の期待値と分散で与えられる。MATH 一様サンプリングの場合は、$\rho =1/V=$一定の場合である。

この表式は$\rho(x)$の選び方によらない。一様サンプリングの方法を改善するために、もっとも効率のよい$\rho (x)$の選び方を求めよう。そのような$\rho $は分散が与えられた$f$にたいして最小になるように選べばよい。すなわち、拘束条件（eq:ProbCond）をみたし、MATH が最小になるような$\rho $を求めれよい。そのような$\rho (x)$はラグランジュの未定乗数つきの変分MATH で与えられる。これを解くと$\rho (x)$は次の最適解が得られる。MATH

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

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

統計力学と情報処理

Ising model

多数の要素が互いに影響を及ぼしあっている（相互作用している）とき、どのようなマクロな性質を持つかを解明するのが統計力学の重要な課題である。たとえば、水（$H_{2}O$）が、氷、水、水蒸気を温度や圧力などの環境でマクロな性質が大きく異なる状態を取る。このような急激な変化（相転移）を理解する（水の状態を記述する理論ではないが）最も基本的なモデルとして、Isingモデルを導入する。

$1,...,N$までの番号づけれれた格子点(site)に,値を持つ場を考えよう。格子点は、結晶の原子の配置でもよいし、（白黒）画像のピクセルやニューロン（神経細胞）でもよい。ここでは、$S=\pm 1$の２値を持つ、Ising スピン（spin）を対応させたモデルを考えることにする。磁性体のモデルでは、$S=\pm 1$は、原子レベルの磁石（のNSの向き）が、上下のいずれを向いているかをあらわしている。いま、お互いに隣り合ったサイト同士$<i,j>$が相互作用するとして、相互作用のエネルギーとして、$-JS_{i}S_{j}$を考える。$J>0$の時、$S_{i}S_{j}>0$のとき、スピンが揃うときがエネルギーが低くなり安定である。一方、$J<0$のときは、$S_{i}S_{j}<0$のとき、スピンがお互い反対向きのときが、安定となる。また、磁性体のときにおける外部磁場$h$との相互作用 $-hS_{j}$を付け加えた系を考える。統計力学では、系全体のエネルギー取り扱うハミルトンの関数（ハミルトニアン）を用いて系を記述する。MATH Isingモデルにおいては、系の相互作用は最近接相互作用でその強さが共通であるが、相互作用の形がどのようになるかは問題によって異なることを注意しておく。

統計力学の処方箋（詳しくは、統計力学とIisngモデルで解説）によると、ハミルトニアンが与えられると、スピンがとる配位の確率分布は、MATH で与えられる。ここで$T\,$は温度、$k_{B}$はボルツマン(Boltzmann) の定数である。

画像修復の問題

Isingモデルの統計力学の解析の例として画像修復の問題を紹介する。

雑音によって乱れた不完全な画像が与えられたとき現画像を推定する問題は、ベイズの公式を用いて、確率論を用いた画像修復の問題として、統計力学の手法を応用することができる。２次元の白黒の画像を考えると、各ピクセルは、格子点上のイジングスピン$\{\xi _{k}\}$と対応付けることができる。元の画像$\{\xi _{k}\}$は、受信側では$\{\tau_{k}\}$であらわすこととする。画像は、伝送の間に乱され確率 $p$で、ビットが反転すると、$\{\xi _{k}\}$が与えられたときの$\{\tau _{k}\}$の出現確率は次で与えられる。MATH 各サイトでビット反転が独立であるとすると、全体では、MATH である。画像修復では、劣化画像MATHから、原画MATHを推定するもんだいであるので、推定画像を$\{\sigma _{k}\}$と書くとすると、ベイズの公式を用いて次のようにかける。MATH これは、劣化画像MATH及び、事前知識MATHが画像修復に必要であることを示している。しかし、MATHが正確にわからなくても、モデル分布MATHをMATHの代わりに使うこととする。

ここでは、モデル分布MATHを次のように考えて構成する：白の背景の中に孤立した黒の点がある場合は、もともとあった黒であるよりも、白の点がノイズによって反転した可能性が大きいと考える。このことを、隣り合うピクセル対は、違う状態であるよりもおなっ状態である確率が高いと解釈し、次の事前分布のモデルを導入する。MATH ここで、$<i,j>$は隣同士のピクセルについての和であり、$Z(\beta _{m})$は企画化因子である。これはイジングモデルをあらわす。したがって、式（eq:BaysePix）へ事前分布を代入すると、画像修復の確率分布（事後分布）は、Ising模型にランダムな磁場がかかったスピン系の統計力学の問題となった。 MATH 最大事後確率の考え方にたてば、式（eq:MAPPiX）を最大にするMATHが修復画像となる。

  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