[譯] 用 Python 實現馬爾可夫鏈的初級教程

• 原文地址：Markov Chains in Python: Beginner Tutorial
• 原文作者：Sejal Jaiswal
• 譯文出自：掘金翻譯計劃
• 本文永久連結：github.com/xitu/gold-m…
• 譯者：cdpath
• 校對者：jianboy, Kasheem Lew

學習馬爾可夫鏈及其性質，瞭解轉移矩陣，並用 Python 動手實現！

馬爾可夫鏈是通常用一組隨機變數定義的數學系統，可以根據具體的概率規則進行狀態轉移。轉移的集合滿足馬爾可夫性質，也就是說，轉移到任一特定狀態的概率只取決於當前狀態和所用時間，而與其之前的狀態序列無關。馬爾可夫鏈的這個獨特性質就是無記憶性。

跟隨本教程學會使用馬爾可夫鏈，你就會懂得離散時間馬爾可夫鏈是什麼。你還會學習構建（離散時間）馬爾可夫鏈模型所需的元件及其常見特性。接著學習用 Python 及其 numpy 和 random 庫來實現一個簡單的模型。還要學習用多種方式來表示馬爾可夫鏈，比如狀態圖和轉移矩陣。

想用 Python 處理更多統計問題？瞭解一下 DataCamp 的 Python 統計學思維課程！

開始吧……

為什麼要用馬爾可夫鏈？

馬爾可夫鏈在數學中有廣泛使用。同時也在經濟學，博弈論，通訊原理，遺傳學和金融學領域有廣泛應用。通常出現在統計學，尤其是貝葉斯統計，和資訊理論上下文中。在現實中，馬爾可夫鏈為研究機動車輛的巡航定速系統，抵達機場的乘客的排隊序列，貨幣匯率等問題提供瞭解決思路。最早由 Google 搜尋引擎提出的 PageRank 就是基於馬爾可夫過程的演算法。Reddit 有個叫子版塊模擬器的子版塊，帖子和評論全部用馬爾可夫鏈自動生成生成，厲害吧！

馬爾可夫鏈

馬爾可夫鏈是具有馬爾可夫性質的隨機過程。隨機過程或者說具有隨機性質是指由一組隨機變數定義的數學物件。馬爾可夫鏈要麼有離散狀態空間（一組隨機變數的可能值的集合）要麼有離散索引集合（通常表示時間），鑑於此，馬爾可夫鏈有眾多變種。而通常所說的「馬爾可夫鏈」是指具有離散時間集合的過程，也就是離散時間馬爾可夫鏈（DTMC）。

離散時間馬爾可夫鏈

離散時間馬爾可夫鏈所包含的系統的每一步都處於某個狀態，步驟之間的狀態隨機變化。這些步驟常被比作時間的各個瞬間（不過你也可以想成物理距離或者隨便什麼離散度量）。離散時間馬爾可夫鏈是隨機變數 X1，X2，X3 … 的序列，不過要滿足馬爾可夫性質，所以轉移到下一概率只和現在的狀態有關，與之前的狀態無關。用概率數學公式表示如下：

Pr( Xn+1 = x | X1 = x1, X2 = x2, …, Xn = xn) = Pr( Xn+1 = x | Xn = xn)

可見 Xn+1 的概率只和之前的 Xn 的概率有關。所以只需要知道上一個狀態就可以確定現在狀態的概率分佈，滿足條件獨立（也就是說：只需要知道現在狀態就可以確定下一個狀態）。

Xi 的可能取值構成的可數集合 S 稱為馬爾可夫鏈狀態空間。狀態空間可以是任何東西：字母，數字，籃球比分或者天氣情況。雖說時間引數通常是離散的，離散時間馬爾可夫鏈的狀態空間卻沒有什麼廣泛採用的約束條件，還不如參考任意狀態空間下的過程。不過許多馬爾可夫鏈的應用都用到了統計分析更簡單的有限或可數無窮狀態空間。

模型

馬爾可夫鏈用概率自動機表示（相信我它沒有聽上去那麼複雜！）。系統狀態的改變叫做轉移。各個狀態改變的概率叫做轉移概率。概率自動機包括從已知轉移到轉移方程的概率，將其轉換為轉移矩陣。

還可以將馬爾可夫鏈看作有向圖，其中圖 n 的邊標註的是 n 時刻狀態轉移到 n+1 時刻狀態的概率，Pr(Xn+1 = x | Xn = xn)。這個式子可以讀做，從已知狀態 Xn 轉移到狀態 Xn+1 的概率。這個概念也可以用從時刻 n 到時刻 n+1 的轉移矩陣來表示。狀態空間的每個狀態第一次出現是作為轉移矩陣的行，第二次是列。矩陣的每個元素都表示從這一行表示的狀態轉移到列狀態的概率。

如果馬爾可夫鏈有 N 種狀態，轉移矩陣就是 N x N 維，其中（I, J）表示從狀態 I 轉移到狀態 J 的概率。此外，轉移矩陣一定是概率矩陣，也就是每一行元素之和一定是 1。為什麼？因為每一行表示自身的概率分佈。

所以模型的主要特徵包括：狀態空間，描述了特定轉移發生的概率的轉移矩陣以及由初始分佈給出的狀態空間的初始狀態。

好像很複雜？

我們來看一個簡單的例子幫助理解這些概念：

如果 Cj 難得心情不好，她會跑步，或者大吃特吃冰淇淋（譯者注：原文 gooble 應為 gobble），要麼打個盹兒來調整。

根據以往資料，如果她睡了一覺調整心情，第二天她有 60% 的可能去跑步，20% 的可能繼續待在床上，還有 20% 的可能吃一大份冰淇淋。

如果她跑步散心，第二天她有 60% 的可能接著跑步，30% 的可能吃冰淇淋，只有 10% 的可能會去睡覺。

最後，如果她難過時縱情冰淇淋，第二天只有 10% 的可能性繼續吃冰淇淋，有 70% 的可能性跑步，還有 20% 的可能性睡覺。

上面由狀態圖表示的馬爾可夫鏈有 3 個可能狀態：睡覺，跑步和冰淇淋。所以轉移矩陣是 3 x 3 矩陣。注意，離開某一狀態的箭頭的值的和一定是 1，這跟狀態矩陣每一行元素之和是 1 一樣，都表示概率的分佈。轉移矩陣中每個元素的含義跟狀態圖的每個狀態類似。

這個例子應該會幫助你理解與馬爾可夫鏈有關的幾個不同概念。不過在現實世界中如何應用這一理論呢？

藉助這個例子，你應該能夠回答這種問題：「從睡覺狀態開始，2 天后 Cj 最後選擇跑步（跑步狀態）的概率是多少？」

我們一起算一下。要從睡覺狀態轉移到跑步狀態，Cj 有如下選擇：第一天繼續睡覺，第二天跑步（0.2 ⋅ 0.6）；第一天換成跑步，第二天繼續跑步（0.6 ⋅ 0.6）；第一天去吃冰淇淋，第二天換成跑步（0.2 ⋅ 0.7）。算下來概率是：((0.2 ⋅ 0.6) + (0.6 ⋅ 0.6) + (0.2 ⋅ 0.7)) = 0.62。所以說，從睡覺狀態開始，2天后 Cj 處於跑步狀態的概率是 62%。

希望這個例子可以告訴你馬爾可夫鏈網路都可以解決哪些問題。

同時，還可以更好地理解馬爾可夫鏈的幾個重要性質：

• 互通性：如果一個馬爾可夫鏈可以從任何狀態轉移至任何狀態，那麼它就是不可還原的。換句話說，如果任兩個狀態之間存在一系列步驟的概率為正，就是不可還原的。
• 週期性：如果馬爾可夫鏈只有在大於 1 的某個整數的倍數時返回某狀態，那麼馬爾可夫鏈的狀態是週期性的。因此，從狀態「i」開始，只有經過整數倍個週期「k」才能回到「i」，k 是所有滿足條件的整數的最大值。如果 k = 1 狀態「i」不是週期性的，如果 k > 1，「i」才是週期性的。
• 瞬態性和常返性：如果從狀態「i」開始，有可能無法回到狀態「i」，那麼狀態「i」有瞬態性。否則具有常返性（或者說持續性）。如果某狀態可以在有限步內重現，該狀態具有常返性，否則沒有常返性。
• 遍歷性：狀態「i」如果滿足非週期性和正重現性，它就有遍歷性。如果不具有可還原性的馬爾可夫鏈的每個狀態都有遍歷性，那麼這個馬爾可夫鏈也具有遍歷性。
• 吸收態：如果無法從狀態「i」轉移到其他狀態，「i」處於吸收態。因此，如果 當 i ≠ j 時，pii = 1 且 pij = 0，狀態「i」處於吸收態。如果馬爾可夫鏈的每個狀態都可以達到吸收態，稱其為具有吸收態的馬爾可夫鏈。

竅門：可以看看這個網站給出的馬爾可夫鏈的視覺化解釋。

用 Python 實現馬爾可夫鏈

我們用 Python 來實現一下上面這個例子。當然實際使用的庫實現的馬爾可夫鏈的效率會高得多，這裡還是給出例項程式碼幫助你入門……

先 import 用到的庫。

import numpy as np
import random as rm
複製程式碼

然後定義狀態及其概率，也就是轉移矩陣。要記得，因為有三個狀態，矩陣是 3 X 3 維的。此外還要定義轉移路徑，也可以用矩陣表示。

# 狀態空間
states = ["Sleep","Icecream","Run"]

# 可能的事件序列
transitionName = [["SS","SR","SI"],["RS","RR","RI"],["IS","IR","II"]]

# 概率矩陣（轉移矩陣）
transitionMatrix = [[0.2,0.6,0.2],[0.1,0.6,0.3],[0.2,0.7,0.1]]
複製程式碼

別忘了，要保證概率之和是 1。另外在寫程式碼時多列印一些錯誤資訊沒什麼不好的！

if sum(transitionMatrix[0])+sum(transitionMatrix[1])+sum(transitionMatrix[1]) != 3:
    print("Somewhere, something went wrong. Transition matrix, perhaps?")
else: print("All is gonna be okay, you should move on!! ;)")
複製程式碼

All is gonna be okay, you should move on!! ;)
複製程式碼

現在就要進入正題了。我們要用 numpy.random.choice 從可能的轉移集合選出隨機樣本。程式碼中大部分引數的含義從引數名就能看出來，不過引數 p 可能比較費解。它是可選引數，可以傳入樣品集的概率分佈，這裡傳入的是轉移矩陣。

# 實現了可以預測狀態的馬爾可夫模型的函式。
def activity_forecast(days):
    # 選擇初始狀態
    activityToday = "Sleep"
    print("Start state: " + activityToday)
    # 應該記錄選擇的狀態序列。這裡現在只有初始狀態。
    activityList = [activityToday]
    i = 0
    # 計算 activityList 的概率
    prob = 1
    while i != days:
        if activityToday == "Sleep":
            change = np.random.choice(transitionName[0],replace=True,p=transitionMatrix[0])
            if change == "SS":
                prob = prob * 0.2
                activityList.append("Sleep")
                pass
            elif change == "SR":
                prob = prob * 0.6
                activityToday = "Run"
                activityList.append("Run")
            else:
                prob = prob * 0.2
                activityToday = "Icecream"
                activityList.append("Icecream")
        elif activityToday == "Run":
            change = np.random.choice(transitionName[1],replace=True,p=transitionMatrix[1])
            if change == "RR":
                prob = prob * 0.5
                activityList.append("Run")
                pass
            elif change == "RS":
                prob = prob * 0.2
                activityToday = "Sleep"
                activityList.append("Sleep")
            else:
                prob = prob * 0.3
                activityToday = "Icecream"
                activityList.append("Icecream")
        elif activityToday == "Icecream":
            change = np.random.choice(transitionName[2],replace=True,p=transitionMatrix[2])
            if change == "II":
                prob = prob * 0.1
                activityList.append("Icecream")
                pass
            elif change == "IS":
                prob = prob * 0.2
                activityToday = "Sleep"
                activityList.append("Sleep")
            else:
                prob = prob * 0.7
                activityToday = "Run"
                activityList.append("Run")
        i += 1  
    print("Possible states: " + str(activityList))
    print("End state after "+ str(days) + " days: " + activityToday)
    print("Probability of the possible sequence of states: " + str(prob))

# 預測 2 天后的可能狀態
activity_forecast(2)
複製程式碼

Start state: Sleep
Possible states: ['Sleep', 'Sleep', 'Run']
End state after 2 days: Run
Probability of the possible sequence of states: 0.12
複製程式碼

結果可以得到從睡覺狀態開始的可能轉移及其概率。進一步擴充這個函式，可以讓它從睡覺狀態開始，迭代上幾百次，就能得到終止於特定狀態的預期概率。下面改寫一下 activity_forecast 函式，加一些迴圈……

def activity_forecast(days):
    # 選擇初始狀態
    activityToday = "Sleep"
    activityList = [activityToday]
    i = 0
    prob = 1
    while i != days:
        if activityToday == "Sleep":
            change = np.random.choice(transitionName[0],replace=True,p=transitionMatrix[0])
            if change == "SS":
                prob = prob * 0.2
                activityList.append("Sleep")
                pass
            elif change == "SR":
                prob = prob * 0.6
                activityToday = "Run"
                activityList.append("Run")
            else:
                prob = prob * 0.2
                activityToday = "Icecream"
                activityList.append("Icecream")
        elif activityToday == "Run":
            change = np.random.choice(transitionName[1],replace=True,p=transitionMatrix[1])
            if change == "RR":
                prob = prob * 0.5
                activityList.append("Run")
                pass
            elif change == "RS":
                prob = prob * 0.2
                activityToday = "Sleep"
                activityList.append("Sleep")
            else:
                prob = prob * 0.3
                activityToday = "Icecream"
                activityList.append("Icecream")
        elif activityToday == "Icecream":
            change = np.random.choice(transitionName[2],replace=True,p=transitionMatrix[2])
            if change == "II":
                prob = prob * 0.1
                activityList.append("Icecream")
                pass
            elif change == "IS":
                prob = prob * 0.2
                activityToday = "Sleep"
                activityList.append("Sleep")
            else:
                prob = prob * 0.7
                activityToday = "Run"
                activityList.append("Run")
        i += 1    
    return activityList

# 記錄每次的 activityList
list_activity = []
count = 0

# `range` 從第一個引數開始數起，一直到第二個引數（不包含）
for iterations in range(1,10000):
        list_activity.append(activity_forecast(2))

# 檢視記錄到的所有 `activityList`    
#print(list_activity)

# 遍歷列表，得到所有最終狀態是跑步的 activityList
for smaller_list in list_activity:
    if(smaller_list[2] == "Run"):
        count += 1

# 計算從睡覺狀態開始到跑步狀態結束的概率
percentage = (count/10000) * 100
print("The probability of starting at state:'Sleep' and ending at state:'Run'= " + str(percentage) + "%")
複製程式碼

The probability of starting at state:'Sleep' and ending at state:'Run'= 62.419999999999995%
複製程式碼

那麼問題來了，計算得到的結果為何會趨於 62%？

注意 這實際是「大數定律」在發揮作用。大數定律是概率論定律，用來說明在試驗次數足夠多時，可能性相同的事件發生的頻率趨於一致。也就是說，隨著試驗次數的增加，實際比率會趨於理論或預測的概率。

馬爾可夫思維

馬爾可夫鏈教程就到此為止了。本文介紹了馬爾可夫鏈及其性質。簡單的馬爾可夫鏈是開始學習 Python 資料科學的必經之路。如果想要更多 Python 統計學資源，請參閱這個網站。

如果發現譯文存在錯誤或其他需要改進的地方，歡迎到 掘金翻譯計劃 對譯文進行修改並 PR，也可獲得相應獎勵積分。文章開頭的 本文永久連結 即為本文在 GitHub 上的 MarkDown 連結。

---

掘金翻譯計劃 是一個翻譯優質網際網路技術文章的社群，文章來源為 掘金 上的英文分享文章。內容覆蓋 Android、iOS、前端、後端、區塊鏈、產品、設計、人工智慧等領域，想要檢視更多優質譯文請持續關注 掘金翻譯計劃、官方微博、知乎專欄。