平衡状態を探す：マルコフ連鎖 – CFTP

来嶋秀治・松井知己

平成 16 年 8 月 25 日

マルコフ連鎖は確率的な挙動をするシステムを表す基本的なモデルである．本稿ではマルコフ連鎖の平衡状態を完璧に実現してしまう「パーフェクトシミュレーション」について紹介しよう．これは，マルコフ連鎖の無限回の推移を実現してしまう，まさに驚天動地のシミュレーション手法である．

マルコフ連鎖

以下，有限のマルコフ連鎖 M について考える．マルコフ連鎖 M は N 個の状態で構成される状態空間 Ω = {s_{1}, s_{2}, …, s_{N}} と推移確率行列 P を持っているとしよう．すなわち， M が各時刻で取り得る状態は Ω 上のいずれかであり，時刻 t で状態 x ∈ Ω にいるとき，時刻 t+1 の状態が y ∈ Ω である確率は Pr(x から y へ推移) = P_{xy} である．

[figure: s_{1}; s_{2}; s_{3}; 0.4; 0.4; 0.4; 0.4; 0.4; 0.2; 0.2; 0.2; 0.4]

図 1: マルコフ連鎖の絵．状態空間は Ω = {s_{1}, s_{2}, s_{3}} で，矢印に書かれている数字は推移確率を表している．また，定常分布は π = (1/3, 1/3, 1/3) である．

マルコフ連鎖 M の状態空間 Ω 上の確率分布 π = (π(s_{1}), π(s_{2}), …, π(s_{N})) について考える．このとき π は ∑_{x∈Ω} π(x) = 1 および π(x) ≥ 0 (∀x ∈ Ω) を満たす．この π が π = πP を満たすとき，定常分布という．今後，マルコフ連鎖 M は唯一の定常分布 π を持つと仮定する．

マルコフ連鎖を無限回推移させると，その状態は定常分布 π に従うことが知られている．冒頭で述べたマルコフ連鎖の平衡状態とはマルコフ連鎖が定常分布に従う状態を指している．

したがって，マルコフ連鎖を無限に等しい回数推移させると平衡状態に達するのであるが，ここにひとつ問題がある．「無限に等しいとみなせる回数とは何回なのであろうか？」

マルコフ連鎖の定常分布を直接計算することで，その平衡状態を実現することも可能である．しかし，この方法は往々にして実用的ではない．というのも，実用上のマルコフ連鎖では，状態数 N が非常に大きく，推移確率行列 P や定常分布 π を保持することすら難しいのである．そのような例はいろいろあるが [3]，後でイジングモデルと呼ばれる分子運動モデルを紹介する．

1996 年に Propp and Wilson はマルコフ連鎖のシミュレーションを工夫することで，パーフェクトシミュレーションを実現できることを示した [4]．そのアイデアが CFTP (Coupling From The Past) である．このアルゴリズムは「無限回は何回？」という問いに対する答えが分からなくても実行できる，それどころか，有限回で打ち切るアルゴリズムに特有の歪みが一切生じない，まさに夢のようなアルゴリズムである [1, 2]．

CFTP の説明をする前に，準備として更新関数のアイデアを紹介しよう．

更新関数

マルコフ連鎖 M の推移規則は，現在の状態 x と区間 [0, 1) から一様に発生させた乱数 Λ を引数とする関数 ϕ : Ω × [0, 1) → Ω によって記述できる．この関数は更新関数と呼ばれる．更新関数の満たす性質として，乱数 Λ ∈ [0, 1) に対して，

Pr(ϕ(x, Λ) = y) = P_{xy}

が任意の x, y ∈ Ω について成り立つとする．この更新関数があれば，時刻 t での状態が x のとき， [0,1) 上の乱数 Λ を発生させて ϕ(x, Λ) を計算することで，時刻 t + 1 の状態 y = ϕ(x, Λ) を得ることができる^{1}．

[figure: 時刻; t; t + 1; 0.7; s_{1}; s_{2}; s_{3}; s_{1}; s_{2}; s_{3}; \begin{matrix}s_{1} & s_{2} & s_{1} & s_{2} & s_{1} & s_{3} \\ s_{2} & s_{1} & s_{3} & s_{3} & s_{2} & s_{1} \\ s_{3} & s_{3} & s_{2} & s_{2} & s_{1} & s_{3} \end{matrix}; λ :; 0; .2; .4; .6; .8; 1; ∼; ∼; ∼; ∼; ∼]

図 2: 更新関数の絵．左の表は図 1 のマルコフ連鎖を実現する更新関数である．一番左の列は時刻 t における状態であり，表の各列は乱数 λ に対応する時刻 t+1 の状態を表す．右の図では乱数 λ = 0.7 の場合の時刻 t から t + 1 への推移を表す．左の表に照らしあわせると， λ ∈ [0.6, 0.8) なので s_{1} ↦ s_{1}, s_{2} ↦ s_{2}, s_{3} ↦ s_{1} となっていることが確認できる．

更新関数の概念を再帰的に拡張すると，時刻 t_{1} から t_{2} への M の推移も乱数列 λ = (λ[t_{1}], λ[t_{1} + 1], …, λ[t_{2} − 1]) ∈ [0, 1)^{t_{2}−t_{1}} と決定的関数 Φ_{t_{1}}^{t_{2}} : Ω × [0, 1)^{t_{2}−t_{1}} → Ω で表現される．時刻 t_{1}, t_{2} (t_{1} < t_{2}) に対して，乱数 λ を使って時刻 t_{1} で状態 x ∈ Ω から推移した後の時刻 t_{2} の状態を Φ_{t_{1}}^{t_{2}}(x, λ) と表すことにする．すなわち， Φ_{t_{1}}^{t_{2}}(x, λ) =^{def.} ϕ(ϕ(⋯ (ϕ(x, λ[t_{1}]), …, λ[t_{2} − 2]), λ[t_{2} − 1]) である．

CFTP

マルコフ連鎖を本当に無限回推移させれば平衡状態に達するが，それは事実上不可能である．そこで，仮想的に永遠の過去から無限回推移しているマルコフ連鎖について考えよう．すると，そのマルコフ連鎖の現時点（時刻 0）の状態は平衡状態に達している．しかし，そのマルコフ連鎖は初期状態も，推移した回数も分からないので，現在の状態など分かり得ないはずである．

そこで，近い過去の推移から現時点（時刻 0）の状態を推理することを試みよう．このマルコフ連鎖は一様乱数と更新関数で表現されているので，時刻 0 から遡って各時刻での乱数 λ を知ると，「もし時刻 t (t < 0) の状態が x ならば時刻 t + 1 の状態は y = ϕ(x, λ) である」という仮想的な推移が行える．この操作を繰り返すと，「時刻 t (t < 0) の状態が x で，時刻 t から 0 までの乱数列が λ ならば，時刻 0 の状態は y = Φ_{t}^{0}(x, λ) である．」ということがわかる．

乱数の探索を十分過去に遡っていって，「時刻 T (T < 0) の状態がどんな x ∈ Ω であっても，乱数列 λ に対する時刻 0 の状態は同一の y ∈ Ω しかありえない（∃y ∈ Ω, ∀x ∈ Ω, Φ_{T}^{0} (x, λ)）」ということが起これば，それより過去の状態を知らなくても，時刻 t の状態は Ω の要素のいずれかしかありえないので，時刻 0 の状態は y であることがわかる．この状態を収斂（coalsecence）といい，仮想的に無限の時間の推移を経たマルコフ連鎖の現在の状態が y である証拠となる．これこそが CFTP のアイデアなのである！！

アルゴリズム

実際のアルゴリズムとしては，各時刻で一様実数乱数を生成させていくことで，仮想的なマルコフ連鎖の各時刻での実数乱数を得ることができる．このとき，一度生成された乱数は保存されていることに注意してほしい．すなわち，各時刻における乱数を徐々に時間を遡って探索しているのである．

以上のアイデアの下，基本的な CFTP アルゴリズムは以下のように記述される．

アルゴリズム 1 (CFTP アルゴリズム) \underline{Step 1.} シミュレーションの開始時刻を T := −1 とする．空列 λ を用意する． \underline{Step 2.} 乱数 λ[T], を生成し，数列 λ の先頭

1P に対し，この更新関数は唯一ではない．

[figure: 時刻; 0; -1; -2; -3; -4; -5; λ[−5]; λ[−4]; λ[−3]; λ[−2]; λ[−1]; s_{1}; s_{2}; s_{3}; s_{1}; s_{2}; s_{3}; s_{1}; s_{2}; s_{3}; s_{1}; s_{2}; s_{3}; s_{1}; s_{2}; s_{3}; s_{1}; s_{2}; s_{3}]

図 3: CFTP の絵．上の図は，乱数の列 λ = (λ[−5], λ[−4], λ[−3], λ[−2], λ[−1]) = (0.5, 0.1, 0.7, 0.3, 0.5) を用いて，図 2 の更新関数に従って推移するマルコフ連鎖を表す．上の図において，時刻 −3 の状態が s_{1}, s_{2}, s_{3} のいずれであっても，時刻 0 の状態は s_{2} であることが分かる (∀x ∈ Ω, Φ_{−3}^{0}(x, λ) = s_{2} ：図の太線参照) ので，時刻 −4 以前の推移を調べなくて良い．逆に時刻 −2 からの推移を確認すると，収斂していない (Φ_{−2}^{0}(s_{1}, λ) = Φ_{−2}^{0}(s_{2}, λ) = s_{2}, Φ_{−2}^{0}(s_{3}, λ) = s_{3}) ので，時刻 0 の状態が s_{2} と s_{3} のどちらなのかがわからない．

に挿入する．すなわち， λ := (λ[T], λ[T + 1], …, λ[−1]) とする． \underline{Step 3.} Ω のすべての状態について，乱数列 λ を用いて時刻 T から時刻 0 までマルコフ連鎖を推移させる．このとき， a) もし，収斂していれば（すなわち， ∃y ∈ Ω, ∀x ∈ Ω, y = Φ_{T}^{0} (x, λ) ならば）， y を返し，停止する． b) そうでなければ，シミュレーションの開始時刻を T := T − 1 としてステップ 2 に戻る^{2}．

定理 1 (CFTP [4, 1]) CFTP アルゴリズム（アルゴリズム 1）が確率 1 で終了するなら，値を 1 つ返し，得られる値は M の定常分布に厳密に従う確率変数の実現値である． □

CFTP アルゴリズムの反復回数は確率変数である．すなわち，いつ終わるかは確率的にしか分からず，運が（ものすごく）悪いと，永遠に終わらない．このようにいつ終わるか確率的にしか分からないが正確な答えを返すものをラスベガス型の乱択アルゴリズムという．

したがって CFTP アルゴリズムを何らかの理由で途中で止めてはいけない．例えば，停止までに時間がかかり過ぎるなどの理由で（停止した際に得られたはずの）サンプルを捨ててしまうと一般には歪みが生じ，不正確になってしまうからである．

単調 CFTP

こうして無限時間のマルコフ連鎖のシミュレーション結果を（確率的に）有限時間で得ることができた．しかし， CFTP アルゴリズムをそのまま使おうとすると， Step 3. での「Ω の全状態からのシミュレーション」が大きな障壁となる．前にも述べたとおり，マルコフ連鎖のシミュレーションでは，往々にして状態数 N が非常に大きいからである．

そこで次の条件を満たすマルコフ連鎖を考えよう．状態空間 Ω 上の半順序関係 ⪰ が存在し，その最大元 x_{max} および最小元 x_{min} が存在する．また，マルコフ連鎖の更新関数はすべての x ⪰ y (x, y ∈ Ω) に対して， ϕ(x, λ) ⪰ ϕ(y, λ) がどんな λ に対しても成り立つとする．このようなマルコフ連鎖を単調なマルコフ連鎖であるということにする．

単調なマルコフ連鎖であれば以下が成り立つ．半順序つきの状態空間 (Ω, ⪰) では任意の状態 x ∈ Ω に対して， x_{max} ⪰ x ⪰ x_{min} である．単調の定義から，任意の λ に対して ϕ(x_{max}, λ) ⪰ ϕ(x, λ) ⪰ ϕ(x_{min}, λ) が成り立つ．再帰的に任意の t < 0 と任意の数列 λ に対しても Φ_{t}^{0}(x_{max}, λ) ⪰ Φ_{t}^{0}(x, λ) ⪰ Φ_{t}^{0}(x_{min}, λ) が成り立つ．

実は，単調なマルコフ連鎖に対して次の定理が成り立つ．

定理 2 (単調 CFTP [4, 1]) 単調なマルコフ連鎖に対して， CFTP アルゴリズム（アルゴリズム 1）は確率 1 で終了し， ∃y ∈ Ω, ∀x ∈ Ω, y = Φ_{T}^{0} (x, λ) と Φ_{T}^{0} (x_{max}, λ) = Φ_{T}^{0} (x_{min}, λ) は等価である． □

したがって単調マルコフ連鎖に対してはアルゴリズム 1 の Step 3. a) を \underline{Step 3′.} 状態 x_{max}, x_{min} ∈ Ω について，乱数列 λ を用いて時刻 T から時刻 0 までマルコフ連鎖を推移させる．このとき， a′) もし ∃y ∈ Ω, y = Φ_{T}^{0} (x_{max}, λ) = Φ_{T}^{0} (x_{min}, λ) ならば y を返し，停止する．

2反復毎に T は負の小さな数字に更新されることに注意されたい．

で置き換えることができる．すなわち全状態ではなく， 2 状態 x_{max}, x_{min} からのみ推移を行えば良い．このアルゴリズムを単調 CFTP アルゴリズムと呼ぶ．

[figure: 時刻; 0; -1; -2; -3; -4; -5; +1; −1; −1; −1; +1; s_{1}; s_{2}; s_{3}; s_{4}; s_{1}; s_{2}; s_{3}; s_{4}; s_{1}; s_{2}; s_{3}; s_{4}; s_{1}; s_{2}; s_{3}; s_{4}; s_{1}; s_{2}; s_{3}; s_{4}; s_{1}; s_{2}; s_{3}; s_{4}]

図 4: 単調 CFTP の絵． s_{1} ≻ s_{2} ≻ s_{3} ≻ s_{4} として，確率 1/2 でひとつ上の状態へ，確率 1/2 でひとつ下の状態へ推移する（それぞれ動けない場合には自分自身にとどまる）マルコフ連鎖．時刻 −5 における状態 s_{1} （最大元）と s_{4} （最小元）からの収斂（図中の太線）を確認すれば，全状態からの収斂を確認したことになる．

イジングモデル

単調 CFTP の応用例として k × k （直交）格子イジングモデルを紹介しよう．イジングモデルは物質の磁性を温度とスピンから説明する，平衡統計物理では非常に基本的で重要なモデルである [4]．また，このモデルは 2 値（白黒）画像の表現としてもよく知られる．

このモデルでは k × k 格子点上の分子スピンが −1 か +1 かという状態を考える． k^{2} 個の格子点の集合を V で表し， k × k 格子点上の状態 {−1, +1} の割り当て方全体を {−1, +1}^{V} で表すことにする．

イジングモデルに対するマルコフ連鎖 M_{i} は {−1, +1}^{V} を状態空間として持ち，現在の状態 X ∈ {−1, +1}^{V} からの推移は一様実数乱数 Λ ∈ 図 5: 15×15 格子上のイジングモデルの絵（x ∈ {−1, +1}^{V} ）． ∘ は −1， • は +1 を表す．

[0, k^{2}) と更新関数 ϕ : {−1, +1}^{V} × [0, k^{2}) → {−1, +1}^{V} を用いて定義される．表記の簡単のため， X′ = ϕ(X, Λ) で表す．いま， k^{2} 個ある各格子点には 0 から k^{2} − 1 までの番号が振ってあるとしよう．まず， Λ ∈ [0, k^{2}) に対して，その整数部分 ⌊Λ⌋ に対応する頂点 v ∈ V を選ぶ．次に Λ の小数部分を用いて，もし Λ − ⌊Λ⌋ < \frac{e^{2βk_{+}(v,X)}}{e^{2βk_{+}(v,X)}+e^{2βk_{−}(v,X)}} ならば X′ の頂点 v の値を +1，そうでなければ −1 にする．ただし， k_{−}(v, X) および k_{+}(v, X) はそれぞれ，状態 X での v の隣の −1， +1 の個数を表す． v 以外の頂点 w については X′(w) = X(w) とする．このマルコフ連鎖はギブスサンプラーと呼ばれ，定常分布がギブス分布になることを示すことができるが，紙面の都合上，割愛する．

さて，このマルコフ連鎖が単調であることを確認しよう．まず，状態空間 {−1, +1}^{V} 上に半順序関係を定義する．状態の対 x, y ∈

{−1, +1}^{V} に対して x ⪰ y ⇔^{def.} x(v) ≥

y(v), ∀v ∈ V とする．ただし， +1 ≥ −1 とする．この半順序集合 ({−1, +1}^{V} , ⪰) において，明らかに x_{max} ∈ {−1, +1}^{V}, x_{max}(v) =^{def.} +1 (∀v ∈ V) は最大元であり， x_{min} ∈ {−1, +1}^{V}, x_{min}(v) =^{def.} −1 (∀v ∈ V) は最小元である．

次に更新関数が単調であることを確認しよう．すなわち， x ⪰ y ならば任意の λ で ϕ(x, λ) ⪰ ϕ(y, λ) が成り立つことを確認する．このとき，半順序の定義から任意の v ∈ V に対して k_{+}(v, x) ≥ k_{+}(v, y) および k_{−}(v, x) ≤ k_{−}(v, y) であることを考慮すると，

\frac{e^{2βk_{+}(v,x)}}{e^{2βk_{+}(v,x)} + e^{2βk_{−}(v,x)}} ≥ \frac{e^{2βk_{+}(v,y)}}{e^{2βk_{+}(v,y)} + e^{2βk_{−}(v,y)}}

が得られる．したがって， ϕ(y, λ)(v) = +1 ⇒ ϕ(x, λ)(v) = +1 がいえる．同様に ϕ(x, λ)(v) = −1 ⇒ ϕ(y, λ)(v) = −1 もいえる．このことは任意の λ ∈ [0, k^{2}) に対して， ϕ(x, λ)(v) ≥ ϕ(y, λ)(v) (∀v ∈ V) を意味する．したがって， ϕ(x, λ) ⪰ ϕ(y, λ) が成り立つ．

以上のことからイジングモデルのギブスサンプラーは単調であることがわかる．すなわち任意の x ∈ {−1, +1}^{V} に対して，どんな乱数列 λ ∈ [0, k^{2})^{−T} (T < 0) が与えられても， Φ_{T}^{0} (x_{max}, λ) ⪰ Φ_{T}^{0} (x, λ) ⪰ Φ_{T}^{0} (x_{min}, λ) となる．これはある時刻 t (T < t ≤ 0) において，頂点 v ∈ V の値が Φ_{T}^{t} (x_{max}, λ)(v) = Φ_{T}^{t} (x_{min}, λ)(v) ならば ∀x ∈ {−1, +1}^{V}, Φ_{T}^{t} (x, λ)(v) = Φ_{T}^{t} (x_{max}, λ)(v) = Φ_{T}^{t} (x_{min}, λ)(v) である．

おわりに

CFTP によるパーフェクトシミュレーションはセンセーショナルなアルゴリズムとして大きな波紋を与えた．文中でも述べたとおり，このアルゴリズムの厳密性は非常に繊細で，一見良さそうに思える改良が厳密性を崩してしまう恐れがあるので，改良には十分な配慮が必要である [1, 2]．現在，アルゴリズムに対してさまざまな改良が試みられ，いくつかの非常に面白いアルゴリズムがある．興味のある読者へ Wilson の read once アルゴリズム, Fill のアルゴリズム, 連続状態マルコフ連鎖への拡張の 3 つの話題をお勧めする．また， [5] ではパーフェクトシミュレーションの動く様を描いた美しいグラフィックスを見ることができる．

参考文献

[1] X. K. Dimakos, “A guide to exact simulation,” International Statistical Review, 69 (2001), pp. 27–48.

[2] O. Häggström, Finite Markov Chains and Algorithmic Application, London Mathematical Society, Student Texts 52, Cambridge University Press, 2002.

[3] S. Kijima and T. Matsui, “Polynomial Time Perfect Sampling Algorithm for Two-rowed Contingency Tables,” METR 2003-15, Mathematical Engineering Technical Reports, University of Tokyo, 2003. (available from http://www.keisu.t.utokyo.ac.jp/Research/techrep.0.html)

[4] J. Propp and D. Wilson, “Exact sampling with coupled Markov chains and applications to statistical mechanics,” Random Structures and Algorithms, 9 (1996), pp. 232–252.

[5] Web Site for Perfectly Random Sampling with Markov Chains, http://research.microsoft.com/ ˜ dbwilson/exact/