マルコフ連鎖の数値解法のソースを表示

'''【まるこふれんさのすうちかいほう (numerical method for Markov chain)】'''

=== 概要 ===

マルコフ連鎖の定常分布や過渡的分布を数値的に計算するための方法. 定常分布の計算は線形方程式系を解く問題に帰着され, 解法として消去法や状態縮約法などの直接法, ヤコビ法, ガウス・ザイデル法, 状態縮約/非縮約法などの反復法がある. また, 構造化されたマルコフ連鎖に対しては行列幾何形式解を利用する方法も有力である. 一方, 過渡的分布に対しては, 離散時間ではべき乗法, 連続時間ではランダム化を利用する方法などがある.

=== 詳説 ===

　[[マルコフ連鎖]]をマルコフ連鎖の[[マルコフ連鎖の数値解法|数値的]]に解析する際の中心的な対象は[[定常分布]]である. [[状態空間|有限状態空間]] <math>{\mathcal S}=\{ 1, 2, \ldots, N \}\, </math> 上の[[既約]]で非周期的な (つまり[[エルゴード的マルコフ連鎖|エルゴード的]]) マルコフ連鎖を考え, その[[推移確率行列]]を<math>\boldsymbol{P}=(p_{ij})\, </math>, 定常分布を <math>\boldsymbol{\pi}=(\pi_1, \pi_2,\ldots,\pi_N)\, </math> とする. [[一様化]]により, 連続時間マルコフ連鎖の定常分布は, 離散時間マルコフ連鎖の定常分布として計算できるので, 以下では離散時間マルコフ連鎖に限定して考える. 

　エルゴード的なマルコフ連鎖では, 定常分布は



<table align="center">
<tr>
<td width="100">(平衡方程式)　</td>
<td><math>\pi_j = \sum_{i=1}^N \pi_i p_{ij}, \quad j=1, 2, \ldots,N, \qquad
\, </math></td>
<td  width="100" align="right"><math>(1)\,</math></td>
</tr>
<tr>
<td width="100">(正規化条件)　</td>
<td><math>\sum_{j=1}^N \pi_j = 1 \, </math></td>
<td  width="100" align="right"><math>(2)\,</math></td>
</tr>
</table>


を満たす一意の解として与えられる (式(1)) の解は定数倍に関して一意でないため, 式 (2) で正規化する). したがって, 定常分布の計算は, 原理的には線形方程式系を数値的に解く問題に帰着される. 状態数 <math>N\, </math> が大きくなければ, 消去法や[[状態縮約法]]などの直接法 (反復計算を伴わない方法) でも解を求めることは可能だが, 一般にマルコフ連鎖によるモデル化はモデルが複雑になるに従って状態数が急激に増加する傾向があるため, そのような場合は計算精度などを考慮して反復法を用いることが多い. 


'''ガウス・ザイデル法'''　反復法では, 反復回数 <math>k \to \infty\, </math> のときに定常分布 <math>\boldsymbol{\pi}\, </math> に収束するような近似値の列 <math>\boldsymbol{\pi}^{(k)} = (\pi_1^{(k)}, \pi_2^{(k)}, \ldots,\pi_N^{(k)})\, </math> を構成する. 例えば, ヤコビ法 (Jacobi method) では, 適当な初期分布 <math>\boldsymbol{\pi}^{(0)}\, </math> からスタートして


<center>
<math>
  \pi_j^{(k)} = \sum_{i=1}^N \pi_i^{(k-1)} p_{ij}, \quad
  j=1, 2, \ldots,N
</math>
</center>


によって分布列 <math>\boldsymbol{\pi}^{(k)}\, </math> を構成し, <math>\boldsymbol{\pi}^{(k-1)}\, </math> と <math>\boldsymbol{\pi}^{(k)}\, </math> が十分近くなった時点で収束したと判断する. エルゴード的なマルコフ連鎖に対しては, ヤコビ法は計算誤差を除けば必ず収束するが, 一般に大きな <math>N\, </math> に対してはあまり収束は速くない. これに対して, [[ガウス・ザイデル法]] (Gauss-Seidel method) では


<center>
<math>
  \pi_j^{(k)} = \frac{ \sum_{i=1}^{j-1} \pi_i^{(k)} p_{ij}
              + \sum_{i=j+1}^{N} \pi_i^{(k-1)} p_{ij} }{1-p_{jj}},
  \quad j=1, 2, \ldots,N
</math>
</center>


によって分布列を構成する. この方法では, <math>k\, </math> 回目の反復で既に更新されている値を逐次利用するため, ヤコビ反復法に比べると一般に収束が速くなることが多い. また, 推移確率行列がブロック構造を持つ場合には, ブロックごとに更新された値を利用する[[ブロックガウス・ザイデル法]] (block Gauss-Seidel method) も有効である. さらに収束を加速する手段として[[過剰緩和法]] (overrelaxation method) の利用がある. 過剰緩和法では, 緩和 (または加速) 係数を <math>\omega\, </math> として


<center>
<math>
  \pi_j^{(k)} = \frac{ \omega \sum_{i=1}^{j-1} \pi_i^{(k)} p_{ij}
              + (1-\omega) \pi_j^{(k-1)} p_{jj}
              + \omega \sum_{i=j+1}^{N} \pi_i^{(k-1)} p_{ij} }{1-p_{jj}},
  \quad j=1, 2, \ldots,N
</math>
</center>


によって <math>\boldsymbol{\pi}^{(k)}\, </math> を計算する. <math>\omega>1\, </math> のときには, 外挿により <math>\boldsymbol{\pi}^{(k-1)}\, </math> から <math>\boldsymbol{\pi}^{(k)}\, </math> を計算しており, 適切な <math>\omega\, </math> を選ぶことで収束を加速することが可能となる. なお, ガウス・ザイデル系の方法では, 初期分布 <math>\boldsymbol{\pi}^{(0)}\, </math> が (2) を満たしていても, 途中の計算でこの制約が満たされなくなるため, 計算の最後に (2) が満たされるよう正規化することが必要である. 


'''状態縮約/非縮約法'''　一方, 複数の状態をまとめて1つの状態と見なした状態数の少ない確率過程に対して反復計算を行う方法に[[状態縮約/非縮約法]] (aggregation/disaggregation method: AD法) がある. 例えば, 状態空間を <math>L\, </math> 個の部分空間 <math>{\mathcal S}_1, \ldots, {\mathcal S}_L\, </math> に分割し, <math>{\mathcal S}_{\alpha}\, </math> には <math>d_\alpha\, </math> 個の状態 <math>(\alpha,1), \ldots, (\alpha,d_\alpha)\, </math> が含まれる場合を考え, 推移確率を <math>\boldsymbol{P}=( p_{(\alpha,i)(\beta,j)} )\, </math>, 状態 <math>(\alpha,i)\, </math> の定常確率を <math>\pi_{\alpha,i}\, </math>, 部分空間 <math>{\mathcal S}_\alpha\, </math> の定常確率を<math>\textstyle \tau_\alpha=\sum_{i=1}^{d_\alpha} \pi_{\alpha,i}\, </math> とする. いま, <math>k-1\, </math> 回の反復で近似値 <math>\pi_{\alpha,i}^{(k-1)}\, </math>, <math>\tau_\alpha^{(k-1)}\, </math> が求められているとしよう. <math>k\, </math> 回目の反復計算のうち, まず縮約フェーズでは, 部分空間 <math>{\mathcal S}_\alpha,\; \alpha = 1,\ldots,L\, </math> をそれぞれ1つの状態 <math>s_\alpha\, </math> に縮約した <math>L\, </math> 状態の確率過程を考え, それをマルコフ連鎖と見なして (特殊なケースを除いて縮約した確率過程はマルコフ連鎖とならない) 推移確率, 例えば <math>s_\alpha\, </math> から <math>s_\beta\, </math> への推移確率を<math>\textstyle q_{\alpha,\beta}^{(k)}=\sum_{i=1}^{d_\alpha} \sum_{j=1}^{d_\beta}\pi_{\alpha,i}^{(k-1)} p_{(\alpha,i)(\beta,j)} / \tau_\alpha^{(k-1)}\, </math> によって定める. このマルコフ連鎖 <math>\boldsymbol{Q}^{(k)}=(q_{\alpha,\beta}^{(k)})\, </math> の平衡方程式を解いて, 更新された定常確率 <math>\tau_\alpha^{(k)},\; \alpha=1,\ldots,L\, </math> を求める. 次に非縮約フェーズでは, 1つの着目した部分空間はそのままで他のすべての部分空間を1つの状態に縮約した確率過程を近似的にマルコフ連鎖と考える. 例えば, 部分空間 <math>{\mathcal S}_\alpha\, </math> に注目した場合には, <math>{\mathcal S}_\alpha\, </math> 内の推移確率は元のままで, <math>{\mathcal S}_\alpha\, </math> 内の状態 <math>(\alpha,i)\, </math> から縮約された状態への推移確率は <math>\textstyle \sum_{\beta \ne \alpha}\sum_{j=1}^{d_\beta} p_{(\alpha,i)(\beta,j)}\, </math>, 逆に縮約された状態から <math>(\alpha,i)\, </math> への推移確率は <math>\textstyle \sum_{\beta \ne \alpha}   \sum_{j=1}^{d_\beta}\pi_{\beta,j}^{(k-1)} p_{(\beta,j)(\alpha,i)}/(1-\tau_{\alpha}^{(k)})\, </math> で与えられるマルコフ連鎖を考え, その定常分布を計算し <math>\pi_{\alpha,i}^{(k)}, \; i=1,\cdots,d_\alpha\, </math> を得る. この計算を, 注目する部分空間を <math>{\mathcal S}_1\, </math> から <math>{\mathcal S}_L\, </math> まで変えながら行えば, 更新された定常確率を求めることができる. この縮約/非縮約の手続きを, 値が収束するまで反復すればよい. 


'''無限状態と過渡的分布'''　状態数が無限のマルコフ連鎖に対しては, 状態空間を適当な有限サイズで打ち切って数値計算を行うが, 打ち切るサイズによって計算時間と計算精度の間にトレードオフが生じるので注意が必要である. 構造が入っている場合 (後述) は, 上の方法を用いるにしてもその構造をうまく利用することによって, 少ない計算量で精度良い解が計算できることが多い. 

　定常分布に比べると, 過渡的分布 (各時点における推移確率) の計算方法はそれほど多くないが, 離散時間マルコフ連鎖に対してはべき乗法, 連続時間マルコフ連鎖に対してはランダム化を利用する方法などが知られている. 


'''構造化されたマルコフ連鎖'''　確率モデル, 特に待ち行列モデルから派生するマルコフ連鎖には, 何らかの構造を持つものが多いため, その構造を利用した数値計算法が開発されている. 代表例として, [[相の方法|相型待ち行列]]に対する[[行列幾何形式解]]を考えよう. [[到着過程]]や[[サービス時間分布|サービス過程]]}に[[マルコフ型到着過程]]や[[相型分布]]を導入することで, 広い範囲の待ち行列モデルは準出生死滅過程 (quasi-birth-and-death process) を含むGI/M/1型, あるいはM/G/1型マルコフ過程などの構造化されたマルコフ連鎖で表現することができる. このうち, GI/M/1型マルコフ連鎖は, レベル <math>n\; (=0,1,\ldots)\, </math> と相 <math>i\;(=1,\ldots,d)\, </math> の組 <math>(n,i)\, </math> によって状態が表されるマルコフ連鎖で, 1回の[[推移]]では高々1つ上のレベルまでしか推移せず, またレベル <math>n\, </math> の状態からレベル <math>m\; (m\le n+1)\, </math> の状態への推移確率 (または推移速度) がレベルの差 <math>m-n\, </math> と各状態の相によって決まる性質を持っている. レベル <math>n\, </math> の状態の定常確率ベクトルを <math>\boldsymbol{\pi}_n=(\pi_{n,1}, \ldots, \pi_{n,d})\, </math> で表すと, 行列幾何形式解より, [[公比行列]] <math>\boldsymbol{R}\, </math> を用いて


<center>
<math>
  \mathbf{\pi}_n = \mathbf{\pi}_0
  \mathbf{R}^n, \quad n=1,2,\ldots
</math>　　　　　<math>(3)\,</math>
</center>


と表される. <math>\boldsymbol{R}\, </math> は推移確率行列の要素を係数とする非線形行列方程式の非負最小解として与えられ, 逐次代入法などで計算することができる. また <math>\boldsymbol{\pi}_{0}\, </math> は境界条件に相当する線形方程式を解いて求められる [2]. この方法は, 本来無限次元の定常分布を有限次元のベクトルと行列で表せるという特徴を持つが, 高速化のためには <math>\boldsymbol{R}\, </math> の計算方法がポイントとなる. 

　なお, M/G/1型マルコフ連鎖は行列幾何形式解を持たないが, やはりその構造を利用したさまざまな方法が考えられている [2]. 



----
'''参考文献'''

[1] D. P. Heyman and M. J. Sobel (eds.), 伊理, 今野, 刀根監訳, 『確率モデルハンドブック』, 朝倉書店, 1995. 

[2] M. F. Neuts, ''Matrix Goemtric Solutions in Stochastic Models - An Algorithmic Approach'', Johns Hopkins Univ. Press, 1981.

[3] M. F. Neuts, ''Structured Stochastic Matrices of {rm M/G/1} Type and Their Applications'', Marcel Dekker, 1989.

[4] W. J. Stewart (ed.), ''Numerical Solution of Markov Chains'', Marcel Dekker, 1991.

[5] W. K. Grassmann (ed.), ''Computational Probability'', Kluwer Academic Publishers, 2000.

[[category:確率と確率過程|まるこふれんさのすうちかいほう]]