iT邦幫忙

1

【隨機模型】馬可夫鏈(Markov Chain)簡介 #附numpy程式實作

馬可夫鏈是為狀態空間中經過從一個狀態到另一個狀態的轉換的隨機過程,
並且下一個狀態的機率分佈只能跟前一個狀態有關

例子

譬如說有A, B, C三個城鎮,
三個城鎮的居民預都可能移動到其它城鎮,
假設每年的人口遷移狀況為:

  • 原本住A鎮的,有0.5的人口會移動至B鎮,有0.5的人口會移動至C鎮
  • 原本住B鎮的,100%會移動至C鎮
  • 原本住C鎮的,100%會移動至A鎮

那麼,我們可以畫一張圖來表示人口移動的狀況:
https://ithelp.ithome.com.tw/upload/images/20200705/20117114GUYeS54Ra2.png

我們發現,
當(A, B, C)的人口數為(0.4,0.2,0.4)的時候,
下一年(A, B, C)的人口數還是(0.4,0.2,0.4),
這時便形成了一個穩定態(stationary distribution)

穩定態的定理

有一個定理說:
馬可夫鏈滿足finite, irreducible, aperiodic,
則它有唯一的穩定態

  1. finite就是說狀態數是有限的,比如說有A, B, C三個城鎮,是有限多個
  2. irreducible是說每個狀態都可以在若干步後抵達其它每個狀態(包括自己)
    譬如說A城鎮可以去B, C兩個城鎮,
    可以走A->C->A就回到自己城鎮了
    B城鎮可以去C城鎮,
    要去A城鎮就先去C再去A

每個城鎮都可以去到其它城鎮,滿足irreducible

  1. aperiodic比較複雜,
    aperiodic是「periodic」的反義詞,
    periodic是說有某個state只能走d的倍數(d>1)那麼多步回到自己,
    譬如說:
    https://ithelp.ithome.com.tw/upload/images/20200705/20117114wQiVrhg6rx.png

狀態0要回到自己,只能走3的倍數那麼多步回到自己,那就是「periodic」

上一個例子的A城鎮可以走兩步回到自己城鎮(A->C->A),
也可以走三步(A->B->C->A),
那就是「aperiodic」

numpy程式求穩定態

本文的目標是學如果馬可夫鏈有穩定態,
用程式把它求出來

我們先把上例的A,B,C改成0,1,2,方便寫數學式
https://ithelp.ithome.com.tw/upload/images/20200705/201171142rBFXj3e4X.png

先用手算說明狀態如何轉移,
假設第一年(0,1,2)的人口數為[a,b,c],
則第二年的人口分布為[c*1, a*0.5, b*1+a*0.5]

馬可夫鏈的狀態轉移,
我們可以簡單用矩陣來寫,
定義轉移矩陣T, T的「第i列第j行的元素」即為狀態i轉至狀態j的比例

譬如本題的轉移矩陣為:
https://ithelp.ithome.com.tw/upload/images/20200705/20117114yHgAAJsi4p.png

要由前一個狀態[a,b,c]得到下一個狀態,只要做矩陣乘法就可以:
https://ithelp.ithome.com.tw/upload/images/20200705/20117114WleyZ8IlwF.png

程式碼實作

如果我們已經從理論上知道馬可夫鏈有唯一的穩定態,
那一開始的「起始狀態」隨便設都可以,
反正都會收斂

import numpy as np

# 現在狀態(Current state)
I = np.array([[1, 0, 0]])

# 轉移矩陣(Transition Matrix)
T = np.array([[0, 0.5, 0.5],
              [0, 0, 1],
              [1, 0, 0]])
              
# 迭代夠多次自然會到穩定態
for i in range(60):
    I = np.dot(I, T)
    print(I)

看懂的話應該很簡單吧?
只要「初始狀態」、「轉移矩陣」設定好,
跑for迴圈一直做乘法就好

<擴充功能>
如果你想要找穩定態最接近的有理數(也許方便寫學校作業/images/emoticon/emoticon39.gif),
可以用Fraction模組

import numpy as np
from fractions import Fraction

# 現在狀態(Current state)
I = np.array([[1, 0, 0]])

# 轉移矩陣(Transition Matrix)
T = np.array([[0, 0.5, 0.5],
              [0, 0, 1],
              [1, 0, 0]])

# 迭代夠多次自然會到穩定態
for i in range(60):
    I = np.dot(I, T)
    print(I)
    # 找最接近的有理數解
    for ele in I:
        for e in ele:
            print(Fraction(e).limit_denominator(), end=' ')
    print()

尚未有邦友留言

立即登入留言