A Blog Entry on Bayesian Computation by an Applied Mathematician
$$
$$
1 概要
相互作用する分子系からなる物質の状態方程式などの性質を調べるために使える汎用手法（ fast computing machine にぴったり！）を提案する．
この手法は配置空間上の，提案分布を正確にした「修正版 Monte Carlo 積分法」だと捉えられる，としている．
現代的には，１度に一つの粒子ずつを動かすことで，棄却率の低下を防いだ Metropolis-within-Gibbs サンプラーと理解できる．1
MANIAC を用いて，２次元剛体円板の液相転移のシミュレーション結果を付した．その結果を，自由体積状態方程式と，４項ビリアル係数展開と比較した．2
一方，(Wood and Parker, 1957) ではこの手法を３次元の Lennard-Jones ポテンシャルに適用しており，相転移の痕跡を掴んでいる．
この手法を Ising 模型に応用した (Glauber, 1963) は (random scan) Gibbs サンプラーの先駆けとみなされている．
2 背景
2.1 サンプリング法としては見られていない
むしろ cell method などの計算手法との類似で捉えており，(Hastings, 1970) に取り上げられるまで全く別個の，広く抽象的に応用されるポテンシャルを持つサンプリング法であるという抽象的な観点とは離れた泥臭い探索を感じる．
まだ Markov 連鎖という単語も使われていない．
2.2 分子シミュレーションの延長としての Monte Carlo 法
話の展開は，分子系のシミュレーションの問題を，Gibbs 分布に関する積分計算の問題に帰着し，Monte Carlo 法を用いるところまで還元するが，そこで不適当な提案分布による重点サンプリングを実行するよりかは，系の Gibbs 分布を直接シミュレーションすれば良い，というものである．
よってここでいう修正 Monte Carlo 法とは，Monte Carlo 法（ここでは重点サンプリングと同義，（一様）乱数を用いた数値計算，くらいの意味）による乱数生成過程を，少しだけ MD 法の考え方を取り入れて，より Gibbs 分布に則した乱数生成をする，くらいの提案である．
後世では次のようにも説明されている始末である：
The Monte Carlo method addresses the sampling problem more abstractly than molecular dynamics, as it samples (obtains samples \(x\) from) the distribution \(\pi_{24}(x)\) without simulating a physical process. (Tartero and Krauth, 2023)
もはや，愚直な Monte Carlo 法が MD 法で，運動方程式を解かない MD 法が Monte Carlo 法，という具合である．事実，気体の状態方程式は，気体分子の運動の力学には無関係に成り立つので，必ずしも正しい力学に従ったシミュレーションが必要というわけではないのである．3
なお，Monte Carlo 法の提案者として Joseph Edward Mayer と Stanisław Ulam の名前を挙げている．
ここで Berni Alder の名前が上がっており，やはり MCMC の開発は MD と極めて深い関係にあることが伺える．
This method has been proposed independently by J. E. Mayer and by S. Ulam. Mayer suggested the method as a tool to deal with the problem of the liquid state, while Ulam proposed it as a procedure of general usefulness. B. Alder, J. Kirkwood, S. Frankel, and V. Lewinson discussed an application very similar to ours.
2.3 Monte Carlo 法の起源について
John von Neumann は Edward Teller と同郷であり，戦時中の ENIAC の開発にも関わっていたことから，偉い人たちを説得して最初の ENIAC のテストに熱核融合反応の計算問題を用いることにした．
計算可能なモデルの構築をしたのが Metropolis である．これが完成し，実際にテストが行われたのは戦後の 1946 年であったが．
そのお披露目会に居合わせた Stanislaw Ulam が，統計的サンプリング技術を電子計算機で復活させることを提案し，Johnny がすぐさまその重要性を理解した．これがモンテカルロ法の始まりとなった，という (Metropolis, 1987)．
Metropolis は Ulam が着想を得た理由として，数学的な背景を持っていたために，統計的サンプリング技術が計算の難しさのために歴史に埋没したことを知っており，ENIAC のポテンシャルを見てこれと関連づけることに成功したのではないかと示唆している．
そして Johnny の熱の入り用が周りも刺激した．1947 年には統計的サンプリングの中でも特に中性子の拡散問題を取り上げて当時の Los Alamos の理論部リーダーであった Robert Richtmyer に手紙を送った (Eckhardt, 1987)．こうして周りを巻き込んで大ごとになっていった．Monte Carlo 法という命名も 1947 年だったという
It was at that time that I suggested an obvious name for the statistical method—a suggestion not unrelated to the fact that Stan had an uncle who would borrow money from relatives because he “just had to go to Monte Carlo.”
On a less grand scale these events brought about a renascence of a mathematical technique known to the old guard as statistical sampling; in its new surroundings and owing to its nature, there was no denying its new name of the Monte Carlo method. (Metropolis, 1987)
特に戦時中の関心もあり，核分裂時の 中性子の拡散 のシミュレーションが問題であった．Monte Carlo 法と呼んでいるが本質的に MD 法チックであり，ここの中性子の散乱・吸収・分裂の系譜をシミュレートすることで全体の統計的性質がわかる，というだけの話であった．
その後 1952 年には後続機の MANIAC が開発され，nucleaer cascade と状態方程式も射程に入った．
この状態方程式を取り扱う際に (Metropolis et al., 1953) がさらに効率的な「モンテカルロ法」を発明したのである．これは Gibbs 分布を直接シミュレーションできるというブレイクスルーであり，「一般の確率分布からサンプリングできる」という今の理解とは大きく異なる文脈の中で発見されたと言うべきである．
2.4 Monte Carlo 法とはなんだろうか
思うに，「ランダムな方法を使って計算する」というのは外道に思えるかもしれない．
だが，実はランダムな系の \(\mathcal{P}(E)\) 上のダイナミクスの決定論的な計算になっているのかもしれない．
そう思わせるだけの透徹性が測度論にはある．
ただ，(Metropolis, 1987) は Monte Carlo 法を 実験数学 (experimental mathematics) と呼んでおり，極めて物理学的な見方で評している：
At long last, mathematics achieved a certain parity–the twofold aspect of experiment and theory–that all other sciences enjoy.
2.5 他のコメント
Note that (Metropolis et al., 1953) move one particle at a time, rather than moving all of them together, which makes the initial algorithm appear a primitive kind of Gibbs sampler! (Robert and Casella, 2011)
この点は，Metropolis 法と MD 法の違いとも通じる．
Metropolis の方法は１つの粒子を順番に動かすが，MD 法では全ての粒子を同時に動かす．
加えて，Metropolis 法は平衡状態のシミュレーションのみを想定しているという点で，本質的にはサンプリング法であるが，MD 法は非平衡系の時間発展も同等に扱うシミュレーション法である．4
本質的に Metropolis 法がサンプリング法であるならば，MD 法は \(N\)-体問題に対する数値解法である．
3 本論
3.1 設定
古典統計を仮定し，２体間相互作用のみを考え，ポテンシャルは球対称であるとする（流体力学では通常の仮定である）．だが，温度や密度には全く仮定を置かない．5
実際の計算のために，粒子数 \(N\) は several hundred に取る．そして正方形の中にいれ，境界条件を最小化するために同様の系が２次元に無限に連なっているとする．２つの粒子 \(A\) の他の粒子 \(B\) との最短距離を \(d_{A,B}\) とし，これのみが粒子 \(A\) にかかる主な力になるとする．
仮に \(N=1\) だとしたら，これは cell method と呼ばれるモデルでもある．こうして粒子を増やすことで，単一相のシステムに対するより良いモデルになるだろうが，二相以上のシステムには限界がある．
以上の仮定から，系のエネルギーが次のように与えられる： \[ E=\frac{1}{2}\sum_{i\ne j\in[N]}V(d_{ij}). \]
この系の平衡状態の性質を計算するには，Gibbs の正準分布を利用し，計算したい物理量 \(F\) に対して \[ \overline{F}=\frac{\int Fe^{-\frac{E}{kT}}d^{2N}pd^{2N}q}{\int e^{-\frac{E}{kT}}d^{2N}pd^{2N}q} \] を計算すれば良い．ただし，\(d^{2n}pd^{2n}q\) は \(4N\) 次元相空間上の体積要素である．
加えて，ここではポテンシャル \(V\) は位置のみの引数としているから，\(2N\) 次元上でのみ計算すれば良い．
このような数百次元上での積分を数値的方法で実行するのは明らかに実行可能でないから，Monte Carlo 法に頼らざるを得ない．と言っても，決定論的な点で値を計算する代わりに，ランダムに点をうつ，というだけの違いではある．
最も簡単な実装としては，ランダムに \(N\) 粒子を配置してエネルギーを計算し（重点荷重），これにウェイトをつけて足していくということが考えられる（重点サンプリングだ！）．しかし，高エネルギーの配置もたくさん生成してしまうから，これによる効率の低減が避けられない（提案分布が悪いのだ！）．
そこで我々は modified Monte Carlo method を考える．そもそも確率 \(e^{-\frac{E}{kT}}\) からサンプルを生成し，荷重を一様にすることを目指す．
3.2 アルゴリズムの記述
これは次のようにする．まず適当に初期分布を決める（格子点上に \(N\) 粒子を配置するなど）．そしてこれをアップデートしていく： \[ X\mapsto X+\alpha\xi_1 \] \[ Y\mapsto Y+\alpha\xi_2 \] \(\alpha>0\) は一度にどれくらい動かすかを調節するパラメータであり，\(\xi_1,\xi_2\in(0,1)\) は一様乱数とする．6
すなわち，\((X,Y)\) を中心とした一辺 \(2\alpha\) の正方形の中で，新たな位置をランダムに決めるのである．7
この動きによるエネルギーの変化量 \(\Delta E\) を計算し，\(\Delta E<0\) ならばこれを実行するが，\(\Delta E>0\) ならば確率 \(e^{-\frac{\Delta E}{kT}}\) によって採択する．
仮に棄却されたとしても，そのポジションから新たな Monte Carlo 標本を取り，最終的に \[ \overline{F}=\frac{1}{M}\sum_{j=1}^MF_j \] を推定量とする．
3.3 アルゴリズムの有効性の検証
まずこの系はエルゴードであると主張しているが，その論証は「任意の粒子が任意の位置に行くポテンシャルがあるため，この手法はエルゴード的である」で終わっている．エルゴードという単語を「任意の状態からもう一つの任意の状態に遷移可能である」という意味で使っている．これは現代的には既約性という．8
続いて，この系をたくさんコピーしてアンサンブルを考えたとき，状態 \(r\) にいるアンサンブルの数 \(\nu_r\) は \[ \nu_r\,\propto\,e^{-\frac{E_r}{kT}} \] を満たすことを示したい．その論証は，上の比率から崩れていたら，平衡に至る方向へ移動が起こるということを具体的に議論している．
以上の２点から，提案されたアルゴリズムは正準分布に収束することの根拠としている．
このアンサンブルによる考え方は極めて直感的に訴える．実際，この語彙を用いて，棄却された場合は元々の状態を Monte Carlo サンプルとしてダブルカウントすべきであることを説明している．これをしなければ，低エネルギー状態のアンサンブルの数を不当に低く評価してしまう，という説明である．
ただし，このアンサンブルによる考え方は自然に我々の思考を詳細釣り合い条件に絞っている．
3.4 附言
収束の速さについて注意喚起しているのみで，ステップサイズ \(\alpha>0\) は大きすぎても棄却率が高まり，小さすぎても攪拌が遅くなるということ以外具体的なことは触れていない．
4 実験
In the case of two-dimensional rigid spheres, runs made with 56 particles and with 224 particles agreed within statistical error. For a computing time of a few hours with presently available electronic computers, it seems possible to obtain the pressure for a given volume and temperature to an accuracy of a few percent. In the case of two-dimensional rigid spheres our results are in agreement with the free volume approximation for \(A/A_0< 1.8\) and with a five-term virial expansion for \(A/A_0> 2.5\). There is no indication of a phase transition.
16 step を焼き入れとし，48-64 いてレーションを実行するのに，MANIAC で 4-5時間かかったという．
References
Footnotes
(Faulkner and Livingstone, 2024, p. 13) 4.1 節も参照．↩︎
このシミュレーションでは相転移は全く観察されなかった．これはこの方法が極めて収束が遅いことによる．その理由は (Peters and de With, 2012) が導入部で説明する通りである．↩︎
(Faulkner and Livingstone, 2024, p. 15) 4.3 節も参照．↩︎
さらに，Lennard-Jones ポテンシャルについての２次元のケースも考えており，次のレポートで報告される予定であるという．↩︎
ここで乱数生成法に関して注記されている．これは middle square process で生成する，としている．↩︎
なお，periodic assumption をしているため（長方形の系の外には同様の系が無数に並んでいるとしたため），境界の外に出ようとした場合は衝突するのではなく，反対側の辺から入ってくるものとする．↩︎
(Robert and Casella, 2011) も指摘している．だが，(戸田盛和 et al., 2011, p. 5) にも同じ意味で「エルゴード的」という単語を使っている記述がある．↩︎