川……f…l…‖‖＝＝‖‖‖＝‖＝‖‖＝＝‖‖＝＝＝‖＝＝＝＝＝‖‖‖＝＝‖‖＝‖‖‖＝＝仙＝＝ll…＝‖＝‖＝＝＝‖‖＝‖‖‖＝‖‖＝‖‖‖＝＝‖＝＝＝‖＝＝‖＝‖＝‖＝‖‖＝＝‖‖＝＝‖＝‖‖＝‖＝‖＝＝＝＝‖‖＝＝＝‖＝＝‖＝＝‖‖‖‖‖＝＝‖‖＝‖‖‖‖‖‖＝‖＝＝＝＝‖＝＝‖＝‖＝＝削l

複数の藤生虚刹を措肇確率過程の■∵ 二∵∵㍉．二、一、・
MarvinK．Nakayama（訳：か沢 利久）
1…ll‖‖‖＝‖‖＝‖＝‖＝‖＝‖＝‖‖‖‖＝‖‖‖‖＝＝‖＝‖‖＝＝‖‖‖‖‖‖‖＝‖＝‖＝‖‖＝＝‖‖＝＝‖‖‖………L………ll＝＝‖‖＝‖＝＝‖‖＝‖‖＝‖州‖＝川＝lL…＝‖‖＝‖＝‖‖＝＝＝＝‖‖＝＝川…l………＝‖‖＝‖＝川＝‖‖‖＝‖＝‖＝‖‖‖‖＝＝＝‖‖‖川‖

SAVE（theSystemAVailabiiityEstimator）はその−−． ‥．．、∴■・「・∴成功事例のひとつといえる。SAVEは，複雑なコン実システムを表現した数学モデルの多くは大規模か

ピュータシステム，特にフォールトトレラントンステ

つ複雑であり，解析的な手法によって評価することはムの信頼性を様々な角度から評佃するためのシミュレ

一般に難しい￠ そのような場合によく用いられるのがーションかシミュレーションである。シミュレーション結果を正

トトレラントシステムとは，非常に高い信頼性が要求しく解釈するためには，出力データを統計的な手法にされるシステムのことであり，航空機の管制，銀行のよって分析する必要があり，伝来酎大間の導出はそのひ

業務処理9 運送会社における荷物の追跡などに使われとつとなっている。

ている∴嵐軋 これらのシステムでは冗長構成を取る

ソフトウェアパッケージである。フォール

精度が漸近的に保証される信頼区間を得る方法としことによって高伝来酢巨を実現している。例えば，ディて9 両生過程（regemerative process）を対象としたスク損傷によるデータの損失を防ぐために，複数のデ

府年過程法（regenerative method）がある

ィスクに岡じデータを保存しておいたり，プロセッサ［11〜13，15，16］。再句ミ過程とは，プロセスが確率的にの故障を考慮していくつかの予備プロセッサを備えて再開始するような時刻（再生点）を無限個含む確率過いたりする。フォールトトレラントシステムは非常に程のことであり，多くのモデルがこの性質を満たす。

大規模かつ複雑であるため，解析的な方法による分析

例えば，安定な単一一一サーバ待ち行列モデルでは，シスは困難であり，シミュレーションによる評価が主とな

テムが空の時に到着した客の到着時刻が再生点に対応っている∴一般に′義子部品の故障時間は指数分布に従する。また9 既約な有限マルコフ連鎖では，特定の状うとされている。よって，修理時間も指数分布に従う

態への訪問時刻がそれにあたる8 二再生過程の見本路と仮定すると，対象としているシステムのモデルはマ（sample path）は番生魚の列（再生点列）によってルコフ連鎖となり，SAVEパッケージを用いて再生独立で同一の分布（i。iud）に従う複数の部分（再生過程法により評価尺度の信頼lズ間を求めることができサイクル）に分割される。マルコフ連鎖の場令は，ある（SAVEパッケージの詳細についる状態を訪れた時刻から再びその状態を訪れるまでのび［21〜23］を参照）。

ては文献［1，2］及間がひとつの府生サイクルとなる。府生過程法では，ところで，既約な有限マルコフ連鎖では，どの状態再生過程が持つこの性質を応用して性能計佃りく度の信に対しても，その状態への訪問時刻の列は再生点列と頼区間を求める。再隼過程法については2章でより詳なるため9 複数の再生蔦列が存在する曲 他の多くの再しく触れる。また，文献［27］も参考となる。

生過程もこの性質を持っている。しかし9 再生過程法

再生過程法は，生産や在棒，輸送，通信，金融といでは，ひとつの内生点列しか利用しておらず，このよった幅広い分野に応用可能であ県］旧Mトーマス。

うな付加的な構造を活かしていない。本論文の主題は，ワトソン研究センター（ニューヨーク）で開発された

このように複数の再生点列を持つ確率過程のシミュレーションについて論じていくことにある。得られる推定量の多くは，標準的な 方法によるものよりも小さな

Marvin K．Nakayama

DepartmentofComputerandInformationScience
NewJerseyInstitute ofTechnology
Newark，NJO7102i982，U．S，A．
（翻訳）おざわ としひさ駒澤大学 経営学部

閑鳩（10）

分散を持つことが確かめられており，その意味で，分散減少法（variance reduction method）のひとつと見ることもできる（文献［24］の11章を参照）。さらに，
これら手法は計算掻の増加を抑え，効率的に実装する© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

オペレーションズ￠リサーチステムの評佃では，芽をシステム状態の時間的推移ことも可能である．

以下では，複数の再生点列を持つ確率過程をシミュ
レーションするためのふたつのアプローチについて説明していく．ひとつめは，3章で述べる置換再生過程

推定量（purmutedregenerativeestimator）である．
もうひとつは4葦で述べる準再生過程法（semi−
regenerative method）である．次章では予備知識として，標準的な再生過程法について説明する．これら説明の中では，議論を簡単にするため，対象となる確率過程を有限な状態空間5＝（0，1，2，・・・，∫）上の既約な

離散時間マルコフ連鎖X＝〈先：ノ＝1，2，…）に限定する．もちろん，取りあげる手法は複数の再生点列を含むより一般的な再生過程にも通用可能である．

とし，βをシステムが動作可能である状態の集合にとる．

再生過程法では次の式を基にγを推定する．
且〃［∑た狂）￣1′（＆）］
プ′＝

（1）

居び［れ（1）］

ここで，且びは確率過程が状態〃から開始された場合の期待値を表す［12］．∽個のれサイクルから再生過程法を使ってγを推定するためには，まず，y（々）＝∑た某社1）′（先）と置き，ア（∽）＝（1／∽）∑欝＝1y（点）によって且〃［∑たどト1′（先）］の値を推定する．次に，r（々）＝r（々）−れ（々−1）と置き，F（∽）＝（1／∽）∑芳＝1
r（々）によって且〃［r（1）］の値を推定する．最後に，∽個のれサイクルを基にした，γの再生過程推定量を2．再生過程法

次で与える．

マルコ7連鎖＿方がある固定された状態〃∈Sから開始したとする（品＝〃）．時刻（停止時刻）の列れ戸（∽）＝

推定量ダ（椚）は次の中心極限定理を満たす．

＝（れ（0），れ（1），れ（2），…）を，マルコフ連鎖ガが状態〃∈5にあった時刻の列とする．例えば，n（0）＝0∽1／2戸（∽）−γヱⅣ（0，1）♂rは初めて状態〃にあった時刻，れ（1）は次に状態〃にこ こ で，♂要＝（且［y（1）2］＋γ2月、［丁（1）2］

あった時刻である．uを復帰状態（returnState）と−2プ宜［y（1）丁（1）］）／（且［㌻（1）］）2であり，ヱは分布の意

呼び，時刻1（点−1）とれ（点）の間におけるマルコフ

味での収束，Ⅳ（0，1）は平均0，分散1の正規型確率

連鎖ガの挙動を々番目の（再生）れサイクルと呼ぶ．

変数を表す［27］．γの信頼区間は中心極限定理より以

今，マルコフ連鎖Xを時刻0からn（∽）までシュミ下の手順で得られる．まず，E［y（1）2］の推定量をレー卜し，∽個の再生れサイクルが得られたとする．

Ⅵ（∽）＝（1／∽）∑袈1y（々）2で，且［㌃（1）2］の推定量を

マルコフ性より，各々＝0，1，2，…に対して，確率過程

佑（∽）＝（1／∽）∑悸＝1r（々）2で与える．次に，（方r両）＋ノ：ノ＝0，1，2…）は同一の分布に従う（確率過且［y（1）㌃（1）］の推定量を 佑（∽）＝（1／∽）∑㍑＝1

程芽が時刻n（0），r（1），れ（2），＝・で確率的に再開始

y（々）㌻（点）で与える．最後に，♂γの推定量を∂γ（∽）2する）ことが分かる．さらに，各れサイクルも独立＝（析（∽）＋戸（∽）2tち（研）−2ダ（∽）惰（∽））／テ（∽）2で与で同一の分布に従うことが分かる．

える．これらを用い，100（1−∂）％信頼区間（0＜∂＜′を状態空間5上の報酬関数とし，マルコフ連鎖1）は，十分大きなサイクル数∽に対して方が状態Jを訪れた時の報酬を／（∬）とする．ここで［ダ（材トz（∂）∂γ（∽）ニ′…＼．Z（∂）∂r（研）は，定常状態における平均報酬γ＝1im卜∞（1／才）∑三三占／（先）について考える．例えば，方を安定な離散時
間待ち行列モデルの待ち行列長過程とし，報酬関数／を／（∬）＝Jで定義すると，γは定常状態における待

∽1／2

ダ（∽）＋

］（2）

∽1／2

で与えられる．ここで，Z（∂）はP（Ⅳ（0，1）≦z（∂））＝1−∂／2を満たす値である．例えば，95％信頼区間では，
z（0．05）＝1．96となる．

ち行列長の期待値となる．また，ある状態部分空間

再生過程法は，定常状態における平均報酬以外の量β⊂5に対して′（ズ）＝1（∬∈β）ととれば，γは確率を推定するのにも利用できる．例として，ある部分空過程方が状態空間β内にある定常状態確率となる

間F⊂Sに到達するまでの時間の期待値cを考える．

（1卜）は，条件（・）が真であれば1，そうでなければ高信頼システムでは，システムが故障となる状態の集0の値をとる関数である）．後者の例は，フォールト合をFにとることで，システムが故障するまでの平

トレラントシステムの定常状態における稼働率を求め

均時間cを評価することができる．確率過程方が初る場合にも用いられる．稼働率とはシステムが動作可

めて状態空間Fに到達した時刻を了1とすると，一能である時間の割合を指す．フォールトトレラントシ

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

（11）175

C＝に増加することを示した。よって，長い見本路に対し且む［min（苅（1），7宣）］見直（7÷く弟（1））］て重点標本抽出法を適用すると，推定量の分散が非常で与えられる［23］n cに関するこの式を，再生過程法
では次のようにして利用する。まず，点番口の再生
苅サイクル（彪＝1，2，…，∽）に対し，そのサイクル開始後，マルコフ連鎖が状態空間ダに初めて人った時に大きくなる可能性がある。しかし，再生過程法を川いることで，見本路を複数の再生サイクルに分解して分析することができ，結果として分散を小さくすることが吋能となる。

刻をrF（点）＝min（ノ＞苅（点−1）：Ⅹメ∈椚とする。
g〃［min（弟（1）∴持）］の推定量をⅣ（椚）＝（1／椚）∑㍑＝1

min（右（ゑ）∴持（ゑ））で与え，g〃［1（7宣＜苅（1））］の推定量をβ（刑）＝（1／椚）∑漂＝11‡二号（々）＜苅（点））で与一える。
これらを川いて，Cの再生過程推定量は

一方∴再生過程法にも幾つかの欠点がある。まず，

この方法でも，初期偏差とは異なる種類の偏り（bias）が生じる。例えば，ふたつの標本平均の比の期待値とそれら標本−、ド甚」の期待値の比は基本的に異なるため，庸雄過程法で求めた其朋寺値の比の推定量も一般には偏りを持つ。また，モデルによっては再生サイ∂（研）＝器で与えられる。∂（椚）を基にしたcの信頼区∬附こついても，式（2）で示したγの信束酎ズ間を求めるのとi†てJ様

クルが非常に良くなるという問題が牛じる。例えば，
多くのノードを持つ持ち行列ネットワークのような大規模システムをシミュレートすると，再4一たの間隔はな手順に従って得ることができる凸再生過程法が適用できる他の興味ある評価量として，
確率過程Xの時間平均分散係数62（time−aVeragevariance constant）［27］や無限期間別引コストの期待値［17］などがある9 時間平均分散係数♂2とは，確率過程斉の時間平均度f＝（1／f）∑j三占先に対して中心極

限定理を適応する時に現れる分散係数であり，克を基にしたγの信頼区間を求めるのに必要となる［27］。

非常に長くなるであろう。これは，得られる再生サイ
クル数が少なくなることをホしており，推定量の精度保証が難しくなることを意味している（詳細については又献［封のp。99を参照）。以下で説明する方法は，
拍隼過程法が持つこの欠点を改善したものとなっている。

3。置換再生過程推定量

割引コストの方はファイナンスなどの分野に現れるl‡りCalvinとNakayamaは，ふたつの再生J〔㌔（列を含む題である。

再生過程法は幾つかの興味ある特徴を持っている8

まず，塊平均法（methodofbatchmeans，文献［24］のpp．528529を参照）など，定常状態におけるrlそ均値を推定する幾つかの方法が持つ初期偏差（initial−
ization bias）の問題を含まない点が挙げられる。初期偏差とは，定常状態をあまり代表しないような状態からシミュレーションを開始した場合に生じるものであり，推定量の精度に影響を与える。例えば，非常に混雑している待ち行列ネットワークのシミュレーショ

ンを，系内に客がいない状態から開始したとすると，開始してしばらくの間は，定常状態にあるシステムでは稀にしか起きないような状況をシミュレートすることになってしまう。その他の特徴としては，時間平均分散係数の推定量の中で最も速い収束率を持つことが挙げられる8 また，重点標本抽出法（importancesampling）を用いる場合においても大変役に卓二つ（重点標本抽揖法については文献［3］の第2章を参照）。
Glynn［20］は9 重点標本抽出推定量の構成要素である
見本路尤度比の分散が9 見本路の長さに従って指数的確率過程を対象とした置換再生過程推定量を文献［6］で初めて提案し，続いて文献［7］でその方法を任意の数の再隼点列を含む確率過程へと拡張した。今までと同様に，時点列苅＝（1（0），右（1）∴石（2），…）をマルコ

フ連館仁方が状態ル∈5を訪れた時刻の列とし，時点
列右＝〈茄（0），7ち（1），7ち（2），…）を他の状態弘∈5を訪れた時刻の列とする。この確率過程艮∵tで定められた性能評価尺度αの推定を考える。αの例としては，方が状態空l悶Fに到達するまでの時間の期待値cや時間軒均分散係数げ2が挙げられる。

ふたつの両性点列を含む場合を例に，置換二再生過程
推定量の基本となる考え方を説明する，図1はそのための見本路である申 簡単化のため，連続な状態空間5上の連続な見本路とした。プ「点列は状態〃への訪問
峠亥Ijの列，7ち一−ごさく列は状態紺への訪問時刻の列である。

確率過程方のシミュレーションを行い，匝†1の卜方の見本路を得たとする。て打点列については，例えば，

時刻れ（0）から右（1）までが最初の再4サイクル，時
刻右（1）から了1（2）までが次の再生サイクルとなる。
この見本路には刀て＝5仰の再甘∵右サイクルがあり，耶帽（12）© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

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

計算することなしに∂＊の値が得られる陽な表現式を与えており，計算量の増加も予想されるほど大きくな
0riginal
SamplePathい．この表現式では，シミュレーションの実行中に幾つかの付加的な値を保持しておくことで，記憶容量や7て（（l）

7て（＝1（2）γ（3） T（4） 右（5） r

計算量の増加を最小にできる．文献［6］での数値実験
では，計算時間の増加は最大でも5％であったが，分∫散は8分の1にまで減少した．置換再生過程法は非常ll−に汎用的であり，他の分散減少法と細▲み合わせて利用

A Permuted
SamplePath

することも可能である．また，文献［8］では，幾つかの置換再生過程推定量に対する中心極限定理を示してV

1て0） γて1） γて2）

雪て3）祁4） 不て5）J

図1ある見本路とその置換見本路

いる。これにより，置換再生過程推定量を基にした信頼区間を求めることが可能となる．

ふたつの再生点列を含む場合について置換再生過程
それを碁にして標準的な再生過程推定量∂を求める

推定量の基本的考え方を説明してきたが，この考え方ことができる．右点列8；は肪＝4個の再生点があり，はもっと多くの再生点列を含む場合へと拡張できる

肱−1＝3個の再生サイクル（右サイクルと呼ぶ）が［7］．また，文献［4］では，特定の置換再生過程推定量存在する．ここで，烏番目の右サイクルとは時刻がある種の意味で最適なことを示した．特に，離散時

右（々−1）から右（点）までの間を指すものとする．図

間有限マルコフ連鎖の場合は，全ての再生点列を用い中では，各右サイクルをそれぞれ異なる線種で描いて構成した置換再生過程推定量が分散最小不偏推定量

た．

となる．
次に，7ちサイクルを置換することによって元の見置換というアイデアがどの程度効果的なのかを示す

本路から新しい見本路を構成する．このような見本路

ために，ここでは，離散時間版アーラン損失システム

は合計で（脆−1）！＝6個作れる．図1のふたつ目のグ（Erlangloss system）のシミュレーション結果を示

ラフはそのような置換見本路のひとつを示しており，

す．こ？モデルには，S個のサーバを持つひとつのス

元の見本路における3番目の右サイクルが1番目に，テーションがあり，そこへ客が一定の率で到着する．

1番目の右サイクルが2番目に，2番目の右サイク

待ち合い室はない．∫個のサーバが全て稼働中の時にルが3番目にきている．結果として得られる再生点の到着した客はサービスを受けることなしに即座に立ち新たな列をそれぞれ7T，了首とした．置換見本路にお去る．このモデルは，コールセンターの分析に用いらけるれサイクルの個数は元の見本路と同じになるが，れるが，そこでは，電話を掛けてくる人が客に，電話れサイクル内での芽の挙動は異なってくる．そこで，回線がサーバに対応する．

（脆−1）！個の置換見本路各々についてαの推定量を計

このモデルは状態空間5＝（0，1，2，…，ざ）上の離散時

算し，それら（脆−1）！個の推定量の平均を新たな推間マルコフ連鎖X＝（X＝ノニノ＝1，2，…）となる．そこ

定量∂＊とすることを考える．文献［6］では，妥を元で，その推移確率行列を，非ゼロ要素をP（Z，オ＋1）＝の見本路，妥（1），妥（2），…，妥（Ⅳ），Ⅳ＝（脆−1）憧置換見s／（2才＋s）＝1−P（Z，才11），0＜ブ＜s，P（0，1）＝P（s，S

本路とした時，各置換見本路ガ（烏）が元の見本路妥と−1）＝1としてP＝（P（J，〝）：∬，〝∈5）で与える．再同じ分布に従うことが示されている．そこで，置換再生点列をある状態への訪問時刻の列に対応させ，ひと生過程推定量をつの再生点列のみを基にして構成した推定量（再生過程推定量）と全ての系列を同時に用いて構成した置換∂＊＝去鼻∂（妥（ゐ））再生過程推定量の両方を求める．実験では，れサイで与えることにする．∂（妥（烏））は々番目の置換見本路
妥（た）から求めた推定量である．

状態〃をいろいろと変えてシミュレーションを行っ

文献［6］は，置換再生過程推定量∂＊が標準的な推定量∂と同じ平均値を持ち，分散はより小さくなることを証明した．さらに，全ての置換見本路を厳密に
2001年4月号

クルを復帰状態〃を訪れた時刻の列とし，その復帰た．

表1，2は，サーバ数が∫＝50とぶ＝100の場合について異なる3組の実験データを用い，時間平均分散係© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

（13）17丁表1アーラン損失モデルにおける時間平均分散の推定（サーバー 数s＝50）

いられている他の方法に共通する。例えば，boot−strapandjackknife［14］といった統計分野における再

標本分散
即

25

比

標準 置換
の。の028

言える．この点はシミュレーションや統計の分野で用摘出法（resamplingmethod）は，ある特別な方法によりデータを再利用して推定量を得ている。関連する

0。0022

1．3

その他の〟法としては，U統計法（U−Statistics9 文

豆0 り。0045 ￠。0022 2，1

献［26］の第5章）やⅤ統計法（V−Statistics，文献㌔」 簿どl畠l 二l∈我冒l．貞≠弥‡ I［25］），置換テスト法（permutation tests，文献［ユ0］など）がある．．制御変量法（controlvariates）と呼表2 アーラン損失モデルにおける時間平均分散の推定

ばれる分散減少法もそのひとつである（文献［24］の11．4節を参照）。

（サーバー数s＝100）

4申 準閏生過程法

標本分散
比

即 標準 置換

複数の内生点列を活屈するもうひとつの方法として，
Calvin，Glynn，Nakayamaによって開発された準4（〕 0。0且且5 0。0の39 2−9頂生過程法（semiーregenerative method）［5］がある。
この子法は9 準再焦過程（semi−regenerative proc−ess）［9］と深く関係しており，名前はそれに由来する。
今，ある確率過程上で定義された時点列（停止時刻の数♂2を推定した結果である。各実験では，独立なシ列）ア＝〈了1（0），r（1），r（2），…）を考え，時刻T（オ）とミュレーションを1，000固実行した。これら1，000回

r（々）での状態が同じ（斉r（り＝ガr（ゐ））であったならば，のシミュレーションから得られたデータを基に，時間

時刻ア（オ）以降における挙動（斉r（腑：ノ＝0，1，2，■・・）

平均分散係数げ2の標準推定量（再生過程推定量）のと時亥りγ（点）以降での挙動（斉r（射＋ノ：ノニ0，1，2，…）が標本分散を求め，それを各表の第2列臼に示した。各同じ分布に従うとする。このような時点列を持つ確率

表の第3列冒は同様にして求めた置換再生過程推定量

過程が準再生過程である。例としては，有限な状態空の標本分散である。最後の列にはこれらふたつの標本

間S Lの既約なマルコフ連鎖ガ＝（先：ノ＝0，1，2，…）分散の比を示したQ 結果が比較可能となるよう，復帰が挙げられる。この例では，ある状態集合』＝（∬l，状態がの値毎に9 個々のシミュレーションで平均的

一方2，…，．r。）⊂Sに対し，芽が集合Aを訪れた時点列を

200，000の状態推移が起きるようサイクル数を調整しTにとればよい。集合Aがひとつの要素からなる場

た。

合，時点列㌻は再生点列となる。以下では，このマそれぞれの表において，最初の行は〃の値が確率

ルコフ連鎖Xを川いて準再生過程法を説明すること

過程ガの定常状態における平均値に等しい場合であにし，呼止列アを集合月＝（Jl，∬2，…，∬。）への訪問時

る匂 他の行では〃の値をその平均から徐々にずらし

亥Ijの列と定義しておく㊥ていった。表より，標準推定量の分散は〃の値が平再生過程法では，ある状態〃∈∫への訪問時別の列

均値から離れるほど大きくなっているのに対し，置換を薄生蔦列とし，それを基に統計データを集めた。こ再生過程推定量ではほぼ同じ値になっていることが分れに対し，準再隼過程法では，集合』への訪問時刻

かる。これは，〃の値が平均値からずれるほど置換にによって分割された見本路を基にデータを集める。以よる分散減少効果が大きいことを示している（数値例

卜では，集合』への連続するふたつの訪問時刻で区では最大で66分のユ となっている）曲 さらに，モテリレ

切られた部分（見本路の−一一部分）を軌道（trajec−のサイズが大きくなるほどその効果も大きくなっていtory）と呼ぶことにする。マルコフ連鎖がある状態

るように見える∂ヱ7■∈』から開始されたとし，r（0）＝0，ア（1），T（2），…

置換二麻生過程法は，置換という操作によってシミュをマルコフ連鎖Xが集合月内にあった時刻の列とす

レーションからより多くの情報を引き出し，データをると，時刻ア（点】1）とT（点）の間におけるマルコフ

再利用することでより良い推定量を得る方法であると連鎖の挙動が々杏仁二】の軌道となる血 確率過程Ⅳ＝

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

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

（航：々＝0，1，2，‥・）をⅥち＝方r（妨 すなわち，l鶴がめており，また，それらの式で示された値を推定する

時刻T（点）でのマルコフ連鎖ズの状態となるように方法についても述べている．さらに，この文献では，定義する．このⅣがマルコフ連鎖となることは直ぐ

状態数7のマルコフ連鎖について，γの標準的な推定に分かる．β＝（斤（∬f，∬ノ）：Z，ノ＝1，2，…，d）をそのⅣ

量と階層化推定量を比較した結果も示している．ただの推移確率行列とする．尺（ェz，JJ）は状態Jz∈Aからし，階層化推定量は最適な重み付けで階層化抽出の方開始された軌道が状態み∈Aで終了する確率である．

法を通用して求めたもの，標準的な推定量は元のマル
コフ連鎖芽のシミュレーション結果から時間平均

まず初めに，第2章で取りあげた定常状態における期待報酬γについて考える．準再生過程法では次の（1／乃）∑告‡先を計算して求めたものである．比較の式を基にしてγを推定する［5］．

結果，層別法を用いることで，分散がほぼ12％減少

γ＝するのが分かった．

∑湖（3）∑ふ∈Aβ（∬ォ）g∬z［r］

上記では層別法を用いて平均を推定する場合を例に

この式において，rは確率過程Xが初めて集合Aに戻るまでの時間，且J∫は初期状態をJgとした時の期待値である．p＝（β（∬ど）：Jォ∈A）はⅣの定常分布
であり，方程式β＝β斤と条件p≧0，∑∬ォ∈。β（∬ォ）＝1を満たすべクトルである．

の場合へも適用可能である［5］．また，準再生過程法
は，時間平均分散係数や平均到達時間，モデルパラメータに関する性能評価尺度の微分（尤度比法を利用）など，さまざまな量の推定にも応用できる［19，28，29］．

式（3）を基にしたγの推定量は，階層化抽出の方法（stratifiedsamplingapproach）を用いて構成することができる．この方法では，i番目の層（stratum）
が状態ごど∈Aから開始された軌道に対応するようにd個の層を定義する．れ勿，…，如を和が1となるよ

うなd個の正数とし，乃を総抽出見積り（total
sampling budget

準再生過程法を説明したが，準再生過程法はそれ以外）とする．このnは生成させる軌

道の総個数となる．まず，Z＝1，2，・・・，dに対して，状態∬fから開始された軌道をb沼」個生成させる（実数∂に対し，L則は∂を超えない最大の整数を表す）．状態∬古から開始された軌道は初めて集合Aへ戻った時
点で終了となるが，且J〜［∑長け（先）］と且∬z［r］の推定量をこの軌道上での標本平均として与える．Ⅳの定常分布βは次のように推定する．まず初めに，lγの推移確率行列βの要素斤（ェz，み）を，状態∬古から開始された軌道が状態ズノで終了する割合として求める．
大きさ乃の抽出見積りから得られた斤の推定量を斤乃としておく．次に，定常方程式β〃＝p〃月乃をp乃について解き，それをpの推定量とする．以上の結果得られた，P，且Jど［∑㍍1／（先）］，且∬7［r］の推定量を式（3）に代入して，γの推定量戸刀が得られる．文献［5］では，戸和が中心極限定理を満たすことが証明されており，その結果を用いて戸乃を基にしたγの信頼区間を求めている．

微分情事鋸ま確率的最適化アルゴリズムにおいて有益なものであり，また，システムのクリティカル部分を特定するのにも役立つ．例えば，高信頼システムの設計では，各装置の故障率に関するシステム平均故障時問の微分を計算し，その微分値が最も大きくなるような装置群を改善対象とすればよい．さらに，文献［5］では，準再生過程法と重点標本抽出法を結び付ける方法
や，それを用いて偏りの少ない推定量を求める方法についても述べている．

準再生過程法を重点標本抽出法と一緒に用いることの利点は，重点標本抽出法の実装を，扱っている問題に合った形へと作り直せる点にある［5］．また，準再生過程法と重点標本抽出法を一緒に用いて生成した軌
道は，再生過程法を一緒に用いて得られる再生サイク
ルよりも一般に短い（すなわち，より少ない状態推移しか持たない）．Glynn［20］が示しているように，尤
度比の分散は状態推移数に従って指数的に増加することから，軌道が短いというこの事実は，準再生過程推定量がより雑音の少ない推定量であることを示唆している．

尤度比微分法（1ikelihood ratio derivative
method）の場合も，軌道からデータを集めた方が有利である．文献［18，28］での分析結果は，尤度比微分法の分散が観測時間の長さ（状態推移の数）に対してこのような層別法（stratincation scheme）には，標準的な再生過程法にはない自由度（層の重みれ勿，…，如の設定）が付け加わる．文献［5］では，モデル線形に増加することを示唆しており，再生過程法を用いるよりも準再生過程法を用いた方が分散がより小さくなると考えられる．

に依存する形ではあるが最適な重み付けを表す式を求

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

（15）179

［5］Caivin，J．M．，P．W．Glynn and M．K．Nakayama：∴ −！・一†、・■．．．

The semトregenerative method of simulation output
analysis．Fbrthcoming（2000）曾

複数の二再生点列を持つ確率過程をシミュレートする［6］Calvin，J．M．and M．K．Nakayama：Usingpermuq

ふたっの方法として，置換再生過程推定量と準再生過tationsinregenerativesimuiationstoreducevariance．

程法について述べてきた。これらの手法は，標準的な∴1（1．1ノハi．・〃ヽ．．り／．りり ノ．ノ／／／．ご√川 ／（、（り〃ノり／／（・Jヾ／−〃‖／／（〟∴り／．バ

再生過程法が無視してきた付加的な構造，複数の内生（1998）153−193．

点列を持つという構造を活周したものであり，結米と［7］Calvim，］．M．andM．K．Nakayama：Simuiationof

してより′』、さな分散を持った推定量が得られる。また，

processeswithmultipleregenerationsequence．P7VbaN

これら手法は非常に少ない計算オーバーヘッドで実行／り■′′｛∴∴∫ノしイ．、ぺ【′−／／．．、／、／■吋（晶ト／′亘りソル′・′′1り／．ヾ．、∴・ノ∫一日∴1・！
できる効率的な方法でもあるゎ これら手法のほとんど（2000）179−2椚．

は，対象とする確率過程が複数の二再生点列を持ってい［8］calvin，J．M。and M．K．Nakayama：Centrallimit

れば適用可能であり，モデルが持つ固有の特性には依theoremsforpermutedregenerativeestimator．（砂e7m−

存しない。その意味で，非常に汎川的なものであり，∠go裾5忍βSeαプでゐ，48（2000）776仙787．

他の分散減少法と一緒に用いることも可能である。

；−・∴：i・、巨．／．‥ ．・．・・＼・．∴ ∴ ／−ト・．．．
（Prentice、Hall，1975）．

本論で示した手法は，標準的な二再生過程法が持つ大［10］conover，W．J∴P用CfブcαJAb坤αm∽efγ才c5■由姑オぬきな欠ノ盤のひとつ，非常に長い滴償サイクルが生成さ

Sgco犯dgd才fわ裾．（Wiley，1980）．

れる可能性を解決するものである。この欠点を持った［‖］crane，M．A．andDonald L．Iglehart二Simulating

モデルでは，非常に少ない数の再僅サイクルしか得ら

Stable stochastic systems，Ⅰ：Generalmultiserver

れない可能性がある爛 再生過程法を碁にした推定量やqueues．′〆∠ゐg月Cガ，21（1）（1974）103ml13．

信根I未聞は，二醇生サイクル数が大きくなるに従って漸［12］Crane，M．A．andI〕onald L．Iglehart：Simulating

近的にその精度が向」上していくものであり，数脂lの再Stablestochasticsystems，II：Markovchain．Jqf the

生サイクルしか含まないようなデータから得られた結ACノ好，21（1）（1974）114123．

果にはあまり良い精度は期待できない。それに対し，

L13］cTane，M。A．and Donald L．Iglehart：Simulating

準二再生過程法では，両生サイクルよりも短い軌道を基

Stablestochasticsystems，III：Regenerativeprocessesにデータを集めるので，より多くの観測値を利用でき

and discreteeVent Simuiations．（フわe7dions Research，る。

23（i975）3345．

［14］Efron，B∴Bootstrap methods：Anotherlook at
参考文献

thejackknife．』花形αゐ〆5出め〜才cs，7（1979）126．

［15］Fishman，G．S∴Statisticalanalysis for queueing［1］Blum，A P．Heidelberger，S．S．Lavenberg，M．K．
simulations．Manqgement Science，20（1973）363T369．
NakayamaandP．Shahabuddin：SystemAvailability
Estimator（SAVE）language reference and usersmanualversion4。0．厨gseαⅣゐ屈ゆ0γオ屈』2ブタ∫ノβノば［16］Fishman，G．S∴Estimationin multiserver queue−
ingsimulations．（努erations Research，22（1974）72−7臥［17］Fox，B．L．andP．W．Glynn：Estimatingdiscountedr∫ 肋ね0死点gseαγCゐCe邦′eγ（1993）．

costs．肋視dge椚β乃f5cオg裾Cβ，35（1989）1297−1325．
［2］Blum，A．，P．Heidelberger，S．S．Lavenberg，M．K．
［181Giynn，f），W．：Likeiihood ratio der・ivative
Nakayama andP．Shahabuddin：Modeling and anal−

ysis of system availability using SAVE 77te237d

estimator：An overview．771e1987 mnter SimulaH／山上ノ〃（イ／∴り／（イ∴†ll招／−t心川．− り〃／▼1州／′／71りレ′イ′リ／（1日ノ／／）／．イ∠オ〃乃Co現存柁邦乙、g（1987）366−375．

［19］Glynn，‡）．W∴Likelihood ratio derivative

オプ曙（1994）137−141．

［3］Bratley，P．，B．L。Fox and L＿E．Schrage：A
to Simulation，Second Edition，（SpringerVerlag，

estimators for stochastic systems．Comm．ACM，33（1990）75−84．

［20］Glynn，P．W∴Importance sampling for Markov

1987）．

［4］Calvim，」しM．，P．W。Glynn and M．K．Nakayama： chains：Asymptotics for the variance．Stochastic〃odeね，10（1994）70卜717．
Onthesmall−SamPleoptimalityofmultipleregenera−tion estimators．me1999 mnter Simukltion Co犯々r−
純化∫‰柳瘍ノ旺玩肌（1999）655−661．

瑠＄の（16）［21］Goyal，A．，W．C。Carter，E．de Souza，S．S．Laven−berg and K．S．Trivedi：The system availability© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

オペレーションズ。リサーチestimator．77ze Six：teenth＄ymposium on Fbult−7bler−

jackkni丘ng and their rolein sequentialanalysis．

の扉Co〃勿〟Jオ乃＆JEEg伽∫5（1986）84−89．

A犯乃αね〆5如才ゐfZcs，5（1977）316−329．

［22］Goyal，A．and S．S．Lavenberg：Modeling and［26］Serfling，R．J．：4坤7VXimation77teo7BmS qfMather

analysis of computer system availabilitylBMJqf
兄彷甜雅滋α肌7β甜以東w怖扉 31，6（1987）651−664．

maticalStattstics．（JohmWiley＆Sons，1980）．
［27］Shedler，G．S．：R留ene7diue Stochastic Simulation．

（AcademicPress，1993）．
［23］Goyal，A．，P．Shahabuddin，P．Heidelberger，Ⅴ．
Nicola and P．W．Glynn：A unified framework for［28］Reiman，M．I．andWeiss，A．：Sensitivityanalysis

Simulating Markovian models of highly dependable

for simulations via likelihood ratios．

systems．1EEE77uns．on Con4）ute73，C−41（1992）36−

斤gsgα作ゐ，37（1989）830−844．

［29］Rubinstein，R．Y．：The score function approach

51．

［24］Law，A．M．and W．D．Kelton：Simuhltion Model− for sensitivity analysis of computer simulation
ing and Anab）Sis，77ziYd Edition．（MaGraw−Hill，

models．Mathematics and Coク74）uterSinSimuhltion，28

1999）．

（1986）351379．

［25］Sen，P．K．：Someinvarianceprinciplesrelatingto

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

（17）柑1