完壁にサンプリングしよう

一第三話 終わりある未来−
来嶋 秀治，松井 知己
Ill‖‖‖‖‖＝‖‖‖‖‖‖‖＝‖‖‖＝‖‖‖＝‖‖‖‖‖＝‖＝‖‖＝＝＝‖＝‖＝‖‖＝‖‖‖‖‖‖＝‖‖‖‖‖‖＝‖‖‖＝‖‖‖＝‖‖＝＝州…ll…………lll……ll……l…lll……ll…………ll…………ll……l……lll…ll…………lll‖‖‖‖‖‖‖‖＝‖‖‖＝‖‖＝‖‖‖＝‖州

1．消せない記憶
前回は，単調CFTPアルゴ））ズムを紹介した．これにより，高速なパーフェクトサンプリングの設計が可能となった．しかし，いぎ（単調）CFTPアルゴ
リズムを実装しようとすると，思わぬ問題点が存在する．それは乱数の記憶領域である．
刻

CFTPアルゴ1）ズムではcoalesceを確認すること

刻

図1時刻Jlから時刻J2までの推移の模式図によって，パーフェクトサンプリングを実現する．

coalesceが確認できなかった場合には，更に過去に遡ってシミュレーションを行わなければならないが，図1は′1からねの間の推移を模式的に表したもの

このとき，生成された乱数を記憶しておかなければなである．図1の左側は，全状態について数列人5を用らない．例えば，単調CFTPアルゴリズムならば，いて時刻Jlから時刻J2まで推移を行い，時刻J2で計算時間に比例する記憶領域が要求されるのである．

coalesceしている様子を表している．右側は，数列

Wilsonはこの点について改良したReadOnceアル＾fを用いると時刻t2ではcoalesceしていない2様子

ゴ1）ズムを提案している［4］．今回はそのRead Onceである．図1は単調マルコフ連鎖の最大元と最小元かアルゴリズムを説明しよう．

らの推移を表していると思って見るとイメージしやすい．ただし次の議論では，計算量の議論の一部を除い2．準備て，更新関数の単調性3の仮定は必要としない．

有限の状態空間nを持ち，更新関数￠：nX［0，1）→nで推移を定義されたエルゴード的なマルコフ連

Read Onceアルゴリズムを紹介する前に二つの準備を行う．

2．1a遡るCFTP鎖誹を考える．

CFTPアルゴ1）ズムの重要なアイデアは，「COa−まず，特殊なCFTPアルゴリズムを導入する．

1esceを確認する」ことである．これまでは時刻t＜0から時刻0の間の推移だけを考えていたが，少し一般アルゴリズム1（α遡るCFTPアルゴリズム）化して，時刻Jlから時刻わ（ムく′2）の間の推移を考え

塑里軋旦：正整数α＞0を決める．

る．時刻′1から時刻f2までの推移において，（全状態

塑壁上：シミュレーションの開始時刻をT：＝0とすからの）colaesceとは，（乱）数列人∈［0，1）t2Jtlに対る．空列A＝（）を用意する．反復回数オ：＝1とする．

して，∃y∈n，∀J∈n，◎き…（∬，ス）が成り立っている弛j：シミュレーションの開始時刻をr：＝r−α事と定義する．次に，マルコフ連鎖の推移を定義してとする．

いる更新関数は，COalesceする確率が正（非零）で塑埋j：乱数列∫￣ヱ■＝（オけ］，…，Jけ十α−1］）を生成あると仮定する1．
きじま しゅうじ，まつい ともみ
東京大学大学院情報理工学系研究科〒113−8656文京区本郷7−3−1
2005年5月号

1第一話の節4．1参照．
2 ∃れ，J2∈n，◎言‡（れ，ん）≠◎き≡（Jz，ん）．
3第二話参照．単調性を仮定すると，「全状態からのcoa・1esce」⇔「最大元と最小元からのcoalesce」が成り立つ．
© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

（31）329

−4α−2α−3α−α0

時刻

図2 a遡るCFTPアルゴリズムの概念図

し，数列∫の先頭に挿入する．すなわち，∫‥＝（∫￣〜；の出力は無限の推移の実現値になっている．

∫￣机；…；∫￣1）＝（Jけ］，Jけ＋1］，…，∫ト1］）とす

アルゴリズム1において正整数αの決め方は重要る．

である．アルゴリズム1の計算時間について考えると，拠：nのすべての状態について，共通の乱数列

時間aでcoalesceする確率をDとすると，アルゴリ∫▼fを用いて時刻rから時刻r十αまでマルコフ連ズムの期待計算時間はα也である．すなわち，COa−鎖を推移させる．このとき，1esceするような乱数列が存在しない7ほどaが小さ（a）もし，時刻丁十aでcoalesceしていれば（すないと，アルゴリズムは停止しない．また，αをcoa・わち，］z∈n，∀∬∈n，Z＝◎ぎ＋α（∬，∫￣f）），時刻0の1escence時間の期待値E［71］程度にするとDは定数状態〟＝◎冥∬，∫）を出力し4，停止する．となる8ので，例えば単調CFTPならば期待計算時間（b）そうでなければ，反復回数をオ：＝ダ＋1としては0（E［71］）となり，標準単調CFTPと同程度の計Step2に戻る．

算時間となる．

2．2 前向きcouplingのcoalescence時間
ここで，前向きcouplingのcoalescence時間を議

図2はアルゴリズム1が3反復（g＝3）で終了した論しておく．前向きcouplingのcoalescence時間と様子を表している．
アルゴリズム1は，特に二つの点において，これまは，全状態について時刻0からマルコフ連鎖を推移さでに紹介したCFTPアルゴリズムとは異なるようにせて，初めてcoalesceする時刻T＊のことである9．

見える．一つは「ステップバックする時間がα」であ

時刻r＊での状態は，定常分布に従っているとは限らる点，もう一つは，「第∠反複での終了条件が時刻0ないことに，改めて注意されたい10におけるcoalesceではなく，時刻−ia＋aでのcoa−前向きcouplingのcoalescence時間T＊は確率変1esce」という点である．しかし，アルゴl）ズム1は数である．一方，CFTPアルゴリズムのcoalescence

まぎれもなくCFTPアルゴリズムの一種である．そ時間71も確率変数である11．そして，この二つの確のことを説明しよう．

率変数は同一の分布に従う．もう少し正確に言うと，一つ目の点について，CFTPアルゴ1）ズムでは毎回coalesceを確認する必要はなく，したがって，ス補題1確率変数T＊と71に対して次が成り立つ．

テップバックする時間はどのように決めてもよい5∀t＞0，Pr（71＞t）＝Pr（T＊＞0）．

二つ目の点は少しだけトリッキーである．時刻−わからのシミュレーションを考えた時，時刻−わ＋α

証明：CFTPアルゴリズムで，COalescence時間がtでcoalesceしていれば，時刻0でも当然coalesceし以下となるような数列ス∈［0，1）亡全体の集合をぴぃ］ている．CFTPアルゴT）ズムにおいて，COalesceを確認するためのシミュレーション開始時刻はアルゴリ
ズムの出力に影響を与えない6ので，アルゴリズム17つまり♪＝0．

8来：submultiplicativityから導けますね．松：劣乗法性
ね．

9］z∈n，∀∬∈n，Z＝◎㌃（J，人），かつ∃れ，J2∈
4coalesceしているので，yはxに依存しない．
5第一話の節4．3参照．
6第一話の節4．3参照．
330（32）n，◎㌃￣1（れ，人）≠◎J＊￣1（∬2，ス）．
10第一話の節4．2参照．
11第一話の節3．1参照．

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

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

−α2（1α03α4α

時刻

図3 ReadOnceアルゴリズムの概念図とする．すなわち，任意の数列人∈日月は∃〝∈n，をi：＝i＋1とし，Step2に戻る．
図3はRead Onceアルゴ1）ズム（アルゴリズム2）∀J∈n，〝＝◎㌔（J，人）を満たす．このとき，任意のが3反復（オ＝3）で終了した様子を表している．

人∈日月は］〝∈n，∀J∈n，〝＝◎占（J，ス）も満たす．
また逆も成り立つ．確率変数Aは［0，1）f上の一様分定理2 エルゴード的なマルコフ連鎖に対して，ある布に従うとすると，Pr（T＊≦t）＝Pr（＾∈U［t］）＝Pr（71≦t）．

□

3．Read Onceアルゴリズムαが存在して，アルゴリズム2は確率1で有限停止し，アルゴリズム2の返す値は定常分布に厳密に従う．

3．2 ReadOnceアルゴリズムのアイデア
ずばりRead Onceアルゴ1）ズム（アルゴリズム2）

3．1アルゴリズムと定理

以上の準備をふまえて，Read Onceアルゴリズムのアイデアは，a遡るCFTPアルゴ1）ズム（アルゴリズム1）の模倣である．このことを説明しよう．

を紹介しよう．

Read Onceアルゴリズムがr回の反復で終了した
アルゴリズム2（ReadOnceアルゴリズム）

場合を考える．すなわち，アルゴリズム2が停止した

Step O：正整数α＞0を与える．

時の反授回数オ＝rだったとする．このとき，アルゴ

星座pj：乱数列人1＝り［0］，バ［1］，．‥，ス［α−1］）を生成す

リズム2は〟＊∈nを出力し，アルゴリズム2で〟＊をる12．Qのすべての状態について，共通の乱数列人1決定するのに用いた乱数ベクトルは人＝（Al；人2；…；を用いて，時刻0から時刻αまでマルコフ連鎖を推人r）であったとする（人r＋1は〟＊の決定には無関移させる．このとき，係！）．

（a）もし時刻0と時刻aの間でcoalesceしていれば出力された状態〝＊について考える．アルゴリズム（すなわち，∃z∈n，∀J∈n，Z＝◎岸（J，人l）），変数〝2の中で，時刻0から時刻aの間の推移はcoalesceに時刻αの状態を代入する．すなわち〝：＝◎岸（J，人1）している（すなわち，∃z∈n，∀J∈n，Z＝◎g（∬，ス1））とする．反復回数をオ：＝1とする．

ので，時刻0から時刻771の間の推移も当然coalesce（b）そうでなければ，Steplへ戻る．

していて，〝＊＝那d（J，Å）（∀J∈n）である．

塑埋j：乱数列人…＝（ス［由］，バ［才α＋1］，．．．，ス［（才＋1）α用いられた乱数列人は，a遡るCFTPアルゴリズ−1］）を生成する13．nのすべての状態について，共ム（アルゴリズム1）でも，反復回数がγ回で終了し通の乱数列人＝を用いて，時刻gαから時刻（ダ＋1）αた場合に，そしてその場合だけ現れ，同一の状態〝まで推移させる．特に〟′＝◎〜去＋1）α（〟，人i＋1）とする．

を出力する．このことはRead Onceシミュレーショ

（a）時刻iaと時刻（i＋1）aの間でcoalesceしていれンの時刻を−γαほどずらして見ると分かりやすいば（∃z∈n，∀∬∈n，Z＝◎〜去＋1）d（J，人州）），時刻由の（図3および図2参照）．
まず，∫￣〜＝人￣i＋r十1状態〟を出力し，停止する．
（b）そうでなければ，〝：＝〝′と更新する．反復回数
12便宜上人1と記述するが，この数列を記憶する必要はな（オ∈（1，…，γ））とすると，∫＝（∫一r；∫イ＋1；…；人▲l）＝Åである．このとき，∀g∈（1，…，γ−1）について，∫￣≠を剛−て時刻一才αから時刻rhl＋aまで推移を行ってもcoalesceせず14，i＝し＼

13脚注12と同様，数列人ヱ＋lを記憶する必要はない．
2005年51十l∫rについては∫−rを用いて時刻一−γ8から時刻一用© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

（33）33l＋aまで推移を行うとcoalesceする15．したがって，SteD2：（a）′もし時刻（i＋1）aでcoalesceしていれば∫＝人はアルゴリズム1で，反復回数がγ回で終了し〝を出力する．もしⅣ個のサンプルが出力されていた場合に，そしてその場合だけ出現する．また，そのたら停止する．そうでなければ，〟：＝◎〜㌃1）α（J，入行1）

ときアルゴリズム1が出力する値はy＊＝◎竺r。（J，∫）とし，i：＝1と更新してStep2に戻る．

（∀J∈n）となっている．

最後に数列の出現確率を考える．集合U［γ］をアルこのように置き換えられたRead Onceアルゴリズ

ゴリズム2がγ回の反復で終了するような数列人＝ムでⅣ個のサンプルを生成するのに必要な期待計算（Å1；…；スr）全体の集合とし，集合のγ］をアルゴリ

時間は（N十1）abとなるので，a遡るCFTPアルゴ

ズム1がγ回の反復で終了するような数列∫＝（∫￣r；

リズムとほぼ同じ計算時間で，乱数を記憶することな…；∫￣1）全体の集合としよう．これまでの議論から，

くサンプルが得られることが分かる．

U［γ］＝百レ］である．さらに，アルゴリズム2がγRead Onceアルゴリズムの計算時間はaの決め方ステップで終了する確率はPr（人∈Uレ］）＝♪（1−♪）1￣rに大きく依存する．特にαが小さいと，アルゴリズである．一方，アルゴリズム1がγステップで終了

ムは全く止まらなくなってしまう17．しかし，αをどする確率はPr（∫∈百レ］）＝♪（卜♪）1￣rである16．すなのように決めたらよいものだろうか．次の節ではパラ

わち，ReadOnceアルゴリズムで得られる状態y∈n

メータaを全く与えないでRead Onceアルゴリズムの出現確率はα遡るアルゴリズムでy∈nが得られるを実現する，驚愕のTwin RunアルゴT）ズムを紹介

確率に等しく，したがって，アルゴリズム2はパーフしよう．
3．4 Twin Runアルゴ］）ズムェクトサンプリングを実現する．

TwinRunアルゴリズムのアイデアは，サンプリン

3．3 計算時間
以上の議論で，Read Onceアルゴリズム（アルゴグとは別の乱数列を生成してステップバッタする時間

リズム2）がα遡るアルゴリズム（アルゴリズム1）を決めようというものである．まず準備として，アルを模倣したものであることが分かった．アルゴリズム

ゴリズム中で呼び出す二つの手続き，＊，SuCCeSSfulとからも分かるとおり，Read Onceアルゴ））ズムでは＊−failingを導入する．

乱数を記憶する必要がない．

Read Onceアルゴ））ズムの期待計算時間も，a遡手続きA（＊−SuCCeSSful）

るアルゴリズムと同様であることが分かる．ただ，気塑塑＿＿＿！：シミュレーションの終了時刻rをr：＝1をつけなければならないのは，Read Onceアルゴリに設定する．

ズムのSteplで，COalesceした状態を得るための前塑埋j：二つの乱数んけ一1】，んけ一1］を生成する．

処理が行われていることから，正確には，ReadOnce

便宜上，ん＝（ん［0］，ん［1］，…，ん［r−1］），ん＝（ス2［0］，アルゴリズムで一つのサンプルを得るのに必要な期待ん［1］，…，ん［r−1］）と表す18．数列んとんを用いて，計算時間はα遡るCFTPアルゴリズムの2倍程度でそれぞれ推移を実行する．このとき，あることが分かる．

（a）数列んおよぴんについて，ともに時刻rでさらに，Read Onceの良いところは，アルゴ））ズ

coalesceしている19場合，先にcoalesceしていた乱ムの出力する状態〝と，アルゴリズム2の終了判定

数列の添え字をノ∈（1，2）とする20．同時の場合は1／2のために行った推移のcoalesceした状態（Step2（a）の確率でノを1または2とする．数列んに対する時で得られた状態g）は互いに独立，となる点である．

刻rの状態る＝◎訳∬，ん）を返して終了する．

したがって，サンプルが複数個必要な時に，アルゴリ（b）それ以外の場合，T：＝T＋1としStep2に戻る．

ズム2を連続して実行することができる．具体的にはサンプルがⅣ個必要な場合にはアルゴリズム2のStep2（a）を次のように置き換えれば良い．
14∃Jl，∬2∈n，◎＝ぎ冨十〃（∬．，√り≠◎＝…冨十α（裁，∫￣〜）．

15∃z∈n，∀J∈n，Z＝◎＝㌶＋α（∬，∫￣r）．

16来：一種のPASTAですね（文献［8］参照）．桧：美味しい性質だな．
332（34）17節2．1参月軋
18例えば単調マルコフ連鎖ならば，この乱数列を保存せずにアルゴリズムを実行することができる．
19 ∃z．，為∈n，∀∬∈n，Z．＝◎J（J，ん），為＝◎訳J，ん）．
20例えば．∃g∈n，∀∬∈n，g＝◎J￣1（∬，ん）ならノ＝1．
オペレーションズ・リサーチ© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

手続きB（＊−failing）（A￣∫；人￣…；…；ス￣1）＝（ス［r］，ス［r十1］，…，人［−1］）塑型迅：入力された状態をy∈nとする．とする．

些埋ユ：シミュレーションの終了時刻rをr：＝1に廻：nのすべての状態について，共通の乱数列設定する．

人￣iを用いて時刻Tから時刻r十r￣〜までマルコフ

塑埋＿j：二つの乱数んけ−1］，ん［r−1］を生成する．
便宜上，ん＝（ん［0］，ん［1］，…，ん［T−1］），ん＝（ス2［0］，連鎖を推移させる．このとき，

ん［1］，…，ん［T−1］）と表す．数列んとんを用いて，（すなわち，］z∈n，∀J∈n，Z＝◎享十r￣f（J，A▲り），時それぞれ推移を実行する．特にyl：＝◎訳〟，ん），〝2：＝刻0の状態〟＝◎冥エ，人）を出力し，停止する．

◎訳〟，A2）とする．このとき，（b）そうでなければ，反復回数をダ：＝g＋1として（a）数列んとんについて，いずれか一方でも時刻Step2に戻る．

（a）も し，時刻T＋T￣iでcoalesceしていれば22

Tでcoalesceしている場合，COalesceしていない乱数列の添え字をノ∈‡1，2）とする21．同時の場合は1／2

アルゴリズム3と4が本質的に等しい事の確認は，の確率でノを1または2とする．状態〟J∈nを返し

Read Onceアルゴ1）ズムの場合と同様なので読者にて終了する．

任せる事とし，次ではアルゴリズム4の基本的なアイ（b）それ以外の場合，T：＝T十1としStep2に戻る．

デアを簡単に記す．

図4の上半分は，TwinRunアルゴリズムが4反復以上の準備の下，TwinRunアルゴリズムを記述す（才＝4）で終了した様子を，下半分はアルゴリズム4が4反複（g＝4）で終了した様子を表している．図4のる．

TwinRunCFTP（下半分）の上段は，アルゴリズムアルゴリズム3（TwinRunアルゴリズム）

4で出力する状態を決定する推移を示しており，これ

塑qL主：手続きA（＊−SuCCeSSful）を実行し，出力さは本質的にTwinRun（上半分）の上段と同じものでれた値をJとする．1／2の確率で次の（a），（b）のいずれある．図4最下段の時間軸が右から左へ向かっているかを実行する．のは，値r￣iを決定するのに，（別の乱数列を用い（a）J∈nを出力しアルゴリズムを終了する．て）前向きcouplingを行っている事を表している．
実はアルゴリズム4は，α遡るアルゴリズム（アル（b）Step2に進む．
堅塁旦1：入力x∈flを与えて手続きB（＊一failing）をゴリズム1）では固定されていた遡る時間（α）を，反実行し，変数J∈nを出力された値に更新する．1／2

復ごとにStep2で決定するように変更したものであの確率で次の（a），（b）のいずれかを実行する．

る（ぜひ，アルゴリズム1と4を直接比較していただ（a）J∈nを出力しアルゴリズムを終了する．

きたい）．アルゴリズム4は本質的に2本の乱数列を（b）Step2にもどる．

用いているが，出力を決定するのに用いているのはそのうち1本（上段）だけであり，2本目の乱数列（下Twin Runアルゴ］）ズム（アルゴリズム3）は，次のCFTPアルゴリズムを模倣する．

段）は遡る幅を決定する時計の役割にしか用し、ていなしヽ

さらに，確率変数r￣ヱ■の分布関数は前向きcou−アルゴリズム4（TwinRunCFTPアルゴリズム）

plingのcoalescence時間T＊の分布関数に等しく，

塑埋」：シミュレーションの開始時刻をr：＝0とす

これは，1本目の乱数列（上段）を使ったCFTPのる．空列ス＝（）を用意し，反復回数を7：＝1にする．

coalescence時間71の分布関数にも等しい．したが塑埋j：前向きcouplingを実行し，そのcoales−ってアルゴリズム4の各反復において，（Step4（a）

cence時間をT−t＞0とする．シミュレーションの開で）アルゴリズムが終了する確率はいつも1／2である．

すなわちアルゴリズム4は，COalescence時間71

始時刻をr：＝T一丁￣gとする．
些埋＿＿蔓：乱数列人￣〜＝（人【r］，‥・，人け十r￣i−1］）を生を知ることなく，アルゴリズム1を毎反復1／2の確率成し，数列Åの先頭に挿入する．すなわち，人：＝で終了させる驚くべき仕掛けとなっている．このこと
22正確には，COalescence時間がちょうどT￣Z■の場合は確21例えば，∃Jl，J2∈n，◎J（Jl，ん）≠◎J（J2，A2）ならノ＝2．

2005年5月号

率1／2で（b）に進む．

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

（35）333

Twin Run0rlr2

r3

r4 時刻

r■4T￣3r￣2r￣1 0 時刻

Twin RunCFTP

図4 TwinRunアルゴリズムの概念図から明らかなように，アルゴリズム4の反復回数の期として，Fillはinterruptibleアルゴリズム（Fillの待値はたったの2である．

アルゴリズム）を提案している［1］．Fillのアルゴリ

アルゴリズム3と4が本質的に等しいことから，次ズムもcoalesceの確認を利用してパーフェタトサンプリングを実現する．連続状態のマルコフ連鎖に対すの定理が成り立つ．

るパーフェクトサンプリング法としては，Murdoch

定理3 エルゴード的な有限マルコフ連鎖に対して，andGreenが2段階カップ］）ングのアイデアを用いたTwin RunアルゴT）ズム（アルゴリズム3）が確率1アルゴリズムを提案している［3］．

で有限停止し，アルゴリズムの返す値は定常分布に厳密に従う

今回，パーフェクトサンプリング法について紹介したが，マルコフ連鎖を用いたサンプリングを行う際に．

は，マルコフ連鎖の収束時間が大変重要となる．様々

TwinRunアルゴリズムの驚博すべき点は，パラメータαの設定を不要とする巧妙なアイデアに加えて，なマルコフ連鎖の収束時間の算定に関する研究が，現在も行われている．

手続きA，Bを持ち込むことで，記憶領域を殆ど必参考文献要としないという成果を同時に達成している点であろ［1］］．Fill： Aninterruptiblealgorithmforperfectsam−

う23pling via Markov chains，

4．まとめ

Probability，8（1998），131−162．

3回の連載で，CFTPに基づくパーフェクトサンプリングについて述べてきた．CFTPアルゴリズムは非常に面白いアイデアであるが，「パーフェクト」という性質は非常に繊細で，安易な改訂は危険である24．
マルコフ連鎖を用いたパーフェクトサンプリング法
23釆：どうやってこういう事考えつくんでしょう？ 二松：
確率論とアルゴリズムの両方のセンスが少なくとも必要だ

ねぇ．来：発想も天才的ですよね．桧：初めて読んだ時ホント感動したよ．
24釆：以前，僕も改良を考えたんですけど，よく考えたら間違えてました．松：発表する前に気づいてよかったね．
334（36）The Annals of Applied［2］0．Haggstr6m： Finite Markov chains and algor−
ithmic application， London MathmaticalSociety，Student Texts52，CambridgeUnivesity Press，2002［3］D．］．Murdoch and P．J．Green： Exact sampling
fromacontinuousstatespace，

ScandinavianJournalofStatistics，25（1998），483−502［4］D．Wilson：

Howtocouplefromthepastusingaread−OnCeSOurCeOfrandomness，RandomStructures

Algorithms，16（2000），85−113© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

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