マルコフ連鎖の時間発展の数値計算

(1)

マルコフ連鎖の時間発展の数値計算

樋口さぶろお

龍谷大学理工学部数理情報学科

計算科学☆実習

B L07(2016-05-23 Mon)

最終更新: Time-stamp: ”2016-05-23 Mon 18:41 JST hig”

今日の目標

マルコフ連鎖の分布の時間発展を計算するプロ

グラムが書ける

(2)

マルコフ連鎖の時間発展

L05-Q3

TA Prob and Sol:

マルコフ連鎖の時間発展

次の推移確率行列を持つ

,

状態空間

{1, 2}

上のマルコフ連鎖を考えよう

.

M =

(

₁

3

1

2

3

1

2

)

.

1

定常分布をすべて求めよう

.

2

初期分布

⃗

p(0) = (

1

0

)

のとき

⃗

p(t)

を求めよう

.

3

初期分布

_⃗

_{p(0) =}

1

2

(

1

)

のとき

⃗

p(t)

を求めよう

.

略解

(3)

マルコフ連鎖の時間発展 1

推移確率行列

_M

の固有値

_λ

₁

_{, λ}

_2,

固有ベクトル

_⃗

_u

₁

_{, ⃗}

_u

₂

を求めると

_,

λ

1

= 1

≥ λ

2

=

−

1 ₆

,

⃗

u

1

= (

3

4

) s, ⃗

u

2

=

(

₁

−1

)

s

(s

̸= 0).

固有値

λ = 1

の固有ベクトルである確率ベクトルは

1

7

(

3

4

)

のみであ

り

,

これが唯一の定常分布

.

2

_⃗

_u

₁

は確率ベクトルとなるように

_{s =}

1

7 ととる

. ⃗

u

2 は確率ベクトル

でないので適当に

s = 1

ととる

.

他の取り方をしても最終的には同

じ

⃗

p(t)

が求まる

.

⃗

p(t) =M

t

⃗

p(0)

=(U DU

−1

)

t

⃗

p(0) =

(

⃗

u

1

⃗

u

2

) (

λ

t 1

0 0 λ

t 2

) (

⃗

u

1

⃗

u

2

)

−1

⃗

p(0)

(4)

U

−1

⃗

p(0)

さえ計算すればいいわけだが

,

これを

(

a

_b

)

とおいてあとで

決める方が楽

.

⃗

p(t) =

(

⃗

u

1

λ

t1

⃗

u

2

λ

t2

)

(

a

_b

)

=a⃗

u

1

λ

t

1

+ b⃗

u

2

λ

t

2

.

これが

⃗

p(0)

と等しくなるように

,

連立方程式

⃗

p(0) = a⃗

u

1

+ b⃗

u

2 を解

いて

a, b

を決めると

a = 1, b =

4 ₇

.

⃗

p(t) =

₇

1

(

3 ₄

) +

4 ₇

(

₋₁

1

)

(

−

1 ₆

)

t

.

t

→ +∞

でも

⃗

p(t)

が確率ベクトルでありつづけることから

, a = 1

は計算しなくてもわかる

.

時間変化

.

| −

1 ₆

| < 1

なので

, ⃗

p(t)

→

1 ₇

(

3

4

)

(t

→ +∞)

(5)

マルコフ連鎖の時間発展 0 0.25 3/7 0.54/7 0.75 1 0 1 2 3 4 p(x,t) t p(1,t) p(2,t) 3

⃗

p(t) =

1 ₇

(

3 ₄

) +

₁₄

1

(

₋₁

1

)

(

−

1 ₆

)

t

(6)

L06-Q2

Quiz

解答

:

可約なマルコフ連鎖の定常状態

固有値は

λ = 1 (

重解

)

に対応する固有ベクトルは

,

(

₁

0

)

s +

(

₀

1

)

k

(s

̸= 0

または

k

̸= 0).

固有値

1

なので

,

線形独立な確率

ベクトルを

1

組選ぶと便利

(

他の選び方でも最終的な結果は変わらない

)

で

, ⃗

u

1

=

(

₁

0

)

, ⃗

u

2

=

1 ₂

(

₀

1

)

とすると

, s⃗

u

1

+ (1

− s)⃗u

2

(0

≤ s ≤ 1)

は定常

分布

.

固有値

λ =

1 ₃

に対する固有ベクトルは

,

(

₀

+1

−1

)

s(s

̸= 0).

確率ベクトルに

はなりえないので

,

適当に非零ベクトルをひとつ選んで

(

他の選び方でも

最終的な結果は変わらない

), ⃗

u

3

=

(

₀

+1

−1

)

.

⃗

p(0)

が一般の場合の時間発展は

⃗

p(t) = a⃗

u

1

t

+ b⃗

u

2

1

t

+ c⃗

u

3

· (

1 ₃

)

t

.

(7)

マルコフ連鎖の時間発展 1

_⃗

_{p(0) =}

1

2

(

₁

1

0

)

から

a, b, c

を定めて

,

⃗

p(t) =

1 ₂

⃗

u

1

+

1 ₂

⃗

u

2

+

1 ₄

⃗

u

3

· (

1 ₃

)

t

極限分布

⃗

p(

∞) =





1

2

1

4

1

4





は定常分布のひとつ

.

2

⃗

p(0) =

1

3

(

₁

1

)

から

a, b, c

を定めて

,

⃗

p(t) =

1 ₃

⃗

u

1

+

2 ₃

⃗

u

2

+ 0⃗

u

3

· (

1 ₃

)

t

極限分布

⃗

p(∞) =





1

3

1

3

1

3





は定常分布のひとつ

.

3

{1}

と

{2, 3}

が分かれた推移図

_.

すなわち

_,

可約なマルコフ連鎖で

ある

.

(8)

L06-Q3

Quiz

解答

:

マルコフ連鎖

1

推移確率行列

T

の固有値

λ,

固有ベクトルを求めると

,

λ = 1, ω, ω

2

(

₁

1

)

s,

(

₁

ω

2

)

s,

(

₁

ω

2

ω

)

s

(s

̸= 0)

定常分布は

λ = 1

に対する固有ベクトルで

,

確率ベクトルになるよ

うに

s

を定めると

, ⃗

u

1

=

1 ₃

(

₁

1

)

.

2

他の固有値

λ = ω, ω

2 は

,

|ω| = |ω

2

| = 1

を満たす

.

よって

,

一般には

極限分布は存在しない

.

極限分布が存在するのは

, ⃗

p(0) = ⃗

u

1

= ⃗

p(t)

のように最初から定常分

布であったときに限られる

.

(9)

L06-Q4

Quiz

解答

:

周期的なマルコフ連鎖の定常状態

1

固有値

λ = 1

の固有ベクトルは

, (

1

) s (s

̸= 0).

確率ベクトルになるように

s

を選んで

, ⃗

u

1

=

1 ₂

(

1

)

が定常分布

.

なお

,

固有値

−1

の固有ベクトルは

(

₋₁

1

)

s (s

̸= 0).

適当に

s

を選んで

, ⃗

u

2

=

(

₁

−1

)

とおく

,

一般に

⃗

p(t) = a⃗

u

11 t

+ b⃗

u

2

· (−1)

t

.

2

_⃗

_{p(0) =}

(

₁ 2 1 2

)

から

a, b

を定めると

a = 1, b = 0

で

p(t) = ⃗

⃗

u

1

.

極限分

布は

⃗

p(∞) = ⃗u

1 であり定常分布

.

3

_⃗

_{p(0) =}

(

₁ 3 2 3

)

から

a, b

を定めると

⃗

p(t) = ⃗

u

1

+

1 ₆

⃗

u

2

· (−1)

t

.

振動する

ので極限分布は存在しない

.

この場合

, x = 1, 2

は周期的状態である

.

(10)

マルコフ連鎖の時間発展の数値計算マルコフ連鎖の時間発展の数値計算

ここまで来たよ

3 マルコフ連鎖の時間発展

4 マルコフ連鎖の時間発展の数値計算

マルコフ連鎖の時間発展の数値計算

状態数が大きく規則的なマルコフ連鎖の時間発展の数値計算

マルコフ連鎖での母期待値の数値計算

(11)

マルコフ連鎖の時間発展の数値計算 I

状態

x = 0, . . . , m

− 1

の

m

状態のマルコフ連鎖を考える

.

C

で次のように書こう

.

状態番号が

1

ずれる…

分布

⃗

p(t), p(x, t)

→

1

d o u b l e p [m] =

{ 1 . 0 , 0 . 0 , . . . . , 0 . 0 } ; /∗配列 . m は整数 . ∗/

2

/

∗ {p ( 0 , t ) , p ( 1 , t ) , p ( 2 , t ) , . . . p (m−1, t )} ∗/

推移確率行列

M = (

p

11

p

12

p

21

p

22

)

→

1

d o u b l e M [ ] [ m]=

{ { 0 . 1 , 0 . 3 } ,

2

{ 0 . 9 , 0 . 7 } } ; /∗ 2 次元配列 ∗/

{{p

11

, p

12

}, {p

21

, p

22

}}

行列とベクトルの積

⃗

q = M ⃗

p

→ q

y

=

∑

x

M

yx

p

x

.

(12)

時間発展の計算の疑似コード

1 p [ ] を p ( x , 0 ) で初期化 ; 2 p を出力 ; 3 f o r ( t ){ 4 pn=M p ; /∗行列とベクトルの積∗/ 5 p=pn ; 6 p を出力 ; 7 }

ソースコード 1: マルコフ連鎖の時間発展

1 /∗ 2 M a r k o v 連鎖の分布の時間発展の計算

3 Time−stamp : ”2016−05−23 Mon 06:28 JST hig ”

4 ∗/

5 #d e f i n e CRT SECURE NO WARNINGS /∗ Visual C++ に必要なおまじない ∗/

6 #i n c l u d e < s t d i o . h> 7 8 /∗ 状態数 m ∗/ 9 #d e f i n e NS 3 10 11 i n t m u l t i p l y t r a n s ( d o u b l e ∗pn , double ∗p ) ; 12 i n t p r i n t d i s t ( d o u b l e ∗p , i n t t , i n t m) ; 13 14 i n t main ( ){ 15 i n t t , tmax ; 16 i n t x ; 17 d o u b l e p [ NS ] ; /∗ 分布 p ( t ) ∗/ 18 d o u b l e pn [ NS ] ; /∗ 分布 p ( t +1) ∗/ 19 i n t m=NS ; /∗ 状態数 ∗/ 20 21 22 s c a n f ( ”%d ” , &tmax ) ; p r i n t f ( ”#T=%d\n” , tmax ) ;

(13)

マルコフ連鎖の時間発展の数値計算マルコフ連鎖の時間発展の数値計算ソースコード 2: マルコフ連鎖の時間発展 28 連鎖の分布の時間発展の計算に必要なおまじない状態数分布分布状態数初期分布 29 30 f o r ( t =1; t<=tmax ; t ++){ 31 m u l t i p l y t r a n s ( pn , p ) ; /∗ 漸化式で分布を更新 ∗/ 32 f o r ( x =0; x<m; x++){ 33 p [ x ]= pn [ x ] ; 34 } 35 p r i n t d i s t ( p , t ,m) ; 36 } 37 r e t u r n 0 ; 38 } 39 40 /∗ p にMをかけて q=M p を書き込む. ∗/ 41 i n t m u l t i p l y t r a n s ( d o u b l e ∗q , double ∗p ){ 42 i n t x , y ; 43 i n t m=NS ; 44 /∗ 遷移確率行列 ∗/ 45 d o u b l e M [ ] [ NS ] ={{0.5 , 0 . 5 , 0 . 0 } , 46 { 0 . 5 , 0 . 5 , 0 . 0 } , 47 { 0 . 0 , 0 . 0 , 1 . 0 } } ; 48 f o r ( y =0; y<m; y++){ 49 q [ y ] = 0 ; 50 f o r ( x =0; x<m; x++){ 51 q [ y]+=M[ y ] [ x ]∗ p [ x ] ; 52 } 53 } 54 r e t u r n 1 ; 55 } 56 57 /∗ tとpを1行に出力∗/ 58 i n t p r i n t d i s t ( d o u b l e ∗p , i n t t , i n t m){ 59 i n t x ; 60 p r i n t f ( ”%d , ” , t ) ; 61 f o r ( x =0; x<m; x++){ 62 p r i n t f ( ”%f , ” , p [ x ] ) ; 63 } 64 p r i n t f ( ”\n” ) ; 65 r e t u r n 0 ; 66 }

(14)

マルコフ連鎖の時間発展の数値計算状態数が大きく規則的なマルコフ連鎖の時間発展の数値計算

ここまで来たよ

3 マルコフ連鎖の時間発展

4 マルコフ連鎖の時間発展の数値計算

マルコフ連鎖の時間発展の数値計算

状態数が大きく規則的なマルコフ連鎖の時間発展の数値計算

マルコフ連鎖での母期待値の数値計算

(15)

L07-Q1

Quiz(

ランダムウォークの時間発展

)

次の推移確率行列を持つ, 状態空間

_{{x} = {0, 1, 2, . . . , 99} 上のマルコフ連鎖を考}

えよう.

M =







7

10

2

10

0 · · · ·

0

3

10

5

10

2

0 ..

.

0 ₁₀

3 ₁₀

5 . .. ... ...

..

.

0 . .. ...

₁₀

2

0 ..

.

. ..

₁₀

3 ₁₀

5 ₁₀

2

0 · · · ·

0 ₁₀

3 ₁₀

8 





.

1

d o u b l e p [ 1 0 0 ] , q [ 1 0 0 ] ;

で表される ⃗p, ⃗q に対して, 入力 p を受け取り ⃗q = M ⃗p を計算する関数

1

i n t

m u l t i p l y t r a n s ( d o u b l e q [ ] ,

d o u b l e p [ ] ) ;

を書こう. 行列 M を 2 次元配列で表現せず, M の規則性を利用して書くこと.

(16)

これ長さ 100 の区間のランダムウォーク

数学的には長さ

_∞

の区間のランダムウォークも考えられる

_→

無限状態

マルコフ連鎖

計算機上ではとりあえず有限状態にしないとやってられない

.

ウォーカー

が端から外に出ようとしたらどうする

?

境界条件

(

一般用語

)

壁

(

ランダ

ムウォーク用語

)

吸収壁

0

に到達したウォーカーはそれ以降動かない

反射壁

0

に到達したウォーカーは

₋₁

には行かないけど

, +1

に戻る

ことはある

周期境界条件

0

から

−1

に行こうとしたら

m

− 1 = 99

に飛ぶ

ワー

プ壁

上の例はどれ

?

他の壁に改造すると

?

(17)

マルコフ連鎖の時間発展の数値計算マルコフ連鎖での母期待値の数値計算

ここまで来たよ

3 マルコフ連鎖の時間発展

4 マルコフ連鎖の時間発展の数値計算

マルコフ連鎖の時間発展の数値計算

状態数が大きく規則的なマルコフ連鎖の時間発展の数値計算

マルコフ連鎖での母期待値の数値計算

(18)

マルコフ連鎖での母期待値の数値計算

定義

p(x, t) = P (X(t) = x)

から

,

マルコフ連鎖での母期待値

E[ϕ(X(t))] =

∑

x

ϕ(x)f (x)

=

m

∑

x=1

ϕ(x)P (X(t) = x)

=

m

∑

x=1

ϕ(x)p(x, t)

(19)

L07-Q2

Quiz(マルコフ連鎖の母期待値の時間発展)

次の推移確率行列を持つ

,

状態空間

_{{x} = {1, 2}}

上のマルコフ連鎖を考

えよう

.

M =

(

₁

6

1

3

5

6

2

3

)

.

初期分布を

⃗

p(0) = (

1

0

)

とする

1

母期待値

_E[(X(t))

2

_]

を求めよう

_.

2

定常分布を

_⃗

_p(

∞)

とするとき

_{, log}

|⃗p(t) − ⃗p(∞)|

を求めよう

_.

この関

数が初期条件によってどう変化するか考えよう

.

(20)

お知らせ

月昼樋口オフィスアワー

(1-502)

チューター

/Math

ラウンジ月火水木昼

1-614

プチテスト

2016-05-30

月

4

▶

確率統計及び演習 I の外部記憶ペーパーをまとめに使えば?

https://register.math.ryukoku.ac.jp/archive/

統計検定勉強会今回は受験しない人も歓迎

2016-05-26

木

http://www.math.ryukoku.ac.jp/toukei-kentei/

https://manaba.ryukoku.ac.jp

マイページの下の方に

manaba

出席カード

提出