質問<3615>2007/9/21from=ひで次郎「マルコフ連鎖」
状態{-2,-1,0,1,2}上のランダム・ウォークにおいて、時刻nで粒子が状態i(i=-1,0,1)にあれば
時刻n+1でi-1またはi+1に、時刻nで粒子が状態-2または2にあれば時刻n+1には等確率で、
-1,0,1の1つに移動するものとする。このマルコフ連鎖{Xn}に対して次の問に答えよ。
- 推移確率行列Pを求めよ。
- を求めよ。
宜しくお願いします。
★希望★完全解答★
まずは翻訳していってみましょう。
状態をI={-2,-1,0,1,2}とします。
時刻nで粒子が状態i(i=-1,0,1)にあれば時刻n+1でi-1またはi+1に(等確率で)移動する
これをまずは馬鹿正直に数式に翻訳します。
例えば時刻nで粒子が状態i=−1にいて、時刻n+1でi−1(つまりこの場合−2)に移動する、と。
これを条件付き確率の表記法に準じて書くと、
と表現出来ます。同様に、時刻nで粒子が状態i=−1にいて、時刻n+1でi+1(つまりこの場合0)に移動する場合は、
と書き表します。
この作業を取り敢えずは全部の条件に従って書き出さないとならない。じゃないと推移確率行列Pを作る/作らない以前の話になっちゃうんです。
題意に沿って次々と条件付き確率を書き出すと
となりますね。
時刻nで粒子が状態-2または2にあれば時刻n+1には等確率で、-1,0,1の1つに移動する
これもやるこたあ同じです。
書き出していきます。
まず、問題を解く、解かない以前にこういう材料揃えくらい出来てないとお話しになりません。
話はそれから、です。
推移確率行列Pを求めよ。
マルコフ連鎖の行列Pを使った表現は、n+1期の確率ベクトルを、n期の確率ベクトルをと表現すると、
と書き表せます。
詳しくは質問3590を見てほしいのですが、この問題の場合、それぞれの確率ベクトルの内訳は、状態が{-2,-1,0,1,2}の五つしか無い以上、
だと言う事です。
あとは、線形代数の知識を駆使して、左辺の要素と右辺の要素が一致するように(質問3590のエビデンス参照)先ほど求めた要素をただ並べてやれば良い。
具体的には、例えば、
を表してるので、
となるので、推移確率行列の1行目は
同様の作業を全てのエビデンスに関して行うと、マルコフ連鎖の行列表記は
となり、結果推移確率行列Pは
となることが分かります。
(成分を縦に足していくと全て1になる、と言うことに留意。)を求めよ。
今までで分かることは
だと言うことです。
また、
も成り立ってる、と言う事です(と言うより、マルコフ連鎖とはそう言う要請です)。
すなわち、
も成り立っています。
これを一般化すると、
も成り立っていますね。
すなわち言い替えるとが一体何なのか、と言うのが題意であるワケです。が分かればイイ。
さて、オーソドックスな線形代数の知識を使った手段を用いると、PをPの固有ベクトルを用いて対角化し、のカタチを予測すれば一丁上がりです(質問3566参照)。と・こ・ろ・が。
Pのランクのせいで、この対角化は上手く行かないんですね。そ・こ・で・・・・・・。
数学的にはインチキなんですがを対角化してみます。つまり、極限値を取る限り、
を計算しようと結果は変わんないんです(笑)。何故なら、このテの極限値絡みの問題と言うのは、大体のトコ振動しようと何だろうと「収束値を求めよ」と言う質問のパターンが多い、と言うのがまずは一つ。加えて、マルコフ連鎖そのものが実用的にはある意味、「定常分布に収束するのを見越した」部分で設定してるんですね。つまり、これだけややこしい計算しておいて「発散しちゃった」じゃあまり実用的じゃないんですよ(笑)。よって、その手の「統計学的テクニック」を鑑みる限り、結果上の式を使って議論しても結果はそんなに変わらないのです。
さて、本来はシコシコ手計算でやるトコですが、メンド臭いんで、Maximaを使います。まあ、本当だったら練習代わりに手計算して欲しいんですが、「確率の問題」として考える以上、「線形代数」は計算するのには必要なテクニックですが、それやっちゃうと本末転倒なんですね。よってここからはパソコンに丸投げします。Maximaを起動して下さい。
Maximaで問題の推移確率行列Pを定義する。
行列入力の方法に関しては質問3566参照。
Pの定義が確認出来たらを定義する。次のように[INPUT:]に入力すれば良い。
P2:P^^2
なお、行列の冪乗はカレット(^)を二個並べて計算させる。
P2の定義が確認出来たら、eigenvectorsコマンドを使っての固有値、固有ベクトルを求める。具体的には[INPUT:]に次のように入力する。
eigenvectors(P2)
eigenvectorsコマンドに付いても質問3566参照。
固有値、固有ベクトルの計算が終わったら、コピペを上手く使って、Pの固有ベクトルを列成分に持つ行列Vを生成する。
具体的には、matrix(行列)コマンドを使って
V:matrix([1,(sqrt(3)-3)/3,-(2*sqrt(3))/3,(sqrt(3)-3)/3,1]
,[1,-(sqrt(3)+3)/3,(2*sqrt(3))/3,-(sqrt(3)+3)/3,1]
,[1,2,8/3,2,1],[1,0,0,0,-1],[0,1,0,-1,0])
となるが、括弧()の中身は固有ベクトルの計算結果をコピペしてくれば良い。
行列Vの生成が確認出来たら、今度はそれの転置行列Uを作る。
具体的には[INPUT:]欄にtranspose(転置)コマンドを利用して次のように入力する。
U:transpose(V)
なお、転置行列にする理由は、対角化の為の行列は固有ベクトルを列で並べたモノじゃないといけないから。
行列Vは固有ベクトルを行で並べたモノだったので、ここで数学的要請により転置させた。
Uの定義が確認出来たらUの逆行列IUをinvert(ひっくり返す)コマンドを使って生成する。
具体的には[INPUT:]欄に
IU:invert(U)
と入力する。
ものすごい出力が出てくるのでビビる(笑)。
逆行列IUは少々複雑な出力なんで見やすいようにratsimpコマンドで単純化する。
具体的には[INPUT:]欄に、次のように入力してIUを再定義してやる。
IU:ratsimp(IU)
逆行列IUの出力が整理される。
また、[Enter matrix ...]を使って、先ほど求めた固有値をn乗したモノを左上から右下に向かって順番に置いた行列lambdaを定義する。
行列lambdaの生成が確認出来たら、[INPUT:]にU、lambda、IUの三つの行列の積を計算させるように入力する。具体的には
U . lambda . IU
と入力する。
なお、行列同士の掛算の場合、アスタリスク(*)の代わりにピリオド(.)を用いる。
三つの行列の積、の答が出力される。
さて、Maximaが出力した推移確率行列でにあたる位置(最下段の3列目)の数字を調べてみると、
となっています。すなわち、題意に沿うと、
が何なのか、って事ですね。
ところで、nが絡む項は分子と分母を比べれば分子の方が小さい事が分かるでしょうか?従ってnが増加すると第1項、第2項共にゼロへと収束します。
故に、
ってのが答えになります。
以上です。
Maxima Manual: 行列と線形代数