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

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

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

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

想用 Python 處理更多統計問題?瞭解一下 DataCamp 的 Python 統計學思惟課程android

開始吧……ios

爲何要用馬爾可夫鏈?

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

馬爾可夫鏈

馬爾可夫鏈是具備馬爾可夫性質的隨機過程。隨機過程或者說具備隨機性質是指由一組隨機變量定義的數學對象。馬爾可夫鏈要麼有離散狀態空間(一組隨機變量的可能值的集合)要麼有離散索引集合(一般表示時間),鑑於此,馬爾可夫鏈有衆多變種。而一般所說的「馬爾可夫鏈」是指具備離散時間集合的過程,也就是離散時間馬爾可夫鏈(DTMC)。github

離散時間馬爾可夫鏈

離散時間馬爾可夫鏈所包含的系統的每一步都處於某個狀態,步驟之間的狀態隨機變化。這些步驟常被比做時間的各個瞬間(不過你也能夠想成物理距離或者隨便什麼離散度量)。離散時間馬爾可夫鏈是隨機變量 X1,X2,X3 … 的序列,不過要知足馬爾可夫性質,因此轉移到下一律率只和如今的狀態有關,與以前的狀態無關。用機率數學公式表示以下:算法

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

可見 Xn+1 的機率只和以前的 Xn 的機率有關。因此只須要知道上一個狀態就能夠肯定如今狀態的機率分佈,知足條件獨立(也就是說:只須要知道如今狀態就能夠肯定下一個狀態)。bash

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 連接。


掘金翻譯計劃 是一個翻譯優質互聯網技術文章的社區,文章來源爲 掘金 上的英文分享文章。內容覆蓋 AndroidiOS前端後端區塊鏈產品設計人工智能等領域,想要查看更多優質譯文請持續關注 掘金翻譯計劃官方微博知乎專欄

相關文章
相關標籤/搜索