《講義ノート》

物性研究・電子版 Vol. 7, No. 2, 072214 （2018年11月号）

モンテカルロ法の基礎と応用— 計算物理学からデータ駆動科学へ —東京大学 大学院総合文化研究科 広域科学専攻 相関基礎科学系

福島孝治 1

1953 年に Metropolis らによって平衡統計力学におけるカノニカル平均の数値計算方法としてモンテカルロ法は提案された [1]．この方法はマルコフ連鎖の性質を利用することから，一般の確率分布からのサンプリング方法とみなすことができ，その後，物理学の様々な問題に応用されるようになる．特に相転移の研究では，モデルに含まれるミクロな自由度を直接サンプリングする非摂動論的な方法として，相転移描像の解明に重要な役割を果たしてきた．1980 年代になるとベイズ統計の実践的な計算方法としての利用が見出され，統計学でも方法の発展と応用が精力的に行われるようになった．最近では，このベイズ統計を用いて自然科学の実験・計測データ解析を行なう研究が注目されてきている．そこでは計測データを入力として，それを説明しうるモデルに含まれるパラメータがサンプリングの対象となる．このアプローチでモンテカルロ法は理論と実験の中間の位置で仮説検証ループの潤滑油のような役割が期待されている．このゼミではマルコ
フ連鎖モンテカルロ法の基礎から始めて，相転移研究やデータ駆動科学とも呼ぶべきデータ解析への応用例を示しながら，モンテカルロ法でできることやわかることを議論してみたい．

乱数を用いた計算技法としてのモンテカルロ法

1

モンテカルロ法とは乱数を用いた数値計算技法の総称であり，歴史的には 18 世紀の Buﬀon の針と呼ばれる実験に始まると言われている．床に引かれた多数の平行線に針を落し，交差頻度から円周率の近似値を評価した針の実験が有名である．その後，計算機の発明を期にそのアイデアはアルゴリズムとして Neumann らによってモンテカルロ法と名付けられ，20 世紀の十大アルゴリ
ズムの一つにあげられることになる．もっとも素朴に計算機に発生させた疑似乱数を用いて，厳密に解析計算することが困難な関数の数値積分法はその代表的な例である．
ここで議論するモンテカルロ法の重要な思想は，確率分布の大きな値を与える確率変数を効率よくサンプルする “importance サンプリング” にある．サンプリングができれば，そこから対象となる確率分布についての期待値が計算できる．物理学をはじめとして、ベイズ推定などの統計学の様々な問題はこの形式に帰着することが多く．
物性物理での多体問題 有限温度での相互作用系は平衡状態においてはカノニカル分布に従うが，一般に統計力学的な期待値の計算は解析的には困難である．平均場近似や各種摂動計算で物理的理解はできるかもしれないが，その妥当性の検証にはより正確に確率分布のことを知る必要がある．
データ解析におけるベイズ統計 実験によって得られたデータをある有効モデルにあてはめることを考えると，ベイズ推定ではデータによって条件付けられた有効モデルのパラメータの確率分布（事後分布）を調べることになる．いわゆる最小二乗法もベイズ統計の枠組みに入る．

Importance サンプリングの考えに基づくモンテカルロ法は大別して，マルコフ連鎖モンテカルロ
法 (Markov-chain Monte Carlo : MCMC) とポピュレーション型モンテカルロ法 (PMC) がある
1

E-mail: hukusima@phys.c.u-tokyo.ac.jp1

《講義ノート》

物性研究・電子版 Vol. 7, No. 2, 072214 （2018年11月号）

図 1: モンテカルロ法の分類の概略図 (左) と Importance サンプリングのイメージ図 (右)．通常は図で表すことが困難な超高次元確率分布を考えるが，ここでは単純な二次元確率分布を描く．
MCMC ではサンプル列が，PMC ではサンプル集団が近似的確率分布を表す．

(図 1(左))．あるサンプルしたいターゲットとなる確率分布が与えられたとすると，MCMC では
一連の状態 X のサンプル列を生成し，その経験分布がターゲット分布を近似する．一方，ポピュ
レーション型 MC 法では複数の重みを持ったランダムウォーカーが全体として確率分布を近似するように構成される．絵に描くのは容易でも，サンプル列の構成方法やランダムウォーカーの散布方法は自明ではない．この講義では MCMC 法に焦点を当てて，最近の話題まで解説する．図 2には MCMC 法の簡単な年表を記した．

図 2: マルコフ連鎖モンテカルロ法の開発に関する簡易年表

2

マルコフ連鎖モンテカルロ法

MCMC 法の目的は，高次元・多変量の確率分布からのサンプリングである．ここでは，確率変数のベクトルを X ，その確率分布を π(X) とし，ある与えられた π(X) に従う確率分布 X の生成方法を議論する．一般に確率変数は連続値をとりうるが，簡単のために離散値の場合を考える．
2

《講義ノート》

物性研究・電子版 Vol. 7, No. 2, 072214 （2018年11月号）

MCMC 法は確率分布で定式化される問題に適用することができ，サンプリングが必要なときや確率分布による期待値を評価するときの計算技法とみなせる．例えば，統計物理でのカノニカル分
布やベイズ推定での事後確率分布による期待値計算の際に，MCMC 法は有用になる．以下では，問題を特に限定せずに，一般論を展開する．

2.1

マルコフ連鎖の定常条件

MCMC 法の基本的な考え方は，確率変数 X のマルコフ連鎖を構成し，生成された一連のサンプル集合 {X} が対象となる確率分布 π(X) に従うようにする．ここで，簡単のためマルコフ連鎖
は，サンプル集合 X (1) , · · · , X (M ) の中で，X (m) は X (m−1) のみに依存する系列として定義する．
初期状態 X (0) と状態 X から X ′ への遷移確率 W (X → X ′ ) を指定するとマルコフ連鎖は規定される．理論的には確率分布の時間発展を考え，ある遷移確率で規定されるマルコフ連鎖において，遷移の前後で不変な確率分布をマルコフ連鎖の不変分布（定常分布）と呼ぶ．マルコフ連鎖の研究はたびたび与えられた遷移確率の元での不変分布を求めることであったりするが，MCMC 法を構成する問題はサンプリングしたい確率分布 π(X) に収束するように，遷移確率 W (X → X ′ ) の設計問題とみなせる．そのために遷移確率が満たすべき必要十分条件は以下の 2 つである．
(A) つりあいの条件 : 遷移確率は任意の X に対して関係式∑

π(X)W (X → X ′ ) =X′∑π(X ′ )W (X ′ → X),

(1)

X′

を満たす．この関係式は，ある状態 X に対して流出する確率の流れの合計 (左辺) が流入の合計と等しいことを意味する．

(B) エルゴード条件 : 任意の 2 つ状態 X と X ′ の間の遷移確率がゼロではないか，あるいは有限個のゼロでない遷移確率の積で表される．
以上の条件を満たせば，マルコフ連鎖は初期条件に依らず，唯一の不変分布に収束することが示される．
もう少しだけ補足すると，マルコフ連鎖の n ステップのときの確率分布を Pn (X) とすると，離散時間発展のもとでは 1 ステップ後の確率分布 Pn+1 (X) は

Pn+1 (X) =∑Pn (X ′ )W (X ′ → X)

(2)

X′

となる．これはマスター方程式と呼ばれる．不変分布 π(X) は P∞ (X) と一致するためには、

π(X) =∑

π(X ′ )W (X ′ → X)

(3)

X′

∑

を満たす必要がある．一方，確率の保存則 X ′ W (X → X ′ ) = 1 を式 (3) の両辺にかけ合わせると，つりあい条件∑∑
π(X)W (X → X ′ ) =π(X ′ )W (X ′ → X),
(4)
X′

X′

が導かれる．確率行列 W に対する Perron-Frobenius の定理より，最大固有値は 1 であり，縮退はなく，その固有ベクトルは不変分布になることが示される．それ以外の固有値は真に１より小さく，特に第二固有値の大きさが不変分布への収束の速さを決めている．この第二固有値の小さい
遷移確率がこの意味では効率のよいといえるが，どのように設計すれば第二固有値が小さくなるかは自明ではない．

3

《講義ノート》

物性研究・電子版 Vol. 7, No. 2, 072214 （2018年11月号）

図 3: 確率の流れの概念図：(左) 詳細つりあい条件．(右) つりあい条件

2.2

詳細つりあい条件と Metropolis filter

MCMC 法を構成するために，つりあい条件を満たす遷移確率を求めたいが，条件数が少なくて一般的な構成法は知られていない．これまでは，より条件の強い詳細つりあい条件
π(X)W (X → X ′ ) = π(X ′ )W (X ′ → X)

(5)を課すことが慣例的である．非ゼロの確率が存在する遷移 (X → X ′ ) に対して，順方向の確率の流れと逆向きの確率の流れがつりあうことを要請する (図 3)．当然、これを満たされるときにつりあい条件 (1) は満たされる．このとき，生成されるマルコフ連鎖は可逆 (reversible) であるという．
求めたい分布 に対して，詳細つりあい条件を満たすように遷移確率を選択することは多くの場合に容易である．代表的な遷移確率にメトロポリス型の遷移確率 (Metropolis filter) がある：(

WM (X → X ′ ) = min 1,)

π(X ′ )
.
π(X)

(6)

不変分布の確率が大きな状態への遷移は確率 1 で起きて，そうでないときにも有限の確率で実現することを意味している．Importance サンプリングとは言っても, 必ず確率を大きくする方向にマルコフ連鎖が進むわけではない. また，この右辺は不変分布の比のみの関数であるので、確率分布の規格化定数は現れていない．ベイズ統計での事後分布や統計力学でのギブス分布などでは規格化定数が不明なことが多く，その値を知る必要がないことは一つの大きな利点である．一方，規格化定数には周辺尤度や熱力学関数としての重要な意味があるために，評価したい量ではある．
このことについては後ほど議論することにしよう．
詳細つりあいの条件は不変分布への収束性を保証するための必要条件ではなく，必ずしも詳細
つりあい条件は満たさなくてよい．この事実は昔から知られていたが，現実的な構成方法が明示することは Suwa-Todo 法 [14, 15] まで待つことになる．最近，この詳細つりあい条件を破る MCMC法はいくつか提案されているが，その実質的な効率も含めて現在研究が進められている．

2.3

最適化手法としての Simulated Annealing

MCMC 法の構成方法の単純さが広範な応用を可能にする原因の一つと考えられるが，様々な問題に適用する中である共通の困難点が明らかになってきた．それは，マルコフ連鎖において遷移
確率 W (X → X ′ ) が小さくなり，実質的に遷移が行なわれなくなることである．結果として，サンプルが確率変数の空間で偏りが大きくなり，また初期条件に強く依存し，期待値の評価に大き4

《講義ノート》

物性研究・電子版 Vol. 7, No. 2, 072214 （2018年11月号）

な系統誤差を与えることもある．この現象は遅い緩和の問題と呼ばれている．この問題を解決するためのヒントは最適化手法の開発に見つけることができる．
MCMC 法を最適化問題の解法として発展させた手法にアニーリング法 (Simulated annealing)がある [2]．この方法は最適化問題をはじめ，ベイズ統計での事後分布のモードや物理学の問題での最低エネルギー状態の汎用的探索手法と考えられ，一般に確率分布 π(X) の最大値を与える状
態 X を求めるために用いられる．確率分布が多峰的なときには − log π(X) に対する最急降下法
を用いると，初期条件に依存した局所最適解に収束してしまい，大域的な最適解にたどり着かな
いことがたびたび起こる．これは MCMC 法にも共通する困難点である．アニーリング法では，温度パラメータ T を導入して，変形した分布

π 1/T (X)
1/T (X)Xπ

πT (X) = ∑(7)を考え，MCMC 法を用いて X をサンプルする．十分大きな T では πT (X) はほぼ一様分布とみなせ，マルコフ連鎖に現れる状態はランダム抽出となる．一方で，小さな T では π(X) のモード付
近からサンプルされることが期待される．アニーリング法は大きな T での確率的なゆらぎによって，多峰的な分布の局所解に留まり続ける可能性を低くし，徐々にパラメータ T を小さくすることで，最適解の探索を可能にする手法である．パラメータ T を変化させるときに詳細つりあい条
件を破るために不変分布からのサンプル手法ではないが，多峰的な分布でのサンプリングの困難点を回避する考え方は拡張アンサンブル法へ繋がっていく．また，このような温度由来のゆらぎ
ではなく，量子力学起源の量子ゆらぎを利用した量子アニーリング法の研究は近年 D-wave による実装機械の出現に伴い著しく進んできている．

3

拡張アンサンブル法

多峰的な分布に現れる MCMC 法の困難点を解消するために，マルチカノニカル法 (Multicanonical method)[3]，Simulated Tempering(ST)[4] や交換モンテカルロ法 (Exchange MC, multiplecoupled MCMC, paralle tempering)[9] 等の一連の方法群が 1990 年代に主に統計物理の分野で相次いで提唱された [6, 7]．これらの方法は，興味ある確率分布に修正を加えたり，合体した拡張された分布を取り扱うことから拡張アンサンブル法と呼ばれている [6]．ここでは ST 法と交換法の構成方法を見ておくことにする．

3.1

Simulated tempering と交換法

アニーリング法をサンプリング手法に自然に拡張した方法が Simulated tempering(ST) 法である [4, 5]．名前もなかなか洒落ている．ここでは説明のためにエネルギーが E(X) と与えれれているギブス分布を考える:
1
exp (−βE(X)) .
(8)
π(X) =Z(T )ST の特徴は，温度 T も確率変数とする．具体的に T1 , T2 , · · · , TR の値をとりうるとする．拡張された状態は (X, T ) と表し，その不変分布を
πST (X, Tr ) =1
exp (−βr E(X) + gr )
ZST

5

(9)

《講義ノート》

物性研究・電子版 Vol. 7, No. 2, 072214 （2018年11月号）とする．ここで gr は Tr に依存する重み因子であり，ZST はこの確率分布の規格化定数である．ある与えられた Tr の元での X の確率分布はギブス分布と同じである．一方，温度 Tr の周辺分布は、

πST (Tr ) =1 ∑Z(Tr )egr
exp (−βr E(X) + gr ) =
ZST X
ZST

(10)

となる．もし gr = − log Z(Tr ) とおくと，この周辺分布は温度に依らなくなる．一般に、これは系の自由エネルギーに相当し，評価の難しい量である．しかし、それが可能であるとすると，このマルコフ連鎖では温度に関して一様分布になる．すなわち，高温も低温も同じように出現することが期待される．低温度で多峰的な複雑な分布に引っかかったとしても，高温へ遷移することで脱出することができそうである．具体的に詳細つりあい条件を満たすアルゴリズムを示しておこう．

Step 0 適当な初期条件 (X (0) , T (0) ) を決める
Step 1 温度 T を止めて状態 X を更新する．ここは通常の MCMC 法を行う．
Step 2 状態 X を止めて，温度 T を更新する．
1. 現在の状態を (X, Tr ) として，次の温度 Tn を前後の温度 Tr−1 と Tr+1 を確率 1/2 で提案する．温度の境界は適当に定義することにする．
2. その候補 Tn を採択確率(

πST (X, Tn )A(Tr , Tn ) = min 1,
πST (X, Tr )

)

(11)

で採択する．
アニーリング法では温度は一方的に下がるだけであったが，ST 法では下がることもあれば上がることもある．アニーリング法のよいところを引き継ぎながら、全ての過程で詳細つりあい条件を満たしているためにサンプリング手法としての意味を持つ．これでよいように思えるが、gr は事
前には知らないことが問題として残っている．この推定法として徐々に学習する方法が沢山提案されている．
拡張アンサンブル法の中で最も単純なものは交換モンテカルロ法であるが，これは先の ST 法の問題点を回避する方法として考案された．交換法では，R 個の異なる温度 (T1 , T2 , · · · , TR ) を考え，確率変数 {X} = {X1 , X2 , · · · , XR } と R ケのレプリカを扱うことになる．ST 法と異なり温度は確率変数ではない．このとき拡張された同時確率分布 πEX ({X}) は，πEX ({X}) =R∏

πTr (Xr )

(12)

r=1

とする．πTr (X) は温度 Tr で特徴づけられる確率分布である．交換法のマルコフ連鎖は，この同時確率分布を不変にするように二種類の状態更新を導入する．一つは，それぞれの温度の確率分
布 πTr (X) に対して状態 X の更新を行う．もう一つは，異なる温度の状態の交換である．すなわ
ち，2 つの温度 Tr と Tr+1 の状態 Xr と Xr+1 を交換するのである．式 (12) を不変にするためには，交換採択レート REX を

REX =
=

πEX (X1 , X2 , · · · , Xr+1 , Xr , · · · , XR )πEX (X1 , X2 , · · · , Xr , Xr+1 , · · · , XR )πr (Xr+1 )πr+1 (Xr )πr (Xr )πr+1 (Xr+1 )6(13)

《講義ノート》

物性研究・電子版 Vol. 7, No. 2, 072214 （2018年11月号）とし，メトロポリス型ならば交換採択確率を min(1, REX ) とすればよい．ギブス分布の交換レー
トは，
)
(
1
1
−
(E(Xq ) − E(Xq+1 )),(14)
REX = exp
Tq
Tq+1となり，交換すべき 2 つのレプリカの温度差とエネルギー差から計算できる．これらの 2 つの過程により，各レプリカの変数はそれぞれの温度で更新をしながら，同時に温度軸上を動き回り，長時間の振舞はランダムウォークとみなすことができる．つまり，自分自身でアニーリングとヒー
ティングを行っているように見える (図 4)．結果として，低温で多峰的な確率分布のあるピークに留まり続けることがあったとしても，一度高温になり一様分布を介することにより広く探索できることを期待するわけである．ただし，いわゆるアニーリング法と異なり，全ての素過程で詳細
つりあい条件を満たしていることから，それぞれの温度では πT (X) からのサンプリングが実現でき，ある温度 Tr での期待値は通常のように πTr (Xr ) の平均として計算できる．ST 法のように重み因子 gr を導入することもできるが，この交換レートでは比をとることで表に出てくることがな
い．この意味で ST 法の困難点は回避できているが，その代償として複数のレプリカを同時に平行計算する必要がある．その計算量の多さが弱点となっているが，現在の計算機の主流が超並列計算に移行し，大規模計算の分野では並列計算することになる．
(0)

X1

(0)

X2

(0)

X3

(0)

X4

(0)

(1)

X1

(1)

X2

(1)

X3

(1)

X4

(1)

(2)

X1

(2)

X2

(2)

X3

(2)

X4

(2)

(3)

X1

(3)

X2

(3)

X3

(3)

X4

(3)

X5

X5

X5

X5

(0)

(1)

(2)

(3)

X6

X

X

X

T1
T2
T3
T4
T5
T6

MC Step

図 4: R = 6 レプリカ系での交換法の状態更新の様子．

3.2

近似的数え上げと自由エネルギー計算

ここまで，拡張アンサンブル法の効用として，遅い緩和問題の解消方法として述べてきたが，その他に確率分布の規格化定数あるいは自由エネルギーの計算法としての重要な利点がある．通常の MCMC 法では，規格化定数の計算のためには，確率分布をあるパラメータについての対数微分で定義される観測量を計算し，そのパラメータについての積分を施す手順を採る．ギブス分布については自由エネルギーは熱力学積分と呼ばれる方法である．そこでは，自由エネルギーは，∫ βln Z(β) =
0

dβ ′

∫

β
d
′
dβ ′ ⟨E(β ′ )⟩
ln
Z(β
)
=
−
dβ ′
0

(15)

と表され，右辺のエネルギー期待値が MCMC 法での観測できる量である．交換法ではこの熱力学積分のための複数の温度の情報が得られ，マルチヒストグラム法を併用することにより効率よく計算できる．また，マルチカノニカル法では，重み関数に近似的に求めた状態密度から，再重7

《講義ノート》

物性研究・電子版 Vol. 7, No. 2, 072214 （2018年11月号）

0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

1
72
2340
45456
589158538675235826516
176198256
645204321
1758028568
3538275120
5185123200
5409088488
3885146784
1829582496
524514432
81145872
5415552
90176

1e+10

Exact
MC

1e+08
1e+06
10000
100
1
0

2

4

6

8
10 12
# of dimers"

14

16

18

図 5: 二次元格子上のダイマー敷き詰め問題：(左) 典型例．(中) 厳密数え上げによる敷き詰めパ
ターン数のダイマー数依存性 [13]．(右) マルチヒストグラム法を用いた敷き詰めパターン数のダ
イマー数依存性．
率法を用いることで規格化定数相当の量が計算できる．これらは，問題の規模の指数関数以上ある大きな数を数え上げる一般的方法とみなせる．また，この方法は物理では自由エネルギー，ベイズ統計では対数尤度が精度よく評価できる方法であることを意味していて，拡張アンサンブル法を用いた自由エネルギー地形を描く標準的な手法と言える．
ここでは二次元格子上のダイマー敷き詰め問題の例を紹介する. 正方格子上に長さ１の棒を置くことを考える．棒の数を指定して，その置き方の場合の数を数える問題は状態数の数え上げとみなせる．6 × 6 の格子に関しては厳密に数え上げることができる [13]．その結果は図 5(中) に示す．
その右の結果はマルチカノニカル法の結果である．強調したいことは，計算しているモンテカルロ
法のステップ数は高々105 程度であるにもかかわらず，1010 近くの場合の数が精度よく評価できていることである．もちろん，数え上げているわけではない．このような拡張アンサンブル法による数え上げの研究は，魔法陣の数 [16, 17], N クイーン問題 [8], ラテン方陣 [6] などがある．また，ベイズ推定には変数選択と呼ばれる重要な問題があるが，この近似数え上げ法を応用して，変数選択の全ての組合せに対する評価誤差についての状態密度の計算が可能になってきている [25]．

4

最近の話題–詳細つりあいの破れた方法，Lifting，Event-chain MC–

近年，MCMC 法の発展として，少し興味を持たれているのが詳細つりあい条件を破る MCMC法の発展である．前章で議論したように，つりあい条件を満たせばよいが，遷移確率の構成法の容易さから強い条件である詳細つりあい条件を用いることが多かった．一方で，詳細つりあい条
件は状態空間での遷移経路毎に確率の流れがつりあっているために，必然的に拡散的な確率の流れしか期待できないが，もしそれを破ることで確率の流れに方向性のある流れが生じさせることができれば不変分布に速く収束するかもしれないことは想像できる．コーヒーに砂糖を入れたと
き、じーと拡散的に解けることを待つよりも、スプーンで混ぜた方が速く解けそうだというわけである．

8

《講義ノート》

物性研究・電子版 Vol. 7, No. 2, 072214 （2018年11月号）

この流れの最初の重要な仕事が Suwa-Todo により与えられた [14]．そこでは基本的には Peskunの定理的な思考のもとで，棄却率を減少させる遷移確率の構成方法を与えた．その後，Lifting と呼ばれる概念に基づく幾つかの方法が提案されてきた [18, 19, 21, 22, 23]．現時点では，詳細つりあい条件を破る MCMC 法は大きくこの二つに分かれるように思われる．Lifting とは，仮想的に状態空間を拡大して，そこで確率の流れにバイアスをかける考え方を指す．
この Lifting の考えの延長線にある最近の発展に Event-chain Monte Carlo(ECMC) 法がある．
この ECMC 法は最初は剛体球のモンテカルロ法として導入された [21]．周期的境界条件の空間で剛体球は必ず右方向にしか動かないとする特殊な遷移を考える．その定義から左に進む遷移はな
いので，詳細つりあい条件はやぶられている．しかし，うまく玉突き運動を構成することでつりあい条件を満たそうする方法である．運動する球を表す変数を状態変数に加えている点が Liftingの考えに即していると解釈されている．剛体球の密度が大きいときに，拡散的な運動では大きく配置を変えられないが，玉突き運動を続けることで状態が大きく変化しうるわけである．
その後，ポテンシャルで表される相互作用系に ECMC 法は拡張された [22]。この拡張はまった
非自明な技法の導入が必要であった．剛体球のときの衝突に相当する Event を相互作用系で同定することの難しさは少し考えてもわかる．しかし、この拡張によって、かなり広いクラスの確率分布に応用できるようになってきたと思われる．まだ技術的な問題のために，上記の条件を満たす全ての問題に適応できるわけではない．我々もこの方法を発展させて，統計物理の連続スピン系へ応用した [24]．

5

最後にここまでマルコフ連鎖モンテカルロ法を簡単に解説した．私が最初にモンテカルロ法に触れたのは学部 4 年の卒業研究のときであった．最近のみなさんと比較すると，やや奥手かもしれない．
そもそも個人が自由につかえるコンピュータがそれほど身近にあったわけではない時代であった．
最初に読んだ本格的な論文がスピングラスの有名なモンテカルロ計算の論文であった．当時，話題になった素晴らしい研究であったらしいが，それを見た私は絶望的になった記憶が今でも残っている．論文の中身が理解できるようになると，このスピングラスという系ではモンテカルロ法での相関時間が驚くほど長いらしいことがわかった．これでは普通の計算機では研究ができないと思ったわけである．面白い問題が目の前にあるのに，手も足もでない状況であった．そこで，ス
ピングラス系でも研究可能なモンテカルロ法の方法論の開発に取り組んだ．それが私の博士論文のテーマになった．
近年，機械学習が流行っていることもあり，データ解析に注目が集まっている．その中で，モンテカルロ法は一つの有力な解析方法である．この講義でもデータ駆動科学の事例を紹介するが，あくまでも解析手法は方法であって，目的が一番重要である．道具を磨くことには美味しい食材が必要である．それを忘れて，道具だけ磨いているようでは本当の研究の醍醐味や楽しさが消えてしまう気がする．この講義ではみなさんとそこそこ切れ味のよい道具を共有することを目標として，みなさんの興味のある問題がうまく料理できる助けができれば幸いである．

参考文献[1] N. Metropolis, A.W.Rosenbluth, M.N.Rosenbluth, A.H.Teller and E.Teller, J. Chem. Phys.
21, 1087/1091 (1953)

9

《講義ノート》

物性研究・電子版 Vol. 7, No. 2, 072214 （2018年11月号）[2] S. Kirkpatrick, C.D.Gelatt, and M.P.Vecchi, “Optimization by simulated annealing”, Science 220 67/680, 1983.
[3] B. A. Berg and T. Neuhaus: Phys. Lett. B267 , 249 (1991).
[4] E. Marinari and G. Parisi, Europhys. Lett. 19, 451 (1992).
[5] A. P. Lyubartsev, A. A. Martinovski, S. V. Shevkunov, and P. N. Vorontsov-Velyaminow,
J. Chem. Phys. 96, 1776 (1992).
[6] Y. Iba: Extended Ensemble Monte Carlo, Int. J. Mod. Phys. C, 12, 623/652, (2001).
[7] J. S. Liu, Monte Carlo Strategies in Scientific Computing, Springer-Verlag New York (2001).
[8] K. Hukushima, Comp. Phys. Comm., 147, 77/82 (2002).
[9] K. Hukushima and K. Nemoto, J. Phys. Soc. Jpn., 65, 1604/1611, (1996).
[10] 伊庭幸人他, 計算統計 II マルコフ連鎖モンテカルロ法とその周辺，岩波書店，2005．
[11] K. Hukushima, Phys. Rev. E 60, 3606, 1999.
[12] F. Wang, and D.P. Landau, Phys. Rev. Lett. 86, 2050, 2001.
[13] W.Krauth, Statistical Mechanics: Algorithms and Computations, Oxford Master Series inPhysics (2006).
[14] H. Suwa and S. Todo, Phys. Rev. Lett. 105, 120603 (2010).
[15] S. Todo and H. Suwa, J. Phys.: Conf. Ser. 473, 012013 (2013).
[16] K. Pinn and C.Wieczerkowski, Int. J. Mod. Phys. C9, 541 (1998).
[17] A. Kitajima and M. Kikuchi, PLos ONE 10 (2015), DOI: 10.1371/journal.pone.0125062.
[18] K. S. Turitsyn, M. Chertkov, and M. Vucelja, Physica D 240, 410 (2011).
[19] Y. Sakai and K. Hukushima, J. Phys. Soc. Jpn. 82 064003 (2013).
[20] Y. Sakai and K. Hukushima, arXiv preprint arXiv:1511.08100 (2015).
[21] E. P. Bernard, W. Krauth, and D. B. Wilson, Phys. Rev. E 80, 056704 (2009).
[22] E. P. Bernard and W. Krauth, Phys. Rev. Lett. 107, 155704 (2011).
[23] M. Michel, S. C. Kapfer, and W. Krauth, J. Chem. Phys. 140, 054116 (2014).
[24] Y. Nishikawa, M. Michel, W. Krauth, and K.Hukushima, Phys. Rev. E 92, 063306 (2015).
[25] Y. Igarashi, H. Takenaka, Y. Nakanishi-Ohno, M. Uemura, S. Ikeda, and M. Okada, J.
Phys. Soc. Jpn. 87, 044802 (2018)10