壁にサンプリングしよう

完―第二話l

天 と地の狭間で―釆嶋 秀治，松井C 「 TP(7@落とし穴・知己

シ を有限時間で実現(2)
温理の基本概念である・ただし式 (2)中の定数アについて紹介しよ，度 T(7@時のスピン状態の存在確率を表す平衡統計物

する CFTP アルゴリズムを紹介した・今回は CFTP
アルゴ リズム をより効率よくするためのアイテ

卜(T)= T)e-H(PZ)で表わされる・ この分布は GjbbS 分布と呼ばれ，「前回は，無限のシミュレーショ｜・

Z(T)

つ些
エ三、2LH(( ― l ・― l}')，有限の状態空間 n 上の エ ルゴード的なマルコフ 連は分配関数と呼ばれ，統計物理において非常に重要な鎖を考える・前回紹介した CFTP アルゴリズムは 状値である・ MCMC 法は ，分配関数の計算法として古態 空間 Ql と 更新関数1) づ n で定義されたマルコフ連鎖に対して，過去のある時点における全状くから知られる有効な手段である態からのシミュレーションを実行し，現時点での間の大きさは川である

イジングモデルは

2 値 画像

coaIescence を確認するというものであったの表現として用いられると述べたが，例えば，

図 1のめ: nX[o，(F 対するマルコフ連鎖 ) の状態 空

イジングモデルしかし・・・「会状態からのシミュレーションを実行」することは，

ように 15X15 という小さな画像の 例 ' に対してさえ，実務上の多くの場合，状態数が膨大になるため大変

困

状態空間のサイズはChain MOnte

1 2 2 元分割 表

割難である・ 実際に MCMC(MarkoV・，，，=10"7 となる・・271@
分割 表は， 与えられた行和と 列 和を満たし，CarIo) 法がよく用いられているイジングモデルと分表のランダムサンプ りング を例に挙げて説明しよ つ1 12( 図 2) ・分割表に対する要素が非 負 整数の表である、イジングモデル

各

統計学的興味として，付和と 列 和の相関の検定があげイジング (ISing) モデルは，物質の磁性を温度とられる・正確検定法に必要な 少 値の計算に，マルコスピンから説明する，平衡統計物理では非常に基本的で重要なモデルである( 白黒 )[4]

また，このモデルは

・[7]

画像の表現としてもよく知られる

2値

000
●0000
00000
●●●
000
0oooo@oUOUo@@@ooooO●●●●●●●●
o0oo@@oo●●●●●●●●
O0o@oo0o@@@@OOoO
0 ●●●
0 000
●●
o0Oo@@o・

頂点数乃り グラ 7 G 二 ((V,E) 上のイジングモデルを考えよう

イジングモデルでは，各頂点上の・( ― 1 十 U が与えられた，状態 ガ三 ( ― 1 十 n，，●●●についギ) - 牛。・，ズ (")2(めO
●●
O
O ●●●
●●
● ●●
●
●● ●●
o(1)

0・(ょT( び )G ( ― 1 十 1) は状態 ょ三 ( ― 1，+1) 。における頂点 ひ三 げのスピン f ― l 十 1) を表す，

状態 ェの 存在確率 Pr( ガ) は エネルギ関数に対して，

図丁，，・0 o0

しゅう じ，まっいともみ

東京大学大学院情報理工学系研究科
十 113-8656文京区本郷 7-3-1

264 (46)o o0 0 ●00
●● ●●

●

イ●ジ
十 各子上の。工H X 1 5DLの 糸八 お三 丁ノレ、 ン十 グモ テす哀"）ノを表よO
1[十1は・―?

' 来 : 外部磁場はないのですか
きじま00U00000
00000
0O0OO00000● ● ● ● ● ●●● ● oo●●
0 0 0 0●00oOOUOo@o000
00 ● ●●●●●●●●●● ●● ●●OUU OO ●●●●● ●0て議論する・温度 T ノ O において，状態 T はエネルを持つ，ただし，

●●●

● ●●
●●

0 0 0 ● ●●● ● ● 00
O O 00 O

ス ビン。

● ● ●●●●●@@o0o@@@@OoOooUOUU
0ooooo00 0Ooo@o0ら省略・

松 : 式が煩雑になるか，松 : 丸が画像の ピ クセル で グラフの頂点だ ・来 : 各項
点 はその上下左右の頂点と隣接しているのが普通ですね，© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

オ ベレーションズ・りサーチ

（左），およびそれを満たす分割表（右）の例

フ連鎖を用し−て分割蓑を一様に生成し，モンテカルロ法を用いる手法がよく知られている．しかし，与えら

j≡

さ：

： さ

t t＋1時刻

t t＋1時刻

図3 単調マルコフ連銀

れた周辺和を満たす分割表の個数（すなわち状態空間の大きさ）は一般に膨大となり，さらにそれを数える

三∴∴∴∵図2 与えられた行和と列和ロ［司 巨］巨］□

ことは行数を2に固定した場合でさえ困難3であることが知られている．

2．単調CFTP
それではいよいよ，効率的なCFTPアルゴリズムを設計するための「とっても巧妙な更新関数」を導入しよう．

一5−4 −3 −2−10 時刻

図4 単調マルコフ連鎖のcoalesceの様了一①⑪④④

定義1更新関数が単調であるとは，次の二つの条件つとする．このとき，次の命題が導かれる．

を満たすことである．

套牲4：状態空間n上に半順序関係 ≧が存在す

命題2 更新関数￠で定義されたマルコフ連鎖は半順る．

序集合（n，≧）に関して単調であり，∃LrmaX，∃trmIn∈

垂堕j：任意の推移は半順序を保存する．すなわち，

fl，∀x∈n，．rmaX≧r≧．rmi。とする．この時，T＜0おJ≧牒を満たす任意の状態対（J，〟）∈n2および任意のよび人∈［0，1）Tlに対してcoalesceしている（すなわ実数ス∈［0，1）に対して，￠（J，バ）≧￠（ダ，人）が成り立つ．

ち］〟∈n，∀∬∈n，y＝◎冥J，人））ことの必要十分条
件は◎冥Jmax，ス）＝郵拉mi。，人）である．

単調な更新関数で定義されたマルコフ連鎖のことを，例 図4は図3のマルコフ連鎖を用いて時刻−5から単調マルコフ連鎖と呼ぶ．

時刻0まで推移させた例である．図4では時刻−5に例 図3は単調マルコフ連鎖の例である．このマルコおける全状態からのサンプルパスが時刻0でcoa−

フ連鎖は状態空間（sl，S2，S3，S。）を持ち，半順序（実はlesceしているが，このことは，図の太線，すなわち全順序）をsl＞s2＞s3＞∫4と定義する．マルコフ連鎖

時刻−5における最大元51と最小元∫。からの2本のの推移はコインフリップによって行われ，コインがおサンプルパスを確認するだけでよい．

もて（回）の場合は図3の左の図のように，うら（□）の場合は図3の右の図のように推産する．図3
において，推移を表す矢印が交差することはなく，こ以上の議論から，標準的な単調CFTPアルゴリズ
ムは次のように定義される．

れは推移に対して半順序が保存されることを意味する．
すなわちこのマルコフ連鎖は単調である．

アルゴリズム1（単調CFTPアルゴリズム）
阜土塁旦ユ：シミ

いま， 半順序集合（n，≧）に対して，唯一の最大元ュレーションの開始時刻をr：＝一1とする．空列人＝（）を用意する．

∬max∈nおよび最小元Jmi。∈nが存在するとしよう．

型些j：乱数甘n…，凡打2」＋1］を生成し，数列人

すなわち，任意のJに対してJmax≧J≧Jmi。が成り立の先頭に挿入する．すなわち，人：＝（ス［T］，ス［T3来：＃P完全ですね．松：そう．「ナンバーピー」完全だ−ト1］，…，スト1］）とする．
Step3：最大元J。aXと最小元Lrminについて，共通のね．

2005年4月号© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

（47）265

乱数列スを用いて時刻Tから時刻0までマルコフ連を確認しよう．まず，状態空間卜1，十1）γ上に半順鎖を推移させる．このとき，序関係 ≧ を導入する．状態∬，〝∈卜1，＋1）yに（a）もし，COalesceしていれば（すなわち，∃y∈n，〟＝離抱max，A）＝引抱ml。，ス）），時刻0の状態〝を対して出力し，停止する．とする．ただし，十1≧−1である．この半順序集合（b）そうでなければ，シミュレーションの開始時刻（（−1，＋1）y，≧）において，あきらかに最大元およびをT：＝2TとしてStep2に戻る．

最小元x。aX，．rmi。∈ト1，＋1）Vが存在して，それぞれ∬≧〃笛エ（〃）≧y（〃）（∀〃∈Ⅴ）∬max（〃）＝十1（∀〃∈Ⅴ），∬ml。（〃）＝−1（∀〃∈Ⅴ）を
定理3 エルゴード的な単調マルコフ連鎖に対して，満たす．

アルゴリズム1は確率1で有限停止し，アルゴリズム次に，マルコフ連鎖滅を実現する更新関数￠：

1の返す値は定常分布に厳密に従う．

（−1，＋1）Vx［0，乃）→ト1，＋1）yを次のように定義する．いま，Ⅴの各項点には0から乃−1までの番号が3．単調マルコフ連鎖の例振ってあるとしよう．実数ス∈［0，乃）に対して4，その定理3から，単調なマルコフ連鎖を設計できたなら整数部分しス」の番号を持つ頂点を〃∈Ⅴとする．このば，最大元と最小元の2状態からのみシミュレーショ

とき，次の時刻の状態z＝￠（∬，ス）を，実数スの小数ンを行えば良いので，状態空間の大きさに依らず部分八一しス」を用いて，coalescence時間に比例する時間でパーフェクトサンプリングが行える．前の節で非常に簡単な単調マルコz（〃）＝（フ連鎖の例を示したが，この節ではより実践的な例と−1，（0≦バーしス」＜♪（〃，ズ）），＋1（カ（〃，ガ）≦バーL人」＜1），z（棚）＝エ（棚），ぴ∈Ⅴ＼（〃），

して，イジングモデルと2元2行分割表の一様サンプと定義すると，￠はマルコフ連鎖誹を実現している．

リングに対する単調マルコフ連鎖を紹介しよう．

更新関数￠の単調性を示す．いま∬≧〃として，3．1イジングモデル

実数人∈［0，乃）に対応する推移を考える．表記の簡単のため∬′＝￠（∬，ス）および〟′＝￠（肌ス）と記し，L人」の番号を持つ頂点を〃とする．このとき，更新関数に単純無向グラフC＝（Ⅴ，E）と温度Tに対するイジングモデルを表現するマルコフ連鎖誹を次のように定義する．マルコフ連鎖滅は状態空間卜1，＋1）yを持ち，現在の状態ズ∈卜1，十1）yから次の時刻の状態ガ′への推移は次のように定義される．まず，Cの頂よって〃以外の頂点のスピンは保存されるので，〃以点〃∈Ⅴを一つ，一様ランダムに選ぶ．頂点びに隣外の任意の頂点紺に対して，∬′（紺）≧y′（ぴ）が成り立っ．次に，頂点〃について々＿（〃，〝）≧々一（〝，∬）を考慮
すると，カ（〃，∬）≧♪（〃，〝）が成り立つ．このことは更接する頂点のうち，状態Xにおいてスピンが＋1の新関数の定義から，頂点の個数を々．（〃，ズ）で表し，−1の頂点の個数を∬′（〃）＝−1⇒ 〟′（〃）＝−1，々＿（ぴ，ズ）で表す．また，〝′（〃）＝＋1⇒J′（〃）＝＋1如，X）竺

が成り立つことを意味する．したがってJ′（〃）≧〝′（〃）

e与棉ズ）
（3）

e与烏−（仇ズ）＋e与々十（む・ズ）となる．とする．このとき，確率〆〃，X）でズ′（〃）＝−1とし，確率1一♪（〃，X）でX′（〃）＝十1とする．その他の頂定理4 更新関数￠は単調である．すなわち，もし∬点緋∈Ⅴ＼（〃）のスピンについてはズ′（抑）＝ズ（紺）と≧〟ならば任意の実数ス∈［0，乃）に対して￠（J，ス）≧する．有限マルコフ連鎖誹は既約で非周期であるので，￠（肌ス）が成り立つ5．

エルゴード的である．また詳細均衡方程式が成り立つ
ことから，定常分布はGibbs分布（2）となる．このマルコフ連鎖J4tはM−H鎖（Metropolice−Hastingschain）と呼ばれる．

−1の2状態からのみシミュレーションすればよい．
4桧：一様実数乱数のドメインを【0，1）から［0，乃）に変換

上記で定義されたマルコフ連鎖滅が単調であること

266（48）

上記より，アルゴ］）ズム1ではxmax＝1とxml。＝したところに注意ね．来：メモしときます．
5桧：外部磁場があっても単調性の証明は同じだよね．
来：式が煩雑になるけどいいんですか？松：それはイヤ．
オペレーションズ・リサーチ© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

（4）

図5 周辺和（左上），および推移の例）

UL（⊃〃丁ズ／×V八和◇方／ 和円／／和以和侃〝ノ10 1 2 3 4 5 6 t
図6 累積和ベクトルの例

3．2 2元2行分割表なのでエルゴード的である．また，詳細均衡方程式を行和ベクトルr＝（rl，巧）および列和ベクトルβ＝満たすので，マルコフ連鎖の定常分布は一様分布である．

（∫．，…，S乃）を満たす2×乃分割表の集合を≡とする．

状態空間≡上に半順序関係を導人しよう．任意の

ただし，rl十r2＝51十‥・＋s乃＝Ⅳである（図5）．以降，≡1二の一様分布を定常分布に持つマルコフ連鎖ノ彫を定ズ∈冨に対して，累積和ベクトルム聖（ム（0），義する．

ム（1），…，ム（〃））を分割表X∈≡の汗り列の値をズ［g，ノ］で表す．状（g＝0），0郁）聖‡態空間≡を持つマルコフ連鎖誹の推移は次のように定義される．まず隣り合う2列（ノ，ノ＋1）（ノ∈（1，．‥，〃

ガ［1，1］＋…＋ガ［1，〜］（g∈（1，．．．，乃）），と定義する．このムはズの一意表現になっている．

−1））を一様に選ぶ．次に，ノ列目およびノ＋1列日に任意の対X，y∈≡に対して，X≧yの必要十分条件をムー／y≧0とする．2項関係 ≧が＝上の半順序

十々 一々 一々 十々（5）であることは明らかである．図6の折れ線は累積和ベクトル（の線形汎化関数）である．半順序関係の大小の形の部分表を加えたとき各要素が非負となる表に一様に推移する．このとき，整数烏は負の値も取り得関係は折れ線表現での上下関係に対応する．
この半順序集合における最大元と最′ト元について考ることに注意されたい．図5は推移の例である．状態える．二つの特別な分割表風，丸∈≡を北西隅の規方に対して（3，4）列が選ばれたとき，非負性を満たす則および北東隅の規則で得られる表とする．与えられ表はH，坑，坑の三つであり，それぞれに確率1／3た半順序集合（≡，≧）におし、て，任意の方∈≡に対しで推移する．マルコフ連鎖｝瀦は有限で既約で非周期的て，晶≧方≧丸が成り立つ．このことは，図6にお2005年4日圭© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

（49）26丁いてすべての分割表X∈≡がぶ」と＆の間に挟まれることを表す．

この定理はマルコ7連鎖のmixing時間の算定に用いられるcoupling法を利用して示すことができる．

このように半順序を定義すると，マルコフ連鎖滅に

定理5はマルコフ連鎖ノ紺に対する任意の更新関数（に対して単調な更新関数を定義できる．詳細は文献［1］おけるcoalescence時間）について成り立つことに注を参照いただきたい．

意が必要である．

4．coa］escence時間とmixing時間
以上の議論から，COalescence時間の短い単調マル定理5より，マルコフ連鎖に対して適当な更新関数を設計して，そのcoalescence時間を知ることができればマルコフ連鎖のmixing時間すなわち収束の速さコフ連鎖を設計することで，計算効率の良いパーフェの上限を押さえることができる．例えば単調な更新関クトサンプリング法が得られることが分かる．一方，数が設計できれれば，COalescence時間の期待値を計

通常MCMC法で用いられている近似サンプリングで算機で推定することはそれほど難しくない．

も，マルコフ連鎖の収束の速さは計算効率の重要な鍵

さらに次ではマルコフ連鎖ノ紺に対する更新関数￠となる．実はcoalescence時間とマルコフ連鎖の収束は半順序≧に対して単調であり，また半順序集合の速さの間には密接な関係がある．次では，この関係（n，≧）は最大元∬maxと最小元∬m．。を持つ場合を考えについて議論しよう．

よう．ここで，∬maxと∬ml。を結ぶ（≧に関する）極以降，有限の状態空間nを持ち，推移確率行列P大鎖の長さをβとする．このとき，定理5とは逆におよび定常分布方を持つエルゴード的なマルコフ連

COalescence時間の期待値がmixing時間によって上鎖．〝を考える．準備としてマルコフ連鎖の収束スピーから押さえられる．

ドを評価するための尺度を導入しよう．（同一の）有限状態空間n上の二つの確率分布レ1とレ2が与えられ

定理6 更新関数￠を用いたアルゴリズム1のcoa−た時，レ1とレ2の間の総分布距離は1escence時間71はdTV（レ1，レ2）竺‡黒レ庸一レ2（∬）l現れ］≦21n（2e（β−1））丁（‡）と定義される．マルコフ連鎖滅のmixing時間6は，を満たす．

任意の正数E＜1に対してT（E）竺max（min（tl∀s≧t，dTV（7r，Pf）≦E））
J∈nと定義される7．ただし，だは時刻0での状態J∈nとした際の時刻sでの状態の確率分布とする．mix−ing時間はマルコフ連鎖に対して固有の定数であり，更新関数には依存しないことに注意されたい．
coalescence時間は更新関数に依存する確率変数ですなわち，マルコフ連鎖の収束がある程度の精度（E≦1／e）で要求される場面において，もし単調な更新関数が設計できれば，定理5と定理6からmixing時間の上下限がcoalescence時間の期待値によってE［n］≦r（‡）≦r（ど）≦E［れk−121n（2e（β−1））のように押さえられるのである．

ある．関数￠はマルコフ連鎖一瀦を実現する更新関数であるとし，確率変数二nを更新関数￠に対するcoalescence時間とする．

5．まとめ
連載第2回はパーフェクトサンプリングを実現する「単調CFTP」を紹介した．これを用いると，単調で定理5 マルコ7連鎖一瀦を実現する更新関数￠の収束の速い更新関数を設計できれば，COalescence時coalescence時間の期待値E［71］に対して，マルコフ間に比例する時間で完壁なサンプリングができる．ま連鎖．4tのmixing時間T（E）はた節4の議論から，この時間は信頼のおけるサンプリ

0＜∀E＜1，丁（ど）≦E［71］ど￣1ングを行うために最低限必要な時間と同程度であるこを満たす．とも分かった．
今回紹介したもの以外の単調マルコフ連鎖として，

6釆：訳は撹拝時間でしょうか？ 松：読めるけど書けん．

タイリングのランダム生成に対するものがある．

7松：つまりmixing時間は「初期状態に寄らず，総分布Proppのホームページ中にある文献［7］では，この単距離がど以下になるために必要な時間」っていう定義だね．
261（50）オペレーションズ・リサーチ© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

調CFTPアルゴリズムの動く様子が見られる．しかムは効率的になった．しかし，いぎ単調CFTPアル

し単調マルコフ連鎖の設計は容易ではなく，あまり多ゴリズムを実装してみようとすると，「乱数の列を保くは知られていない．

存する領域」の確保が困難であることが分かる．次回

節2で紹介したオーソドックスな単調CFTP以外にもCFTPアルゴ））ズムが効率的に実行できるもの

はこの点を克服する「絶妙のアルゴリズム」を紹介しよう．
（「第三話 終わりある未来」に続く）が提案されている．以降，二つの例を簡単に紹介しよ

つ．

参考文献

5．1単調CFTP以外の例
Propp and Wilsonは，連結無向グラフ上の（根付き）全張木のランダムサンプリングに対して，多項式時間のパーフェクトサンプリング法を提案している［2］．このCFTPアルゴリズムでは，マルコ7連鎖の「最新の推移」が現在の状態に保存されるという性質を用いて，過去からのシミュレーションをすることなく現在の状態を一意に特定する．実はこの間題は，グ［1］S．KijimaandT．Matsui： Polynomialrtimeperfect
sampling algorithm for twoてOWed contingency

tables， METR2003−15，（2003）．
［2］J．Propp and D．Wilson：

Exact sampling withcoupledMarkovchainsandapplicationstostatisticalmechanics，

Random Structures and Algorithms，9（1996），232＝252．

［3］］．Propp and D．Wilson： How to get a perfect

ラフ上の全張木の数え上げを利用した多項式時間の厳

random sample from a generic Markov chain and

密な一様サンプリング法が知られているが，これに比generatearandomspannlngtreeOfdirectedgraph，

べてCFTPアルゴリズムはおどろくほどシンプルで

JournalofAlgorithms，27（1998），170r217．
［4］伊庭幸人： 統計学者・数理I二学者のための統計物理ある．

入門一改訂版（3．02版）一桁f一スピン模型とマルコフ

小沢，高橋，高矯は待ち行列ネットワークに対して連鎖モンテカルロ法を中心にして−， ISM Reseach

定常分布からのパーフェクトサンプリング法を提案しMemorandum，No．635，February，1997．ている［5］．彼らのアルゴリズムは広い意味で単調［5］小沢利久，高橋成風高橋幸雄： 吼洞なマルコフ連銀とネットワークモデルヘの応札 ト1本オペレーションズ

CFTPアルゴリズムであるが，リサーチ′学会2004年春季研究発表会アブストラクト集

・マルコフ連鎖が非エルゴード的である，・最大元からのシミュレーションのみでcoalesceを確認できる，（2004），104−105．
［6］間瀬茂，武田純： 空間データモデリングー空間統汁、γ‥の応用， データサイエンス・シリーズ7，共、Ⅲ1版，2001．という二つの大きな特徴を持つ．

［7］http：／／www．math．wisc．edu／propp／tiling／www／

5．2 次回予告

applets／

単調CFTPを導入することによって，アルゴリズ

2005年4月号© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

（51）269