使用 Python 的马尔可夫链与模型示例

在线运行和编辑此教程中的代码

马尔可夫链是一种数学系统，通常定义为一组随机变量，它们按照某些概率规则从一个状态转移到另一个状态。这些转移集合满足马尔可夫性质：转移到任一特定状态的概率仅取决于当前状态和经过的时间，而不取决于之前状态的序列。由于这一独特特性，马尔可夫过程被称为无记忆。

在本教程中，您将了解何时可以使用马尔可夫链、什么是离散时间马尔可夫链（DTMC）。您还将学习构建（离散时间）马尔可夫链模型所需的组成部分及其一些常见性质。接下来，您将使用 Python 以及其 numpy 和 random 库实现一个简单模型。您还将学习用状态图和转移矩阵等方式来表示马尔可夫链。

想用 Python 探索更多统计主题？请查看 DataCamp 的《Python 中的统计思维》课程！

让我们开始转移……

为什么使用马尔可夫链？

马尔可夫链在数学中用途广泛。它们被广泛应用于经济学、博弈论、通信理论、遗传学和金融。在统计学，尤其是贝叶斯统计和信息论情境中也广泛出现。谈到现实问题，它们被用来提出研究方案，例如研究机动车巡航控制系统、机场排队系统、货币汇率等。为谷歌搜索引擎最初提出的 PageRank 算法就基于马尔可夫过程。Reddit 的 Subreddit Simulator 是一个全自动的子版块，使用马尔可夫链生成随机帖文和评论，超酷！

马尔可夫链

马尔可夫链是具有马尔可夫性质的随机过程。随机过程（也常称为随机性或随机过程）是定义为一组随机变量的数学对象。马尔可夫链要么具有离散的状态空间（随机变量可能取值的集合），要么具有离散的指标集（通常表示时间）——因此存在很多马尔可夫链的变体。通常，“马尔可夫链”一词保留用于时间集合为离散的过程，即离散时间马尔可夫链（DTMC）。

离散时间马尔可夫链

离散时间马尔可夫链涉及一个系统：在每一步处于某个状态，状态在各步之间随机变化。这些步通常被看作时间点（当然也可以指物理距离或任何其他离散度量）。离散时间马尔可夫链是具有马尔可夫性质的随机变量序列 X₁、X₂、X₃、…，其性质是：转移到下一状态的概率仅取决于当前状态，而不取决于之前的状态。用数学概率公式表示为：

Pr( X_n+1 = x | X₁ = x₁, X₂ = x₂, …, X_n = x_n) = Pr( X_n+1 = x | X_n = x_n)

如您所见，X_n+1 的概率只依赖于它之前的 X_n 的概率。这意味着，仅需知道上一状态就足以确定当前状态的概率分布，满足条件独立的规则（换句话说：您只需要知道当前状态即可确定下一个状态）。

X_i 的可能取值构成一个可数集合 S，称为链的状态空间。状态空间可以是任何东西：字母、数字、篮球比分或天气状况。虽然时间参数通常是离散的，但离散时间马尔可夫链的状态空间并没有广泛一致的限制，而是指在任意状态空间上的过程。然而，许多马尔可夫链的应用采用有限或可数无限的状态空间，因为它们更易于进行统计分析。

模型

马尔可夫链可用概率自动机来表示（听起来复杂，其实并不！）。系统状态的变化称为转移。与各种状态变化相关联的概率称为转移概率。概率自动机会把给定转移的概率包含在转移函数中，从而形成转移矩阵。

您可以把它看作一系列有向图，其中第 n 个图的边用从时间 n 的某一状态到时间 n+1 的其他状态的概率标注，即 Pr(X_n+1 = x | X_n = x_n)。这可以读作：在给定状态 X_n 的值时转到状态 X_n+1 的概率。同样的信息也可以用从时间 n 到时间 n+1 的转移矩阵来表示。状态空间中的每个状态作为一行又作为一列出现一次，矩阵中的每个单元格告诉您从该行状态转移到该列状态的概率。

如果马尔可夫链有 N 个可能状态，则矩阵为 N x N 矩阵，其中条目 (I, J) 表示从状态 I 转移到状态 J 的概率。此外，转移矩阵必须是随机矩阵，即每一行的各项之和必须恰好为 1。为什么？因为每一行代表其自身的概率分布。

因此，该模型由状态空间、描述特定转移概率的转移矩阵，以及跨状态空间的初始状态（由初始分布给出）来表征。

一大串术语，对吧？

让我们看一个简单的示例来理解这些概念：

当 Cj 伤心时（这并不常见）：她要么去跑步，要么大吃冰淇淋，要么小睡一会儿。

根据历史数据，如果她在难过的一天选择睡觉消磨时光，第二天她去跑步的概率为 60%，继续宅在床上的概率为 20%，大吃冰淇淋的概率为 20%。

当她伤心去跑步时，第二天继续跑步的概率为 60%，大吃冰淇淋为 30%，而睡觉只有 10% 的概率。

最后，当她在难过的一天沉迷于冰淇淋时，第二天继续吃冰淇淋的概率仅为 10%，去跑步的概率为 70%，而睡觉的概率为 20%。

在该状态图中展示的马尔可夫链有 3 个可能状态：sleep、run、icecream。因此，转移矩阵将是一个 3 x 3 的矩阵。请注意，从某一状态发出的箭头之和总是恰好为 1，同样，转移矩阵中每一行的各项之和也必须恰好为 1——表示一个概率分布。在转移矩阵中，各单元格与状态图中的箭头起到相同的作用。

现在您已经看过这个例子，应该对与马尔可夫链相关的不同概念有了直观认识。但是，这些理论在现实生活中如何以及在哪里可以使用呢？

结合上述示例，您现在可以回答类似这样的问题：“从状态：sleep 出发，在一个难过的 2 天周期结束时，Cj 处于跑步（状态：run）的概率是多少？”

我们来算一下：为了从状态：sleep 转到状态：run，Cj 可能在第一次（或第一天）保持在状态：sleep，然后在第二次（第二天）转到状态：run（0.2 $\cdot$ 0.6）；或者第一天转到状态：run 并在第二天仍停留在那里（0.6 $\cdot$ 0.6）；或者她在第一次转到状态：icecream，然后在第二次转到状态：run（0.2 $\cdot$ 0.7）。因此概率为：((0.2 $\cdot$ 0.6) + (0.6 $\cdot$ 0.6) + (0.2 $\cdot$ 0.7)) = 0.62。也就是说，如果她从状态：sleep 开始，在难过的两天后转到状态：run 的概率为 62%。

希望这能让您了解使用马尔可夫链网络可以回答的各种问题。

同时，在这一点弄清楚之后，也更容易理解马尔可夫链的一些重要性质：

可约性：如果从任一状态到达任一状态都是可能的，则称马尔可夫链是不可约的。换言之，如果任意两个状态之间存在具有正概率的一系列步，则该马尔可夫链是不可约的。
周期性：如果链只能在某个大于 1 的整数倍数时返回到某一状态，则马尔可夫链中的一个状态是周期的。因此，从状态 “i” 出发，链只能在周期 “k” 的倍数时返回到 “i”，且 k 是这样的最大整数。如果 k = 1，则状态 “i” 为非周期的；如果 k > 1，则为周期的。
暂留性与常返性：若从状态 “i” 出发，存在非零概率永不返回 “i”，则称状态 “i” 为暂态。若不是暂态，则状态 i 为常返（或持久）态。若期望在有限步数内返回，则常返态称为正常返，否则称为零常返。
遍历性：若状态 “i” 既非周期又为正常返，则称其为遍历态。如果一个不可约马尔可夫链中的所有状态都是遍历态，则称该链为遍历链。
吸收态：如果不可能离开某一状态 i，则称其为吸收态。因此，当且仅当 p_ii = 1 且当 i ≠ j 时 p_ij = 0，状态 “i” 为吸收态。如果每个状态都能到达某个吸收态，则该马尔可夫链为吸收马尔可夫链。

提示：如果您还想看马尔可夫链的可视化讲解，请访问这个页面。

Python 中的马尔可夫链

让我们尝试用 Python 为上面的示例编写代码。尽管在实际应用中，您可能会使用一个更高效封装马尔可夫链的库，但这些代码能帮助您入门……

先导入您将用到的一些库。

import numpy as np
import random as rm
技术支持

现在定义状态及其概率：即转移矩阵。请记住，矩阵将是 3 X 3，因为您有三个状态。同时，您需要定义转移路径，也可以用矩阵来实现。

# The statespace
states = ["Sleep","Icecream","Run"]

# Possible sequences of events
transitionName = [["SS","SR","SI"],["RS","RR","RI"],["IS","IR","II"]]

# Probabilities matrix (transition matrix)
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 可能不是。它是一个可选参数，用于为采样集合指定概率分布，在这里就是转移矩阵。

# A function that implements the Markov model to forecast the state/mood.
def activity_forecast(days):
    # Choose the starting state
    activityToday = "Sleep"
    print("Start state: " + activityToday)
    # Shall store the sequence of states taken. So, this only has the starting state for now.
    activityList = [activityToday]
    i = 0
    # To calculate the probability of the 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))

# Function that forecasts the possible state for the next 2 days
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
技术支持

您会得到一组可能的随机转移及其发生概率，起始于状态：Sleep。进一步扩展该程序，例如在相同起始状态下迭代数百次，您就能看到在任一特定状态结束的期望概率及其估计值。让我们重写 activity_forecast 函数，并新增一些循环来完成这件事……

def activity_forecast(days):
    # Choose the starting state
    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

# To save every activityList
list_activity = []
count = 0

# `Range` starts from the first count up until but excluding the last count
for iterations in range(1,10000):
        list_activity.append(activity_forecast(2))

# Check out all the `activityList` we collected    
#print(list_activity)

# Iterate through the list to get a count of all activities ending in state:'Run'
for smaller_list in list_activity:
    if(smaller_list[2] == "Run"):
        count += 1

# Calculate the probability of starting from state:'Sleep' and ending at state:'Run'
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 中学习统计的资源，请务必查看此页面。

您是否有兴趣用 Python 探索更多统计学的实践案例？看看 DataCamp 的《统计思维案例研究》或《Python 网络分析入门》课程。

主题

Python 数据科学

主题

技术

Python 中的马尔可夫链：新手教程

使用 AI 探索