完壁にサンプリングしよう！
一第一話 遥かなる過去から−
来場 秀治，桧井 知己
l………ll……l…lllll…………l州Il…………‖………l州Il……………………………………ll……11……l……l…l………l……ll……l…………l……………l州Il……ll………………llll…………l……ll…………‖………ll……l川II

MCMC法の利点は状態空間が非常に大きかったり，1．はじめにあるいはその状態数を知ることさえ困難な場合に，局

本連載は3回にわたって「マルコフ連鎖の極限分布所的な情報を手がかりに，集合全体の持つ情報を計算の実現」をテーマとする．この話題は，MCMCすることができる点である．例えば，前者の例として，（MarkovchainMonteCarlo）法との関わりが深い．

イジングモデルの分配関数（正規化定数）の計算が，連載はアルゴリズムの話が中心だが，その奇抜なアイ

また後者の例として，分割表の正確検定が，それぞれデアを楽しんでいただきたい．初回は「無限のシミュ挙げられる．広い意味でのMCMC法は，いまや統計レーションを有限時間で実現」してしまう，夢のよう物理，医療統計，遺伝情報学，数値積分，画像処理，なアルゴリズムの基本アイデアを紹介しよう．

最適化算法など，様々な分野に現れる．
2．3 定常分布の設計

2．マルコフ連鎖の設計

MCMC法を用いる上で，目標の分布を定常分布に2．1マルコフ連鎖持つマルコフ連鎖の設計は重要である．これに関して，以降，状態数が有限のマルコフ連鎖について考える．

目標の分布を設計する方法が知られている．代表的な

いま，マルコフ連鎖誹は∽個の状態で構成される状方法として，M−H（Metropolice−Hastings）法や

態空間n＝（∫1，S2，．．リSm）と推移確率行列Pを持ってGibbsサンプラが挙げられるが，これらが目標の定常いるとしよう．すなわち，誹が各時刻で取りうる状分布を持つことは，次の定理によって保証される．

態はn上のいずれかであり，時刻fで状態J∈nにいる時，時刻f＋1の状態が〟∈nである確率はP（J，定理1エルゴード的マルコフ連鎖誹は有限の状態ダ）である．ただし，P（∬，〝）は行列PのJ，〟成分を空間nと推移確率行列Pを持つとする．状態空間n表す．

上の確率分布万（∑J∈。打（∬）＝1，方（∬）≧0（∀J∈n））マルコフ連鎖誹の状態空間n上の確率分布ガ＝が，任意の状態対（∬，〝）⊂nについて，（打（sl），汀（s2），…，方（sm））について考える．このとき冗方（J）P（J，y）＝打（y）P（y，J）（1）

は∑J∈。方（∬）＝1および方（J）≧0（∀J∈n）を満たす．

を満たす時，確率分布冗はマルコフ連鎖．粛の唯一のこの方が打P＝冗を満たす時，定常分布と呼ぶ．今定常分布である．

後，マルコフ連鎖滅は唯一の定常分布冗を持つ（エルゴード的）と仮定する．このとき，マルコフ連鎖を定理中の式（1）は詳細均衡方程式と呼ばれる．例えば図

無限回推移させると，その状態は定常分布冗に収束

1のマルコフ連鎖は任意の状態対（J，〟）⊂（sl，S2，53）にすることが知られている［1］．
2．2 MCMC法

MCMC法の基本的なアイデアは非常に簡明で，「マルコフ連鎖の極限分布からサンプリングを行い，モンテカルロ法（等）を行う．」というものである．
きじま しゅうじ，まつい ともみ
束京大学大学院情報理工学系研究科〒113−8656文京区本郷7−3−1
2005年3月号

図1マルコフ連鎖の推移グラフ© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

（31）169

人＝0〜喜〜呈〜蔓〜呈〜亘〜1

6

6

6

6

さいころの目‥ 口【司［司［］巨］匡亡

f＋1時刻

図2 更新関数の例
ついてP（J，〟）＝P（〟，J）が成り立つので，定常分布

フ連鎖を何回推移させても，有限回の推移で打ち切っは一様分布になる．詳細均衡方程式は目標の定常分布てしまうとひずみは生じる．

を持つための十分条件であり，マルコフ連鎖の定常分

これに対し，1996年にProppandWilsonは「無限布は必ずしもこれを満たす必要はない．しかし，回のシミュレーションを有限時間で実現」する乱択ア

MCMC法の利点が局所的情報のみで推移できる点にルゴリズムを提案した．このアルゴリズムはCFTPあることを考慮すれば，この局所的な方程式はマルコ（couplingfromthepast：過去からのカップリング2）

7連鎖の設計に非常に有効である．と呼ばれる．それでは，いよいよCFTPについて紹2．4 更新関数介しよう．

マルコフ連鎖の推移は，一様実数乱数A∈［0，1）と，3．過去からのカップリング（決定的）関数￠：nX［0，1）→nをもちいて表現することができる．関数￠は，一様実数乱数A∈［0，1）と3．1アルゴリズム

任意のJ，y∈nに対して，Pr［￠（∬，A）＝〟］＝P（∬，y）

状態空間nと更新関数￠を持つマルコフ連鎖滅にを満たさなければならない．この関数￠は更新関数

対して，CFTPアルゴリズムは次のように定義され（updatefunction）と呼ばれる．

る．

図2の左側の表は，図1のマルコフ連鎖に対する更新関数の例である．簡単のため，乱数の値が区間［（オアルゴリズム1（CFTPアルゴリズム）ー1）／6，才／6）（オ∈（1，2，…，6））に含まれる事象をさいこ塾§旦ユ：シミュレーションの開始時刻をr：＝−1とろの目口だと思うとよい．図2の右の図はさいころのする．空列人：＝（）を用意する．

目が匝］であったときの時刻Jから′＋1への推移を図堅塁旦＿≧：乱数バ［r］を生成し，数列人の先頭に挿入す示したものである．

る．すなわち人：＝（バ［r］，ス［r＋1］，…，スト1］）とする．

更新関数の概念を再帰的に拡張すると，時刻′1か

迦呈：nの各状態∬について，∬を初期状態とし，らゐ（ム＜ち）への推移も乱数列人＝（ス【ム］，バレ1＋1］，…，共通の乱数列Aを用いて時刻rから時刻0までマルス［ね−1］）と下記のように定義される関数◎き…：nX［0，コフ連鎖を推移させる．このとき，1）と2−亡1→nを用いて（a）もし，時刻0でcoa（esce3していれば（∃y∈fl，◎拇，ス）竺両（…如，舶］），…，鵬−2］），鵬−1］）∀∬∈n，y＝◎冥∬，ス）），時刻0の状態yを出力し，停止する．と表すことができる．

もちろん，状態空間nと推移確率行列Pを持つマ
ルコフ連鎖の更新関数の決め方は一意ではない1．

（b）そうでなければ，シミュレーションの開始時刻をT：＝T−1としてStep2に戻る．

2．5 問題点は何か？

目標の定常分布を持つマルコフ連鎖を手に入れた後アルゴリズム1が終了したときのけlの値をcoa・の問題点は，「何回推移させれば定常分布（極限分布）lescence時間と呼ぶ．以降，図3を例にアルゴリズと言えるのか？」ということである．しかし，マルコムを説明する．
2松：この名前長くない？釆：僕は慣れちゃいました．

1松：うまい更新関数が必要になるんだよね．釆：研究者の腕の見せどころですね．
叩0（32）

3松：この訳語は合流でどうかな．釆：収敷くらいが無難じゃないですか？松：「収欽」って漢字書けないだよぉ．

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

オペレーションズ・リサーチ

1反復目

3反復目

：Step2でさいころの目且を得る．Step3でcoalesceを確認するが，やはりcoalesceしていないので（b）に従う．

4反復目：Step2でさいころの目且を得る．Step3でcoalesceを確認する．このとき，時刻0の状態はs3でcoalesceしているので，状態s3を出力しアルゴリ

ー5 −4 −3 −2 −1 0 時刻ズムを終了する．

2反復目

定理2 任意のエルゴード的な有限マルコフ連鎖に対して，ある更新関数4が存在して，アルゴリズム1は［司 巨］確率1で有限停止し，アルゴリズムの返す値は定常分布に厳密に従う．

上記の定理から分かるように，アルゴリズム1によっー5 −4 −3 −2 −1 0 時刻て，定常分布に厳密に従うサンプ1）ング（perfectsampling5）を行うことができる．次節では，このアルゴリズムが何をしているのかを解説しよう．

3反復目

3．2 寓話的解説−マルコフの洞窟

回 回 巨］とある村には「マルコフの洞窟」と呼ばれる洞窟があった．「マルコフの洞窟」の奥ではマルコフ仙人がさいころをふってマルコフ連鎖を推移させている（図2参照）．このマルコフ連鎖はいつ始まって，初期状−5 −4 −3 −2−10 時刻

態が何だったのか，何回推移しているのか，そもそもマルコフ仙人は何のためにこんなことをしているのか，村の誰も知らない．とにかく仙人はこの村ができるよ4反復目

回＋巨］［司＋回りもず一つと苦からさいころをふり続けているので，マルコフ連鎖は定常分布に収束している．だから時刻を決めて（時刻0と呼ぶ），その状態をサンプルすれば定常分布からのサンプリングができるのだが，仙人は少々へそ曲がりで，マルコフ連鎖の「状態」を教えー5 −4 −3 −2 −1 0 時刻てくれない．

ある時，村の若者はこう考えた．仙人の教えてくれ

図3 アルゴリズムの例

るさいころの目から，マルコフ連鎖の状態を推理する
1反復目：Step2でさいころをふり且の目を得る．
Step3ですべての状態sl，S2，S3について，さいころのことはできないだろうか？

若者は仙人に尋ね，時刻−3から時刻0の間に出た目回を用いて時刻一1から0までに対する推移を行う．

さいころの目が匡］，且，［召であることを聞き出した．

このとき，時刻0でcoalesceしていない（状態s2と若者はさっそく推理を始めた．もし時刻−3の状態がs3が存在する）ので（b）に従う．

slだった場合，さいころの目回，瓦，回にしたがって2反復目：Step2でさいころの目瓦を得る．Step3で推移すると時刻0の状態はs2である．同様に，時刻すべての状態について，さいころの目［司，［召を用いて一3の状態がs2だとすれば，時刻0の状態は53，時刻時刻−2から0まで推移を行う．特に，時刻一1から0への推移は1反復目で得られた且を用いているので，4松：定理2の証明は，存在性の証明だけなんだよね．

5桧：この訳語って，完璧サンプリングでいいかな？

1反復目と同一の推移であることに注意する．やはり来：理想サンプリングもありですね．「完壁」は書けるん

時刻0でcoalesceしていないので（b）に従う．

ですか？

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

（33）1Tl

−3の状態が∫3なら時刻0の状態はs2である．だかる）と，厳密性は失われる．

ら，少なくとも現在の状態が∫1でないということは4．いくつかの注意点分かる（図3「3反復目」参照）．しかし，真実はただ一つ．推理にはさらに情報が必要だ．

4．1うまい更新関数

若者はさらに仙人に目をたずねた．仙人はその前の定理2において，確率1でアルゴリズム1が有限停目（時刻−4の目）が回であると告げた．若者は，時止するには うまい更新関数 を考える必要があった．

刻−4の状態が51，S2，53だった場合それぞれについて，実は，この うまい更新関数というのは，COalesce

さいころの日向，且，直し圧で起こりうる時刻0の状する確率が0でないものを指す．もしcoalesceする態を推理した．その結果，時刻−4の状態が∫1，ざ2，ざ3

確率が0の更新関数を使ってしまうと，アルゴリズムの場合のいずれであろうと，時刻0の状態は∫3しか1は停止しない．

ないという結論に達したのである（図3「4反復目」参照）．

そのようなcoalesceしない更新関数の例を挙げてみよう．図4は状態空間（51，52）を持つマルコフ連鎖

こうして若者は思惑通り時刻0の状態を知り，すなわち定常分布からのサンプルを得たのである．

で推移確率は図4の右上のように定義される．このマルコフ連鎖に対して図4左上の表のような更新関数を

しかし，若者はこの一つのサンプルだけでは満足で考える．この更新関数に対して，乱数［0，0．5）の時きなかった．若者はもっと多くのサンプルを必要としの推移は∫1ト◆Sl，∫2ト◆ざ2となり，乱数［0．5，1）の時ていたのだ．残念ながら，マルコフ仙人はたった一人の推移はslト◆S2，S2ト・Slとなり，COalesceしないこしかいない．とが分かる．

ふと，若者は気づいた．マルコフ仙人の教えてくれ

実は，任意のマルコフ連鎖に対して うまい更新関るのは各時刻のさいころの目だけなのだ．その日とい数 が存在することがいえる8．これは状態空間の有うのも，仙人はさいころをふっているだけである．さ限性を考えると比較的容易に示せる．

いころの目を知るだけだったら，何も仙人に聞かなく4．2 なぜ前向きではいけないのか？

たって自分でさいころをふったって一緒のことではな

CFTPアルゴ］）ズムの鍵はcoalesceにあるといえいか！そう，さいころを自分でふって，現時点からる．しかし，ならば，アルゴリズム1のようなまどろ過去に向かってその日を決めていけば，仙人なんて要つこしいことをしなくても，通常のマルコフ連鎖シミらないのだ！

ュレーション同様，「時刻0からの推移を行い，COa−3．3 もう少し数学的な説明

1esceした時刻での状態を返せばよいのではないの

マルコフ連鎖の極限分布からのサンプリングは，片か？」という疑問を持つ方もいるかもしれない．この側無限の乱数列をサンプリングし，対応するマルコフようなアルゴリズムを 前向きアルゴリズム

連鎖の最終状態についての同値類を判定することに他う．次に 前向きアルゴリズム

ならない6

サンプリングが行えないことが簡単に理解できる例とと呼ぼでは，パーフェクト

ここでアルゴリズムのcoalesceの確認とは何かを考えよう．ある時刻での状態に依存することなく最終

入：0 〜．5 〜1

状態が一意に定まっているということは，無限長の乱数列のある時刻より後ろを見て同値類判定を行うことに対応する．このとき，乱数列がマイナス方向に無限であることが肝である．実は， うまい更新関数に対して最終的にcoalesceしている確率はほぼ1とな

0．3（亘）−−−・（室）ることがいえるので，COalesceが確認がされた条件下での同値類の生起確率は定常分布と一致する7．
したがってcoalesceしないからといって，決して中断してはならない．もし中断してしまう（ようにす6松：厳密には ど−Ⅳ の議論が必要になるね

7来：ちょっと難しくないですか？松：難しいねぇ…
1丁2（34）（画一−一一（室）上 土＋1時刻f 亡＋1時刻

図4 coalesceしない更新関数の例
8松： うまい更新関数が存在するというより，COa−1esceしない更新関数が間抜けなんだよね．来：…

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

オペレーションズ・リサーチ

入：0 〜．25〜．5〜．75〜1［司 巨］［司 匝］［召−5 −4 −3 −2 −1 0 時刻

事象：Ⅰ ⅠⅠⅠⅠⅠⅠⅤ
ⅠⅠⅠⅠ

ⅠⅠⅠ

図6 図3の3反復目の後，時刻−5まで遡った例

IVう．
いま，時刻r＊＜0からのシミュレーションを行ったとき，時刻0でcoalesceして状態y∈flであることが確認されたとする．さらに時刻r＊−1の乱数を生成し，時刻r＊−1からシミュレーションを行う．
012

43

5 時刻

図5 前向きではいけない例このとき，時刻r＊の状態は，もちろんnのうちのいずれかなので，時刻0ではやはり同一の〝にcoa−して図5の例を紹介する．この例では状態空間（51，S2，lesceしている．図6は，図3の4反複目の後，さら53，S。）を持つマルコフ連鎖を考える．このマルコフ連に時刻−5に遡って振ったさいころの目が［司だった例

鎖はP（J，ダ）＝1／4（∀（∬，〝）∈n2）であり，定常分布はである．時刻0の状態は図3の4反複目と同じs3で万＝（1／4，1／4，1／4，1／4）である．

ある．

このマルコフ連鎖に対して図5の左上のような更新

つまり，シミュレーションの開始時刻rを時刻関数を与える．もし，時刻0で事象ⅠⅤが起これば当r＊より過去の時刻に設定すれば，時刻T＊からの然coalesceするので，この場合 前向きアルゴリズ

coalesceで得られる状態が得られるのであって，アム，はs4を返す．それ以外の場合，いずれにしても時ルゴ1）ズム1のStep3におけるcoalesceの確認を毎刻1の状態としてsl，S2，∫3を取りうる．そして時刻1反復行う必要はない．そこで次のように改良した標準での事象がⅠⅤなら，やはりs。を返す．それ以外の場CFTPアルゴリズムを作ることができる．カギとな合，時刻2の状態としてはやはリs．，S2，S3のいずれもるアイデアは，Step3（b）で過去へ戻るステップ幅を2取りうる．同様に考えると，この更新関数を用いて倍2倍していくことである．

前向きアルゴリズム を行うと，得られる値は常にアルゴリズム2（標準CFTPアルゴリズム）
姐：シミュレーションの開始時刻をr‥＝−1とs4となってしまうのである．

注意として，この更新関数を使っても，CFTPアする．空列人：＝（）を用意する．

ルゴリズムで得られる値はもちろんパーフェクトサン

阜垣良三：乱数封r］，…，刈LT／2」＋1］を生成し，数列人プリングである．この例は極端ではあるが，一般にの先頭に挿入する．すなわち，人：＝（ス［r］，ス［T前向きアルゴリズム では偏りが起こるので，パー＋1］，…，スト1］）とする．

フェクトサンプリングとしては決してやってはならな

呈堕旦＿蔓：nのすべての状態について，共通の乱数列人い．

を用いて時刻rから時刻0までマルコフ連鎖を推移
4．3 簡単な改良させる．このとき，いま，アルゴリズム1のcoalescence時間を確率変（a）もし，COalesceしていれば（すなわち，∃y∈n，数T＊で表す．各反復でのStep3のシミュレーショ

∀∬∈n，〝＝◎XJ，ス）），時刻0の状態〟を出力し，停ン時間を考えれば，アルゴリズムの計算時間は止する．
（b）そうでなければ，シミュレーションの開始時刻を1＋2＋…＋け＊l＝lT＊l（lr＊けり

2T：＝2TとしてStep2に戻る．に比例することが分かる．

これまでの議論から，ここで，アルゴリズム1でcoalesceを確認しても停止せず，さらに遡ると何が起こるのかを考えてみよ2005年3月号

このアルゴリズムももちろんパーフェクトサンプリングを実現する．このアルゴリ
ズムの計算時間は1＋2＋4＋…＋2「伽2r＊ト1＜4け＊lに比© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

（35）1丁3例するので，COalescence時間のたかだか4倍の時間のシミュレーションで時刻0の状態は得られる．

実用の面において，マルコフ連鎖の収束の速さは大きな問題といえよう．近年，組合せアルゴリズム理論の研究者を中心に様々なマルコフ連鎖の収束の速さの5．まとめにかえて算定に関して精力的な研究が行われている．

連載第1回はパーフェクトサンプリングを実現する「過去からのかソブリング」を紹介した．マルコフ連

5．2 次回予告
CFTPは確かに面白いアイデアである．しかし，鎖の収束に関する研究の中で，1996年のPropp andマルコフ連鎖を用いたサンプリングの最大の利点は，WilsonのCFTPアルゴリズムには多くの研究者が驚全体を知らなくても実行できる点にある．通常，

き，強い関心を持った．英語の紹介記事では amaz＿

MCMC法が状態空間の非常に大きな対象に対して用ing，いられることを考えると，全状態からのcoalesceをsurprizing と紹介されることも多い．以降，簡単にこの分野の研究の動向について述べる．

確認とし−うのはナンセンス極まりない．この点を克服

5．1研究の動向する「とっても巧妙な更新関数」について，次回お話マルコフ連鎖，あるいはランダムウォークは確率モしよう．

デルとして非常にすぐれ，モデル解析，シミュレーシ（「第二話 天と地の狭間で」に続く）ョンに活躍する．さらに近年では，便利な「道具」と参考文献しても注目されている．

統計物理の分野において，MCMC法は20世紀半ばからすでに活躍していた．これをシミュレーションと捉えるか，物理現象を記述するためのモデルと捉えるかは各人に任せるとして，近年，MCMC法を「道具」と捉えた研究が進んでいるのは確かである．一つには計算機能力の向上により大規模な問題に対する計算が可能となり，様々な分野でMCMC法が実用的に用いられる場面が増えてきたためである．

1丁ヰ（36）［1］0．Haggstr6m：
ithmic Application，FiniteMarkovChainsandAlgor−London MathematicalSociety，StudentTexts，52，CambridgeUniversityPress，2002．

［2］J．Propp and D．Wilson‥

Exact sampling with

COupledMarkovchainsandapplicationstostatisticalmechanics，

Random Structures and Algorithms，9（1996），223−252．

［3］大森裕浩： マルコ7連鎖モンテカルロ法の最近の展
開， 日本統計学会誌，31（2001），305−344．

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

オペレーションズ・リサーチ