数式処理ソフト　DERIVE（デライブ） de ドライブ

８８．マルコフ過程（1）（非対称型道中双六：宿場数=3～9）

はじめに

「第87回『道中双六の問題（補筆3）（今後の課題）』の末尾で、『マルコフ過程』について、今後、勉強することを『課題5』として挙げました。
　そこで、今回、『確率過程超！入門』（松原望著：東京図書：2011年10月第1刷）を入手して、（少し）勉強してみました。
　
　これにより、同書に記載の『推移確率行列』の考え方を使って、以下に述べましたように、『非対称型道中双六』の宿場数＝3の場合の平均日数の厳密解を求めることができました。
　得られた結果は、平均日数＝（4p^2－2p＋2）/p^3、となりました。
　従来は、平均日数に関する確率 pに関して閉じた関数形を求めることができず、2m+4 日目に上がる確率を f（2m+4）とするとき、mをm=0～の整数として、
　f（2m+4）=（1－p）^m p^3 Σ（A+B+C=m）（p^（m－A）Comb（A+B，A） Comb（B+C，B））、
　ただし、A、B、Cは、0以上の整数、は、得られていましたので、
　平均日数=Σ（m=0～n）（2m＋4）f（2m+4）、n⇒大として、数値計算で求めていたものです。
　また、宿場数＝4並びに5の場合も、宿場数＝1～3までの結果と数値計算を併せて検討したところ、これらの平均日数の厳密解（の予想式）を得ましたので、それらも記載しました。

　なお、今回述べた方法だけでは、五十三次（宿場数＝53）のように大きい宿場数の場合の一般解を得ることは、難しいと思われますが、以前に考えていたよりも、平均日数の関数形が簡単な可能性が高くなったことにより、今後、より広い範囲で厳密解が求められるのではないかと考えています。　
　※宿場数＝6～9について、予想式を追加しました。（2018/4/6 追記）
　目次へ戻る

（1）推移確率行列

非対称型道中双六と対称型道中双六

「これまで、何回も説明してきたので、読者の方がいい加減うんざりされるじゃろう。
　ここでは、必要なことだけをごく簡単に説明しよう。
　詳しくお知りになりたい方は、第87回の『道中双六の問題の概要』などをご覧ください。
　
　１．駒は、振り出しから、確率1で、1番目の宿場に進むことができる。
　２．駒は、各宿場で、確率 pで、次の宿場に進み、確率 (1－p)で1つ戻るものとする。
　３．駒が最後の宿場から上がりに進んだときは、双六上から消える。（上がりから戻ることはない）
　４．駒どおしの関係は、無いものとする。
　５．振り出しから宿場まで、各宿場間、最後の宿場と上がりまでの所要時間は、同一とし、1（日）と勘定する。
　というルールに従って、駒を進めるとき、平均何日で上がれるかという問題じゃ」

「えーと、読者の方に、わたしから説明するけど、
　おじぃさんが『非対称型道中双六の問題』と、わざわざ、『非対称型』と断っている訳はね、
　前回で『対称型道中双六』について、今後の課題としてあげたからなの。
　ちなみに、対称型というのは、こんな感じね。
　
　どこが違うかっていうと、対称型でも駒は、振り出しからスタートするけど、上がり方向だけでなく、下がり方向にも戻ることができるのね。
　それと、非対称型では、振り出しから進む確率は、1だったんだけど、対称型では、進む確率は、p、戻る確率は、1－p、というように振り出しが普通の宿場と同じになっている点が違うのよ」

「お銀さんかい。
　上掲の本のおかげで、ようやく、胸のつかえがとれましたな。
　宿場数=3の場合の平均日数（の厳密解）が分かりましたぞ」
　目次へ戻る

マルコフ過程

「あんたは、水戸黄門かいな。
　それより、マルコフ過程の『マルコフ』というのは、人の名前かしら？」

「おお、ともちゃんかい。
　マルコフ過程のマルコフとは、ロシアの数学者、アンドレイ・マルコフ氏（1856年～1922年）のことじゃ。
　同氏の生涯・業績のあらましは、ウィキペディアを参照して欲しい。
　上掲の本に書かれているように、
　『マルコフ過程』とは、現在のシステムの状態が一つ前の時刻（連続時間では、微少時間前）における状態にのみ陽に依存する場合のシステムの推移（遷移）を（調べることを）指すようじゃ。
　なお、時間がとびとび（離散的）の場合を特に『マルコフ連鎖』とも言う。
　ただし、広い意味では、1つ前だけでなく、それ以前の状態にも陽に依存するケースも含めることがあるそうじゃが、ここでは、狭い意味のマルコフ過程を考えていこう」

「その、システムの状態というのは、道中双六で言うと、何を指すのかしら？」

「道中双六の問題で、差分方程式を立てたことがあったのを覚えておるじゃろう。
　それを振り出しの位置をゼロからではなく1から数えることにして、振り出しを1とし、最後の宿場をc、上がりをc＋1、とする。
　ここで、c=宿場数＋1、である。
　例えば、宿場数が4の場合（c=5）では、振り出しを1、最初の宿場を2、最後の宿場を5、上がりを6とする。
　さて、一般に、日数をtで表し、t日経過後の各位置にいる駒の確率密度をf（j，t）とすると、t=0、1、2、・・に対して、
　（1）　f(1，0)＝1、初期条件
　（2）　f(1，t)＝(1－p)×f(2，t－1)、振り出し
　（3）　f(2，t)＝f(1，t－1)、最初の宿場
　（4）　f(c，,t)＝p×f(c－1，t－1)、最後の宿場、
　（5）　f(c＋1，,t)＝p×f(c，t－1)、上がり、
　（6）　f(j，t)＝(1－p)×f(j+1，t－1)＋p×f(j－1，t－1)、j= 2～c－2の宿場、
　が成立している。
　ただし、
　宿場数が1（c=2）では、最初の宿場が最後の宿場であるため、（4）及び（6）を使わない、
　宿場数が2（c=3）では、（6）を使わないことに注意するようにな」

「DERIVEで平均日数を計算する関数を作るときには、こんな図を書いて、調べたわね。
　
　上の図は、第85回の『道中双六の問題（補筆2）』に表れたもので、余分な説明を省いたものよ」

「そうじゃったな。
　今回は、確率密度関数を次のようなベクトルで表す。
　初期状態　t=0、
　U0=［f(1，0)，f(2，0)，f(3，0)，f(4，0)，f(5，0)］、これは、c=4、宿場数が3の例、
　一般には、U0=［f(1，0)，f(2，0)，f(3，0)，f(4，0)，f(5，0)，・・・，f(c，0)，f(c＋1，0)］、
　一方、t=m日後の状態は、
　Um＝［f(1，m)，f(2，m)，f(3，m)，f(4，m)，f(5，m)］、これは、c=4、宿場数が3の例、
　一般には、Um=［f(1，m)，f(2，m)，f(3，m)，f(4，m)，f(5，m)，・・・，f(c，m)，f(c＋1，m)］、
　などと書くことにしよう」

「なるほど。
　ベクトルと行列を使おうという訳ね」
　目次へ戻る

推移確率行列

「その通りじゃ。
　U0から1日後の状態は、宿場数が3（c=4）の例でどうなるか考えてみよう。
　なお、状態を表すベクトルは、上の例では、横ベクトルとしているが、左から行列を掛ける普通の記法に合わせると本来は縦ベクトルする必要がある。
　しかし、HP上では、縦ベクトルを書きにくいので、横ベクトルのように書いている。
　すると、
　0日後では、U0=［1，0，0，0,，0］、初期状態
　1日後では、U1=［0，1，0，0,，0］、
　2日後では、U2=［1－p，0，p，0,，0］、
　3日後では、U3=［0，(1－p)^2，0，p^2,，0］、
　4日後では、U4=[(p － 1)(p^2 － 1), 0, － p(2p^2 － p － 1), 0, p^3]、最初の上がりの日、
　などとなる。
　ここで使われているのが、推移確率行列 A じゃ。
　U1=A×U0、などと計算する。
　Aは、DERIVE画面では、こんなかんじじゃな。
　
　行列の要素をA（i，j）で表し、状態ベクトルをU、要素をU(k)、k=1～c＋1、で表すとき、
　A*Uのk番目の要素は、
　Σ（j=1～c＋1）A（k，j）×U（j）、であることは、行列×ベクトルの計算の基本じゃ。
　上で掲げた推移確率ベクトルは、これを見ると、
　j 列目の要素から i 行目の要素に至る確率とする必要がある。
　A（2，1）＝1、
　　1番目（振り出し）から2番目（最初の宿場）には、確率1で進む、
　A（j，j＋1）＝1－p、j＝1～c－1、
　　j＋1番目からj番目には、1－pで戻る
　A（j＋2，j＋1）＝p、j＝1～c、
　　j＋1番目からj＋2番目には、pで進む
　A（c＋1，c＋1）＝1、
　　c＋1番目の要素は、以前の値を保持（してそれに加算する）。
　わかりやすく説明を加えると、次のような図解となる」

　

「すると、2m＋4日後の状態は、mを0以上の整数として、
　A^（2m＋4）×U0 で与えられるということね。
　でも、Aは、行列よね。
　そのべき乗の計算は、面倒でしょ？」

「さすがは、ともちゃん。
　まさに、一を聞いて十を知る女子(おなご)じゃな。
　そこで、行列の固有値と対角化を使おうということじゃ。
　なお、一つご注意を。
　上の状態では、ベクトルの要素の和が常に1となっているが、上がりを状態ベクトルのc＋1個目の要素としているためじゃ。
　たとえば、宿場数が3の場合では、c＝4であり、上の説明では、状態ベクトルと行列の次元は、c＋1＝5となる。
　しかし、今回のテーマである平均値等の計算では、最後の宿場までであればよい。
　これは、上がりの一つ手前の要素にpを乗じたものが2m＋4日後に上がる勘定となることから、
　3日後の最後の要素の値（p^2）にpを乗じて、p^3を作れば、これが最初に、上がりに加算される量となる。
　従って、平均日数＝Σ（2m＋4）×p×（2m＋3日後の状態ベクトルの最後の要素）、で計算できる。
　すなわち、ベクトルと行列の次元は、c＝4でいいことになる。
　上の図解を使えば、次のようになる。

　

　次節で、宿場数1の場合で具体的に説明しよう」　　
　目次へ戻る

宿場数=1（c=2）の平均日数

「前節の最後の注意によれば、状態ベクトルは、2次元で、［振り出し，宿場］、
　すなわち、初期状態ベクトル U0＝［1，0］となる。
　では、ともちゃん、推移確率行列 A は、どうなるかな」

「宿場数が1の場合は、平均日数が2/p となった例ね。
　ちなみに、上がるまでの日数は、2m＋2、m=0～の整数、が必要。
　さて、推移確率行列をAと書けば、
　
　となるわ」

「では、推移確率行列 Aの固有値と固有ベクトルを求めて欲しい」

「行列の固有値と固有ベクトルについては、第63回『再帰関数（2）（数列の母関数、固有値のべき乗法）』で取り上げたことがある。
　固有値は、数式処理ソフト DERIVEでは、
　EIGENVALUES（行列，（変数名））で求めることができる。（変数名は、省略可）
　また、固有ベクトルは、
　EXACT_EIGENVECTOR（行列，厳密な固有値）で縦ベクトルとして、求めることができる。
　今回の固有値は、EIGENVALUES(A))により、([√(1 － p), － √(1－ p)])の2つであることが分かる。

　一つ目の固有ベクトルは、EXACT_EIGENVECTOR(A, √(1 － p))から、
　
　と求められた。
　2つめも同様に、EXACT_EIGENVECTOR(A, －√(1 － p))から、
　
　と求められる」

「2つの固有ベクトルから、推移確率行列 Aを対角化する変換行列をRと書くとき、対角化された行列をBとすると、
　B＝R^(-1) ×A×R、と求められる。（R^(－1) は、Rの逆行列）
　この変換行列Rは、固有ベクトルを並べれば、得られるが、ベクトル同士を結合する必要がある。
　そのために、DERIVEでは、ファイルメニューから、ユーティリティファイルを読み込む必要がある。
　Program Files (x86)\TI Education\Derive 6\Math\VectorMatrixFunctions.mth、
　これにより、APPEND（ベクトル，ベクトル，・・）でベクトルを結合できる。
　※APPEND関数は、DOS版でも同様の機能を持つ関数があるが名称が変わったので注意。

　
　ただし、元になっている固有ベクトルの表示が縦ベクトルなので、転置（右上のアポストロフィ）により横ベクトルに変換している。
　
　こうして得られた行列を再度、転置すると、下図のRが得られる。
　
　では、これが変換行列になっていることを確かめてご覧」

「Rの逆行列は、R^(－1)なので、
　R^(－1)×A×Rを作ると、
　
　となって、確かに、対角化されていることが分かるわ」

「これで準備が整ったな。
　まず、対角化された行列 Bの（2m＋2－1）＝（2m＋1）乗を求める。
　
　これが、B^(2m＋1)、だな。

　次に、R×B^(2m＋1)×R^－1、を作る。
　

　最後に、［1，0］’を右から掛けると、（初期状態ベクトルU0を転置して縦ベクトルとしている）
　
　これは、［振り出し，宿場］なので、最後の要素は、確率pで上がりに進む。
　そこで、最後の要素を取り出して、p×（2m＋2）を掛けて、m=0～nまで加えると、平均日数が得られる。

　平均日数＝Σ（m=0～n）（2m＋2）×p×（1－p）^m、
　
　ここで、n→∞とすれば、
　平均日数＝2/p、の式が得られるという訳じゃ」

「なるほど。
　これまでの回では、極限をとる前の式は、
　平均日数＝2/p － 2(1 － p)^(n + 1)(np + p + 1)/p 、だったので、
　上の結果との差し引きすると、下図のようにゼロになるため同じ結果であることが分かった」
　

「DERIVEを使う際、2つほど注意するところがある。
　変数 pの範囲を0＜p＜1、
　変数mの範囲を正の整数とする。
　これは、「入力」メニューから変数の変域から行う。
　
　上の図は、mを正の整数としているところじゃ。
　さて、次節では、宿場数が2の場合を計算してみよう」　
　※ 上記のように、n次正則行列を対角化可能なのは、一次独立なn個の固有ベクトルが得られる場合です。（2018/4/9 追記）
　目次へ戻る

（2）宿場数＝2　（c=3）の平均日数

前節と基本的な内容は同じなので、項目毎に結果のみを示すことにします。

初期状態ベクトル、推移確率行列 A

初期状態ベクトルは、3次元なので、［振り出し，宿場1，宿場2］、となり、
　転置をつけた状態で、Ｕ０＝([1, 0, 0]`)、と表せます。
　推移確率行列 A は、上がりを省略した状態では、3次行列となり、
　
　となりました。
　目次へ戻る

固有値

Aの固有値は、3つあり、
　([0, √((p + 1)(1 － p)), － √((p + 1)(1 － p))])、でした。
　
　目次へ戻る

固有ベクトル

3つの固有値に対応した固有ベクトルは、次の通りです。
　ゼロに対しては、
　
　残りの2つに対しては、
　
　が得られました。
　目次へ戻る

変換行列 R

B＝R^-1 ×A×R、ここで、Bは、対角化された推移行列とする変換行列は、
　
　と求められました。
　確認のため、R^-1 ×A×R を計算すると、
　
　と確かに対角化されていることが分かります。
　目次へ戻る

平均日数

平均日数＝Σ（m=0～n）（2m＋3）×p×（（R×B^(2m+2)×R^(-1)）×U0の最後の要素）、
　上の式で、宿場数1の例と比較して、2m＋1 → 2m＋2、2m＋2→ 2m＋3、と変わっていることに注意してください。
　まず、B^(2m+2)は、
　
　となる。
　（（R×B^(2m+2)×R^(-1)）×U0の最後の要素は、
　
　となるので、
　平均日数＝∑(p^2(2m + 3)(－ √((p + 1)(1 － p)))^(2m)/2 + p^2(p － 1)^m(2m + 3)(－p － 1)^m/2、
　
　これは、(－ (2np^2 + 3p^2 + 2)((p + 1)(1 － p))^(n + 1)/p^2)、であり、
　
　最後に、n→∞とすれば、
　平均日数＝1＋2/p^2、が得られました。
　これは、以前に求めていた平均日数の式と等しいことが分かります。　
　目次へ戻る